光谱多普勒效应是指深空探测器在远离(或接近)导航光源的过程中, 光的频率减小(或增加)的现象。 光谱频率变化反映了探测器与导航光源的相对运动。 当航天器远离太阳时, 使得观测到的谱线波长比静止的谱线波长要长, 表现为谱线朝红端移动一段距离, 叫作多普勒红移。 当航天器相对光源运动时, 会产生多普勒频移, 将频移转化为光谱红移, 就可推导出航天器的相对径向速度。 具体计算过程如下:

假设多普勒红移为z, 航天器相对于太阳的径向速度为v、 太阳光的原始频率和波长分别为f₀和λ ₀、 相对运动后观测到的太阳光频率和波长分别为f和λ , 光速为c, 那么多普勒红移公式为

z=λ-λ0λ0=f0-ff=vc(1)

所以在不考虑相对论的情况下, 航天器远离太阳时, 光谱频率和波长的关系式分别如式(2)和式(3)

f0=1+vcf(2)

λ=1+vcλ0(3)

光谱频移为

Δf=f0-f=vcf(4)

由式(4)可知光谱波长的移动为

Δλ=λ-λ0=vcλ0(5)

基于多普勒效应, 太阳光谱发生多普勒频移。 因此需要建立标准太阳光谱与实测太阳光谱之间的关系方程, 在太阳多普勒速度估计过程中, 以标准太阳光谱为模板, 与实测太阳光谱进行比较, 从而得到它们之间的多普勒频移。 因此, 对于相同的导航源, 由上述公式可知, 根据多普勒频移可以获得多普勒速度。 由式(5)可知无法直接建立光谱的波长的移动量Δ λ 和速度v之间的关系, 所以将其转换为对数形式, 即

log(λ)=log(λ0)+log1+vc(6)

假设模板光谱模型为s(λ ₀), 待测光谱模型为p(λ ₀), 根据式(5)可知

p(λ0)=s(λ0+Δλ)+nλ(λ0)(7)

式(7)中, p和s是模型函数, n_λ 是噪声。 所以式(6)可以表示为

p[log(λ0)]=slog(λ0)+log1+vc+nλ(λ0)(8)

令 τ=log1+vc, 因此只要计算出τ , 就可以得到速度v

v=c(eτ-1)(9)

由式(9)可知, 计算太阳光谱多普勒的波长位移, 并且将之转换为与速度之间的关系, 就可以得到最终结果。

EMD-NDFT算法结构如图1所示, 主要分为两个部分:

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 自适应EMD-NDFT的多普勒速度估计框图Fig.1 Doppler velocity estimation block diagram of adaptive EMD-NDFT

(1)自适应EMD阈值滤波降噪过程。 首先对待测太阳光谱数据进行预处理, 对信号进行经验模态分解, 再对分解的IMF信号根据其自身特点自适应降噪, 最后对降噪后的IMF信号重构得到重构太阳光谱。

(2)用NDFT的低频部分进行特征匹配求解多普勒速度。 对标准太阳光谱和重构光谱分别进行NDFT(非均匀傅里叶变换), 并提取频谱幅值的高能量的低频部分作为计算量, 再进行匹配得到多普勒速度。

1.2.1 EMD滤波降噪过程

经验模态分解(empirical mode decomposition, EMD)是一种数据驱动的时域分解算法, 可自适应地将信号按照频率和幅值大小依次分解成一组本征模态函数(intrinsic mode function, IMF)^[9], 被广泛地应用于非平稳、 非线性信号的去噪处理。

假设标准太阳光谱信号为x₀, 待测太阳光谱信号是x, 采用EMD自适应滤波降噪过程如下:

Step1: 对待测信号进行EMD分解。

[imf1, imf2, …, imfm, …, imfn, r]=emd(x)(10)

EMD将组成原始信号的各尺度分量从高频到低频进行提取, 分解得到的特征模态函数顺序是按频率由高到低进行排列的, 即首先得到最高频的分量, 然后是次高频的, 最终得到一个频率接近为0的残余分量r。 其中n为EMD分解后的层数, 根据信号特征, n的数量是自适应的。

Step2: 设置合适阈值, 并滤除噪声。 详细步骤将在1.2.1.2介绍。

[imf'1, imf'2, …, imf'm, …, imf'n, r]=filter[imf1, imf2, …, imfm, …, imfn, r](11)

Step3: 重构信号。

x'=R([imf'1, imf'2, …, imf'm, …, imf'n, r])(12)

式(12)中, R是重构函数, x'是重构后的光谱信号。

1.2.1.1 经验模态分解

假设太阳光谱信号为x, EMD方法处理光谱信号的具体过程如下:

Step1: 找到信号x所有的极值点;

Step2: 用3次样条曲线拟合出上下极值点的包络线E_max和E_min, 并求出上下包络线的平均值m, 在x中减去它: h=x-m;

Step3: 根据预设判据判断h是否为IMF;

Step4: 如果不是, 则以h代替x, 重复以上步骤直到h满足判据, 则h就是需要提取的IMF信号C_i;

Step5: 每得到一阶IMF, 就从原信号中减去它, 重复以上步骤, 直到信号最后剩余部分r_n就只是单调序列或者常值序列;

经过上述步骤将原始信号x分解成一系列IMF以及剩余部分的线性叠加

$\text{x}=\overset{N}{\mathop{\underset{i=1}{\mathop \sum }\, }}\, {{C}_{i}}+{{r}_{n}}$ (13)

此外, 余项的筛选标准^[10]如式(14)

${{h}_{q}}=mean\left\{ \overset{q}{\mathop{\underset{i=1}{\mathop \sum }\, }}\, \left\{ im{{f}_{i}}\left( t \right)-\text{mean}\left[ \frac{im{{f}_{i}}\left( t \right)}{\text{std}\left[ im{{f}_{i}}\left( t \right) \right]} \right] \right\} \right\}$ (14)

其中imf_i(t)为第i阶IMF, 当h_q偏离0时, 第q阶IMF即为混合模态与余项的区分边界。

1.2.1.2 自适应阈值滤波

本文基于有界噪声辅助分析^[11] 进行改进, 实现自适应EMD阈值滤波。 如图2所示。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 EMD阈值滤波框图Fig.2 EMD threshold filtering block diagram

具体步骤如下:

Step1: 得到EMD分解中噪声能量的理论传播模型^[10]。

一阶IMF信号的能量

${{E}_{H}}\left( 1 \right)=\overset{N}{\mathop{\underset{i=m}{\mathop \sum }\, }}\, imf_{i}^{2}\left( i \right)$ (15)

k阶IMF信号能量

${{E}_{H}}\left( k \right)=\frac{{{E}_{H}}\left( 1 \right)}{{{\beta }_{H}}}\rho _{H}^{-2\left( 1-H \right)k}$ (16)

ρH≈2.01+0.2h-12+0.12H-122(17)

其中N为EMD分解时窗的长度, H为赫斯特(Hurst)指数, 当H=0.5时, 对应噪声为高斯白噪声, β _H=0.719, ρ _H=2.01。

Step2: 获取无噪信号。

将式(15)、 式(16)建立的理论模型与含噪声信号的实际能量传播进行匹配, 若已知一阶IMF信号能量为E_H(1)根据式(16)就可以推算出余下各个阶次IMF信号的能量大小E_H(k)。 将其与实际IMF信号能量对比, 两者差值越大, 说明该信号携带的噪声信息越小。 当给两者取对数后, 差值偏离0很大时, 信号可定义为低频无噪声信号。 将太阳光谱数据作为EMD的输入, 若实际能量在估计能量95%的置信度区间内, 则为含噪信号。

Step3: 获取噪声IMF并滤波。

在一定置信度下与理论能量分布模型匹配的IMF定义为噪声IMF。

由于第一阶IMF信号可被视为噪声, 但是imf₁在某些点携带了极大的能量, 若直接将其视作噪声滤除, 会使重构信号丢失信息。 因此可以先对imf₁进行阈值滤波, 留下能量高的谱线。 滤波过程中采用随机取样imf₁与剩余的IMF构造具有相同信噪比的信号, 对应可以得到多个滤波结果, 对其取平均值即可得到最终滤波结果。 具体阐述如下:

(1)对imf₁信号进行均值滤波, 去除部分极少量高能量信号, 得到imf'₁。

(2)将imf'₁与其余IMF构造相同信噪比的信号: 由假设共有m-2个信号需要进行阈值滤波, 那么就需要构造匹配这些IMF信号的同信噪比的噪声信号, 作为该阶IMF的阈值, 进行滤除, 部分高阶信号可视为无噪声信号, 不作处理, 那么一次滤波后的重构信号表示为

${{x}_{1}}\left( t \right)=\overset{m-1}{\mathop{\underset{i=2}{\mathop \sum }\, }}\, imf{{'}_{i}}\left( t \right)+\overset{N}{\mathop{\underset{i=m}{\mathop \sum }\, }}\, im{{f}_{i}}\left( t \right)+r\left( t \right)$(18)

其中imf'₁(t)是各阶滤波后的本征模态信号。

(3)重复上述过程, 进行M次迭代, 最终滤波结果为

$\hat{x}\left( t \right)=\frac{1}{M}\overset{M}{\mathop{\underset{k=1}{\mathop \sum }\, }}\, {{x}_{k}}\left( t \right)$(19)

1.2.2 NDFT算法

快速傅里叶变换算法(FFT)是信号分析的重要工具, 为了得到多普勒红移速度与光谱波长之间的代数关系, 本文需对太阳光谱进行非均匀采样的傅里叶变换(NDFT), 具体计算过程如下:

傅里叶变换的公式如式(20)

X(jω)=∫0∞x(t)e-jωtdt(20)

那么非均匀采样的太阳光谱信号的傅里叶变换可以表示为

X(jω)=∫0∞x[log(λ)]e-jωlog(λ)d[log(λ)](21)

其离散傅里叶变换为

$\text{X}\left( \text{k} \right)=\overset{N-1}{\mathop{\underset{n=0}{\mathop \sum }\, }}\, x\left[ \text{log}\left( \lambda \right) \right]{{\text{e}}^{-j\frac{2\text{ }\!\!\pi\!\!\text{ }}{{{L}_{\lambda }}}k\text{log}\left( \lambda \right)}}\text{ }\!\!\Delta\!\!\text{ log}\left( \lambda \right)$(22)

其中, k=1, 2, …, N。

根据傅里叶变换的性质可知, 时域上的位移只会导致频域上的相移。 因此待测光谱信号的NDFT结果只与标准太阳光谱的相位不同, 振幅不会发生变化。 那么就可以通过标准光谱和待测光谱的非均匀傅里叶变换相位函数的互相关结果获得频移量, 从而得到多普勒红移, 进而得到径向速度。

x[log(λ)]↔X(k)(23)

$x\left[ \text{log}\left( \lambda \right)-\text{log}\left( {{\lambda }_{0}} \right) \right]\leftrightarrow {{\text{e}}^{-j\frac{2 \pi}{{{L}_{\lambda }}}k\text{log}\left( {{\lambda }_{0}} \right)}}X\left( k \right)$(24)

1.2.3 基于自适应EMD-NDFT的多普勒测速

首先按照1.2.1的步骤对待测太阳光谱进行经验模态分解、 自适应阈值滤波、 重构, 得到重构的待测光谱信号。 结合式(8)和式(24)可知待测光谱和标准光谱的非均匀傅里叶变换关系式如式(25)

$\begin{array}{{l}} {{P}_{k}}=\overset{N}{\mathop{\underset{n=1}{\mathop \sum }\, }}\, x\left[ \text{log}\left( {{\lambda }_{0}} \right)+\tau +{{n}_{\lambda }}\left( {{\lambda }_{0}} \right) \right]{{\text{e}}^{-j\frac{2 \pi}{{{L}_{\lambda }}}k\text{log}\left( {{\lambda }_{0}} \right)}}\text{ }\!\!\Delta\!\!\text{ log}\left( {{\lambda }_{0}} \right) \\ ={{N}_{k}}+b{{S}_{k}}{{\text{e}}^{-j\frac{2 \pi}{{{L}_{\lambda }}}k\tau }}\frac{k}{{{L}_{\lambda }}} \\ \end{array}$(25)

而待测光谱非均匀傅里叶变换也可以由式(22)表示, 那么就可以两者之间的关系求到τ , 由文献[13]可知, 两者的匹配关系为泰勒匹配, 关系式如式(26)和式(27)

$\text{Z}=\overset{N}{\mathop{\underset{k=1}{\mathop \sum }\, }}\, k{{P}_{k}}{{S}_{k}}\text{sin}\left( {{\phi }_{k}}-{{\theta }_{k}}+k\frac{2 \pi}{{{L}_{\lambda }}} \right)$(26)

$\left[ {{\hat{m}}} \right]=\text{min}\left( Z \right)$ (27)

${\hat{m}}$是Z的最小估计下标, 即为τ , 那么代入到式(9)就可以得到径向速度。

从美国国家海洋与大气管理局(National Oceanic and Atmospheric Administration, NOAA)网站^[12]上可以获得太阳黑子活动周期图, 如图3所示, 2018和2020年太阳黑子的数量很少, 而2012年— 2015年是黑子数量多的年份, 在2014年达到了顶峰。 由于要与文献[8]的精度作对比, 所以本文同样选择2018年的太阳光谱作为标准太阳光谱, 并计算多普勒速度估计误差的标准偏差。 太阳光谱数据来自于欧洲南方天文台(European Southern Observatory, ESO)网站^[13], 部分信息如表1所示。 本文中使用的标准太阳光谱如图4所示。 模拟条件如表2所示, 那么光谱图像如图5所示。 图6所示为标准太阳光谱的傅里叶变换结果。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 不同年份太阳黑子数量Fig.3 Number of sunspots in different years

表1

Table 1

表1(Table 1)

表1 太阳光谱数据信息 Table 1 Solar spectral data information 名称 数值 光谱范围 [gte: 378.193, lte: 691.301] 天文台 La Silla Paranal APEX 数据来源 ESO, PHASE3 信噪比/dB 323.4 观测日期 2018-04-08 T21:07:12.957

表1 太阳光谱数据信息 Table 1 Solar spectral data information

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 2018年太阳光谱Fig.4 Solar spectrum in 2018

表2

Table 2

表2(Table 2)

表2 仿真条件 Table 2 Simulation conditions Variable Value 波长/nm 5 381.94~5 431.94 光谱的波长分辨率/nm 0.1 噪声水平/(erg/s/cm2/angstrom) 1 000 速度真值/(m· s-1) 40

表2 仿真条件 Table 2 Simulation conditions

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 标准太阳光谱Fig.5 Standard Solar spectrum

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 傅里叶变换幅值Fig.6 Fourier transform amplitude

根据仿真条件对待测太阳光谱进行自适应EMD分解并重构。 图7所示为待测光谱EMD分解结果, 待测光谱被分解为12个IMF信号。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 待测太阳光谱EMD分解结果Fig.7 EMD results of solar spectrum to be measured

根据文献[10]可知: 在EMD分解中, 各级噪声能量的理论值与各层IMF信号的实际能量值相差越大, 则该信号的噪声信号越小。 在噪声能量理论传播模型95%置信区间内的IMF分量为含噪信号, 那么如图8所示, 第imf₂~imf₄可定义为含噪信号, imf₅~imf₁₂为可定义为低频信号, 即近似不含噪声信号。 图9所示为原信号能量与预测能量的差。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 噪声能量理论模型与实际信号能量对比Fig.8 Comparison between theoretical model of noise energy and actual signal energy

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 理论模型能量与实际信号能量差值Fig.9 Difference between theoretical model energy and actual signal energy

由于在自适应阈值降噪时, 是以imf₁为噪声模板进行阈值滤波的, 那么imf₁幅值大小对降噪效果具有很大的影响。 图10为一阶IMF信号图, 由于imf₁某些点的幅值极高, 若直接将imf₁作为噪声模板可能对结果有影响, 表3所示为不同滤波方式下重构信号的归一化方差, 表4所示为imf₁的处理方式对最终测速精度的影响。 由表3和表4可知, 先将一阶IMF先进行中值滤波后保留高能量部分, 然后把滤除的信号作为噪声模板, 再进行重构后的结果是最优的。 再结合框图2的步骤对待测光谱进行自适应阈值滤波重构结果如图11所示。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 一阶IMFFig.10 IMF1

表3

Table 3

表3(Table 3)

表3 不同滤波方式下重构信号归一化方差 Table 3 Normalized variance of reconstructed signal under different filtering methods 滤波操作 归一化方差 剔除一阶IMF 以imf1作为模板噪声 0.501 0 以imf'1作为模板噪声 0.463 1 保留一阶IMF 以imf1作为模板噪声 1.184× 10-4 以imf'1作为模板噪声 7.092× 10-5

表3 不同滤波方式下重构信号归一化方差 Table 3 Normalized variance of reconstructed signal under different filtering methods

表4

Table 4

表4(Table 4)

表4 不同滤波方式下速度误差 Table 4 Velocity error under different filtering modes 滤波操作 速度误差/ (m· s-1) 剔除一阶IMF 以imf1作为模板噪声 0.341 3 以imf'1作为模板噪声 0.211 3 保留一阶IMF 以imf1作为模板噪声 0.134 1 以imf'1作为模板噪声 0.064 9

表4 不同滤波方式下速度误差 Table 4 Velocity error under different filtering modes

图11

Fig.11

	Figure Option ViewDownloadNew Window
	图11 待测光谱和其自适应阈值滤波重构Fig.11 Spectrum to be measured and its adaptive threshold filtering reconstruction

为了验证EMD-NDFT算法的优越性, 本文做了以下几组实验:

(1)与NDFT-CS方法^[8]进行了光谱测速精度对比。 由文献[8]可知: 采用NDFT-CS方法对太阳光谱进行测速时, 测速精度与噪声通量水平成正比。 与NDFT-CS方法对比, 本文采用自适应EMD-NDFT方法, 使用相同的太阳光谱数据, 在同一实验条件下, 同样做了1 000次蒙特卡洛实验以求解太阳光谱的多普勒速度, 结果如图12所示。 可以看出: 噪声通量水平较小时, 两种方法的精度差距很小。 为了体现本文方法具有更高的抗噪性能, 我们加大了噪声能量, 在噪声通量水平10 000~100 000(erg/s/cm²/angstrom)下, 对比EMD-NDFT和NDFT-CS求解速度的误差精度, 图13为EMD-NDFT算法比NDFT-CS算法的测速精度的优胜率, 结合图12, 可以得到如下结论: 当噪声不断增大时, EMD-NDFT算法所得到的测速精度明显高于NDFT-CS所得结果。

图12

Fig.12

	Figure Option ViewDownloadNew Window
	图12 不同噪声下测速精度Fig.12 Velocity measurement accuracy under different noise

图13

Fig.13

	Figure Option ViewDownloadNew Window
	图13 EMD-NDFT比NDFT-CS的测速精度提升率Fig.13 Higher speed measurement accuracy of EMD-NDFT

(2)对不同光谱波段进行速度测量, 对比EMD-NDFT和NDFT-CS算法的测速精度。 如表5所示, 为噪声通量大小1 000 (erg/s/cm²/angstrom)时, 不同波段的速度误差结果。 经过多次实验对比都可以看出, 本文方法在不同波长范围内的测速精度更优。

表5

Table 5

表5(Table 5)

表5 2018年不同波段速度精度对比(m· s-1) Table 5 Comparison of velocity accuracy of different bands in 2018 (m· s-1) 波长范围/nm NDFT-CS EMD-NDFT 4 281.94~4 331.94 0.900 1 0.639 9 4 781.94~4 831.94 0.150 2 0.110 0 5 381.94~5 431.94 0.080 6 0.063 9 5 781.94~5 831.94 0.061 8 0.040 8

表5 2018年不同波段速度精度对比(m· s-1) Table 5 Comparison of velocity accuracy of different bands in 2018 (m· s-1)

(3)采用其他年份光谱进行测速, 对比两种方法的精度。 由于我们选择的是2018年的太阳光谱数据, 由图3可知, 2018年是太阳黑子较少的年份, 太阳活动较弱, 不能完全保证“ 本文算法精度更高” 这个论点的可靠性, 因此选择太阳黑子高发年份2015年的数据再次作对比实验, 结果如表6所示, 自适应EMD-NDFT方法的测速精度总是高于文献[8]的NDFT-CS算法。

表6

Table 6

表6(Table 6)

表6 2015年不同波段速度精度对比(m· s-1) Table 6 Comparison of velocity accuracy of different bands in 2015 (m· s-1) 波长范围/nm NDFT-CS EMD-NDFT 4 281.94~4 331.94 0.619 8 0.291 9 4 781.94~4 831.94 0.201 3 0.162 0 5 381.94~5 431.94 0.164 7 0.152 1 5 781.94~5 831.94 0.242 3 0.165 7

表6 2015年不同波段速度精度对比(m· s-1) Table 6 Comparison of velocity accuracy of different bands in 2015 (m· s-1)

为了验证EMD-NDFT算法的快速性, 本节从理论和实验两个方面对EMD-NDFT算法复杂度进行分析, 并与NDFT-CS算法进行比较。

2.4.1 算法原理复杂度分析

由于NDFT-CS算法需要使用标准太阳光谱构造测量矩阵Φ (L× N)、 然后再通过测量矩阵与傅里叶变换低频部分的幅值相乘得到零相位矢量Θ =Φ (L× N)· S(N× 1), 再将标准太阳光谱与待测光谱的傅里叶变换低频部分相乘得到观测矢量y=Φ (L× N)· P(N× 1), 然后再进行特征匹配获取多普勒径向速度, 其中N表示光谱的谱线数量, L表示从频域中选取的谱线数量。 在构造零相位矢量Θ 和观测矢量y的过程中, 虽然已经将信号进行了稀疏处理, 但是数据量仍然较大。

而EMD-NDFT算法不需要给分解后的每个IMF信号都进行阈值处理, 只需根据条件给少量噪声IMF进行阈值滤波, 而且无需进行额外的构造矩阵和多次乘积运算, 且仅采用非均匀傅里叶变换的低频部分计算速度时循环次数少, 计算量更小。

2.4.2 程序实际运行结果对比

表7所示为两种算法在光谱分辨率为0.1 nm时的程序运行时间对比, 由实际运行数据可知, 使用EMD-NDFT算法更快速得到结果。

表7

Table 7

表7(Table 7)

表7 程序运行时间对比 Table 7 Program running time comparison 方法 速度误差/(m· s-1) 运行时间/s NDFT-CS 0.080 6 1.631 4 EMD-NDFT 0.063 9 0.396 5

表7 程序运行时间对比 Table 7 Program running time comparison