#### 1. 项目目的:

恒星光谱库根据其来源可分为理论恒星光谱库和经验恒星光谱库, 二者各具优缺点具有很强的互补性, 是天文研究中不可或缺的工具库. 经验恒星光谱库主要的问题之一体现在其恒星在参数空间中的分布是无规则的随机的. 光谱库构建者考虑的是如何尽可能的将参数空间中的''空穴''填补上,却无法控制其参数分布的规则与否. 受限与这种无规则分布情况, 其使用过程中常常采用邻近插值的0阶近似方法, 而不是像理论库中可以采用更高阶的插值计算方法. 其计算精度往往比较粗糙. 这里我们建立了基于径向基神经网络的插值算法用于经验恒星光谱库的插值计算.

#### 2. 算法介绍:

##### 2.1 径向基神经网络结构:

径向基神经网络是有输入输出层和由径向基函数构成的隐藏层组成的三层的神经网络. 具体到这里的光谱插值计算问题时, 恒星的大气参数($\rm T_{eff}, log\,g, [Fe/H]$) 作为网络的输入量 $\mathbf{x}$, 恒星的光谱作为输出层的 $\mathbf{y}$, 隐藏层由核函数$\mathbf{K}$组成. 则有:

$$
y(\mathbf{x})=\sum_{j=1}^N K_j(\mathbf{x}-\mathbf{x}_j)\cdot c_j,
$$

权重系数$\mathbf{c}$可以通过求解线性方程组获得;

$$
\sum_{j=1}^N K_j(\mathbf{x}_i-\mathbf{x}_j)\cdot c_j = y_i \ (i = 1, 2, 3,..., N),
$$

##### 2.2 自适应核函数计算:

这里核函数采用高斯函数的形式:

$$
	K_j(\mathbf{x}_i-\mathbf{x}_j) = {\rm e}^{\frac{-\sum_{d=1}^D(x_{d,i}-x_{d,j})^2}{2\sigma_j^2}},
$$

核函数影响域$\sigma$在经典的径向基网络中通常是人为设定的常数. 但这一设定对于样本空间密度变化较大的情况是不恰当的, 我们借鉴光滑粒子流体算法中核函数的自适应方法引入了对$\sigma$的约束:

$$
\begin{equation}
1 \approx \sum_{j=1}^N \frac{c_0}{(\sqrt{2\pi})^D} \cdot {\rm e}^{\frac{-\sum_{d=1}^D(x_{d,i}-x_{d, j})^2}{2\cdot \sigma_i^2}}  (i = 1,2,3,...,N),
\end{equation}
$$

式中$c_0$作为人为设点的只有参数用来限制核函数影响域的整体大小. 解此方程, 其核函数的影响域在光谱库中恒星样本分布稀疏的地方取较大值, 而在稠密的趋于取较小值.

![MILES光谱库中样本恒星对应的核函数影响域的分布](./fig/isotropy_sigma.png)

上图中给出了MILES光谱库中样本恒星对应核函数影响域的分布, 可以很好的反映出影响域与样本密度的关系. 在确定了核函数之后就可以将其带入到线性方程组中进行计算,求解出系数矩阵$\rm c$ . 在确定$\rm c$之后, 以后的恒星光谱计算就可以通过$\sum_{j=1}^N K_j(\mathbf{x}-\mathbf{x}_j)\cdot c_j$计算得到. 此外对于核函数的各向异性处理以及最佳参数($c_0$以及相应的各向异性参数)的寻优过程可以参考[文献]( https://arxiv.corg/abs/2007.04634 "An upgraded interpolator of the radial basis function network for spectral calculation based on empirical stellar spectral library")中的介绍.

#### 3. 计算流程:

如上所述, 对于径向基插值计算,我们先需要通过求解线性方程组来计算出系数矩阵$\rm c$. 获得$\rm c$之后就可以直接进行插值计算. 需要注意的一点是, 插值计算中仅需要$\rm c$和$\sigma$ , 恒星光谱库的数据在完成初始的$\rm c$和$\sigma$就不在需要. 如此,算法对于某一特定的恒星光谱库进需要完成一次$\sigma$和$\rm c$矩阵的计算即可. 之后的插值计算只需输入插值恒星参数数组就可以批量的计算出其对应的恒星光谱(经验恒星光谱库插值光谱是归一化的, 这里选择对V-band流量进行归一化处理). 以下将以对MILES恒星光谱库的计算为例展示这一过程:

##### 3.1 插值计算案例

为使用方便, 原始的MILES光谱数据与恒星参数被转换并保存成hdf5文件, 对于已经设定好参数$c_0, c_1, c_2$,可以很方便的确定径向基网络:

~~~python
# -*- coding: utf-8 -*-
# @Author: chenglt
# @Date:   2020-01-14 15:26:30
# @Last Modified by:   chenglt
# @Last Modified time: 2020-04-05 12:10:48

import time
from tqdm import tqdm
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np
import sys
import func as rbf

def interpolator(para_input, para_lib, spec_lib, sigma):
    """
    Interpolated spectra calcultion for the input parameter.

Arguments:
        para_input {[narray]} -- [the parameters of stars]
        para_lib {[narray]} -- [the stellar parameters of stellar spctral library]
        spec_lib {[narray]} -- [the coorespendent spectra in library]
        sigma {[type]} -- [the paramter kernel Gaussian function ]
    """

para_lib_rescale, normal_factor = rbf.para_rescale(para_lib)
    para_rescale = para_input*normal_factor
    coeffi_matrix = rbf.coeffi_rbf(para_lib_rescale, sigma, spec_lib)
    spec_output = rbf.sepc_interp_out(para_rescale, para_lib_rescale, sigma, coeffi_matrix)

return spec_output

if __name__ == "__main__":
    para_input = np.loadtxt('./input.dat')   # the stellar parameters for interpolation

filterv = np.loadtxt('./johnsonUBV/GCPD-Johnson.V_ADPS.dat')
    flux_filterv = np.sum(0.5*(filterv[1:, 1]+filterv[:-1, 1])*(filterv[1:, 0]-filterv[:-1, 0]))

# two kind spectral lib here, by different maker 'MILES_mix' or 'MILES'
    miles_lib = pd.read_hdf('./lib/MILES_lib.h5', 'MILES')
    # when use miles original no [:-8] part.
    para_lib = miles_lib[:3].T.to_numpy()[:-8]
    spec_lib = miles_lib[3:].T.to_numpy()[:-8]
    wavelen_lib = np.array(list(map(float, miles_lib.index[3:].to_numpy())))
    # inside or outside the MILES lib stellar area: 
    ind = rbf.area_inside(para_lib, para_input)    
    ############################################################################################################
    # initial paramter for new method:
    c0_sph = 5.98  # this parameter, make the variance minimum in scripts: sigma_anisotropy0.py
    c1_sph = 6.43
    x_sph = 14.09/14.0914

############################################################################################################
    # sigma in new method:
    #sigma = np.load('./sigma.npy')
    sigma = rbf.sigma_rbf(para_lib, c0=c0_sph, c1=c1_sph, ratio=x_sph, anisotropy=1, fig=1)
    plt.close()
    np.save('sigma.npy', sigma)
    # interpolate the sepctra of the input parameters:
    spec_out = interpolator(para_input, para_lib, spec_lib, sigma)
    spec_nom = rbf.spec_normal(wavelen_lib, spec_out, filterv)   # out normalist flux by V-band flux of the spectra.

# show the interpolated spectra
    plt.plot(wavelen_lib, spec_nom[0]+3, 'k-')
    plt.text(5300, 6, r'$T_{\rm eff} = 20000\,{\rm K}, \ {\rm lg}\, g=4.2,\ [Fe/H]= 0.0.$', fontsize=7, style='oblique')

plt.plot(wavelen_lib, spec_nom[1]+2, 'k-')
    plt.text(5300, 5.5, r'$T_{\rm eff} = 10000\,{\rm K}, \ {\rm lg}\, g=4.3,\ [Fe/H]= 0.0.$', fontsize=7, style='oblique')

plt.plot(wavelen_lib, spec_nom[2]+1, 'k-')
    plt.text(5300, 5, r'$T_{\rm eff} = 5000\,{\rm K}, \ {\rm lg}\, g=4.5,\ [Fe/H]= 0.0.$', fontsize=7, style='oblique')

plt.plot(wavelen_lib, spec_nom[3], 'k-')
    plt.text(5300, 4.5, r'$T_{\rm eff} = 4360 \,{\rm K}, \ {\rm lg}\, g=2.0,\ [Fe/H]= 0.0$', fontsize=7, style='oblique')
    plt.xlabel(r'$\rm Wavelength/\AA $')
    plt.yticks([])
    plt.show()
    ############################################################################################################
    # unit replace for the output spectra
    BC_lib =pd.read_hdf('./lib/BC_lib.h5', 'BC');
    p, Rsun, Lsun, vsun, BCsun = 5.67E-5, 6.96e10, 3.83E33, 7.771327556479673e+09, -0.08;
    #unit(erg*s**(-1)*CM**(-2)*K**(-4)),(CM), (erg/s), (erg*s^(-1)*cm^(-2)).
    flux_v_sun = vsun*4*np.pi*Rsun**2
    BC_x = interpolate.griddata(BC_lib.to_numpy()[:, :3], BC_lib.to_numpy()[:, 3], para_input, method='nearest')
    Bol = p*para_input[:, 0]**4

flux_v = Bol*10**(0.4*(BC_x-BCsun)+np.log10(flux_v_sun/Lsun))
    spec_output = spec_nom*np.array(np.mat(flux_v/flux_filterv).T)  #erg*s^-1*cm^-2*A^-1
    
    ## output: in spec_out[:, 0], '0' means that the stellar outside the MILES library area. '1' means inside.  
    # The outside result need to be checked carefully
    mark = np.zeros(len(para_input))
    mark[ind] = 1
    np.savetxt('./spec_out.dat', np.vstack((np.array([-1]+list(wavelen_lib)), np.hstack((mark.reshape(-1, 1), spec_output)))))

~~~

经验恒星光谱库光谱原始的量纲意义不大, 在意行示例的最后一部分是对输出光谱进行量纲还原,最终的输出光谱量纲为$\rm erg\cdot s^{-1}\cdot cm^{-2}\cdot \AA^{-1}$ . 计算中还对输入参数是否属于外插进行了判断,对应于最终输出矩阵的第一列中的值("0"表示外插). 对于外插结果其可靠性较低.

##### 3.2 初始参数的确定

为了能够获得合适的初始参数$(c_0,c_1,c_2)$, 我们基于天牛须算法进行寻优搜索,  以上的计算中可以看出给定参数$(c_0, c_1,c_2)$后我们就能够唯一的确定插值网络. 参数$(c_0,c_1,c_2)$寻优计算的评价函数是通过比较插值计算结果与理论库BaSeL中原始光谱得到的, 这一点在[文献]( https://arxiv.corg/abs/2007.04634 "An upgraded interpolator of the radial basis function network for spectral calculation based on empirical stellar spectral library")中给出详细的介绍. 这里仍然是基于MILES来做演示:

~~~python
# -*- coding: utf-8 -*-
# @Author: chenglt
# @Date:   2019-12-25 10:11:07
# @Last Modified by:   chenglt
# @Last Modified time: 2020-01-05 18:23:00

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import interpolate
import sys
import func as rbf
import random

# import the lib data: for new library, a function will be used to creat the mock library.
para_basel = np.loadtxt('./lib/mock_lib/libpara.dat')[:, 1:]
spec_basel = np.loadtxt('./lib/mock_lib/libspec.dat')[:, 357:553]
wave_basel = np.loadtxt('./lib/mock_lib/libwave.dat')[357:553] * 10

filterb = np.loadtxt('./johnsonUBV/GCPD-Johnson.B_ADPS.dat')
filterv = np.loadtxt('./johnsonUBV/GCPD-Johnson.V_ADPS.dat')

miles_lib = pd.read_hdf('./lib/MILES_lib.h5', 'MILES')

para_lib = miles_lib[:3].T.to_numpy()
spec_lib = miles_lib[3:].T.to_numpy()
wavel_lib = np.array(list(map(float, miles_lib.index[3:].to_numpy())))

# important step
spec_basel = rbf.spec_normal(wave_basel, spec_basel, filterv)

index_inside = rbf.area_inside(para_lib, para_basel, fig=1)
num_inside = len(np.array(index_inside)[index_inside])

# get sigma in old work:
sigma_miles = np.zeros(len(para_lib))
para_miles_nor, _ = rbf.para_rescale(para_lib)
for i in range(len(para_lib)):
    dd = np.sum((para_miles_nor - para_miles_nor[i]) ** 2, 1) ** 0.5
    sigma_miles[i] = sum(np.sort(dd)[:4]) / 3

[c0, c1, x], p_value = rbf.bestc_BAS_3D_smooth(para_lib, spec_lib, para_basel[index_inside], spec_basel[index_inside],
                                           wavel_lib, wave_basel, filterv, num_extract=4528, steps_max=2000,
                                           step_factor_initial=0.5, step_factor_final=0.0001, eta=0.996)
np.savetxt('track.txt', p_value)   # save the trajectory
~~~

以上代码中,我通过天牛须算法取寻找最有的初始参数.可以人为设定最大迭代步数. 这一步的计算较为耗时,一般情况下可以得到一个收敛的结果, 将之带入3.1中就可以完成径向基网路的计算进而得到光谱库对应的插值器. 之后就可以以此插值器方便快捷的进行恒星光谱的计算.

outspacetime/Spectra interpolater for empirical stellar spectral lib

简介

发行版

贡献者

近期动态