基于自适应EMD-NDFT的太阳光谱多普勒红移测速方法

王珍妮, 康志伟*, 刘 劲, 张 杰

1. 湖南大学信息科学与工程学院, 湖南 长沙 410082

2. 武汉科技大学信息科学与工程学院, 湖北 武汉 430081

引 言

太阳光谱蕴含着丰富的信息, 大规模的光谱巡天计划获取的太阳光谱数以亿计, 这对以太阳光谱作为信息源的天文测速导航技术发展有着重要意义。 利用太阳光谱进行速度测量的方法可以追溯到1960年, Franklin等提出通过天体光谱多普勒频移确定航天器速度的构想[1]: 在探测器上搭载一台摄谱仪, 实时获取天体的光谱信息, 解算出探测器相对于天体的径向速度。 2000年, Yim等提出用航天器的分光计测量多普勒频移, 再根据光谱学建模计算径向速度[2]。 2013年, 张伟等提出通过测量光谱红移来计算径向速度矢量的方法[3], 这为基于太阳光谱红移的天文测速导航提供了新思路。

太阳光谱多普勒红移测速是天文测速导航的重要环节, 它以太阳光作为导航源, 通过计算接收太阳光谱相对于标准太阳光谱的多普勒红移反推出航天器和太阳之间的相对径向速度。 由于太阳黑子、 日冕、 耀斑等太阳活动引发的光谱畸变造成了太阳光谱的不稳定, 这将影响着太阳光谱测速和导航的精度。 为此, 2017年, Ning等利用相邻两个观测周期的多普勒径向速度差作为滤波测量信息以减弱频谱失真的影响[4]。 2020年, 许俊等提出采用基于Sage-Husa噪声估计的自适应滤波方法抑制光谱畸变[5]。 2014年, 王永等运用小波变换消除光谱中的混合噪声以获取特征谱线, 通过密度估计法计算光谱红移和径向速度[6]。 另外, 针对深空探测任务具有时延长、 测量精度低、 信号微弱等特点[7], 2021年, Zhang等提出基于镜像NDFT-CS的快速、 高精度太阳光谱测速方法[8], 该方法通过对太阳光谱的NDFT的镜像低频进行稀疏处理, 能有效提高太阳光谱测速精度, 降低计算时间。

为了抑制太阳光谱测速导航中因光谱畸变对光谱测速造成的影响, 提高测量精度和计算速度, 本文提出一种自适应太阳光谱红移测速方法。 该方法旨在通过EMD-NDFT对光谱信号的时域处理与频域处理进行优势组合, 构建一种能对太阳光谱进行分解、 自适应滤波、 特征谱选择、 相关匹配以实现实时、 高精度的太阳光谱多普勒红移测速。

1 原理与方法

1.1 太阳光谱红移测速原理

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

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

(1)

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

(2)

(3)

光谱频移为

(4)

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

(5)

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

(6)

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

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

(7)

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

(8)

v=c(eτ-1)

(9)

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

1.2 EMD-NDFT算法

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

图1 自适应EMD-NDFT的多普勒速度估计框图

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

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

1.2.1 EMD滤波降噪过程

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

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

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

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

(10)

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

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

filter[imf1,imf2, …,imfm, …,imfn,r]

(11)

Step3: 重构信号。

(12)

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

1.2.1.1 经验模态分解

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

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

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

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

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

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

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

(13)

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

(14)

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

1.2.1.2 自适应阈值滤波

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

图2 EMD阈值滤波框图

具体步骤如下:

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

一阶IMF信号的能量

(15)

k阶IMF信号能量

(16)

(17)

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

Step2: 获取无噪信号。

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

Step3: 获取噪声IMF并滤波。

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

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

(18)

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

(19)

1.2.2 NDFT算法

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

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

(20)

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

(21)

其离散傅里叶变换为

(22)

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

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

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

(23)

(24)

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

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

(25)

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

(26)

(27)

2 实验及结果讨论

2.1 仿真条件及数据

从美国国家海洋与大气管理局(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所示为标准太阳光谱的傅里叶变换结果。

表1 太阳光谱数据信息

表2 仿真条件

图3 不同年份太阳黑子数量

图4 2018年太阳光谱

图5 标准太阳光谱

图6 傅里叶变换幅值

2.2 仿真过程

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

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

图8 噪声能量理论模型与实际信号能量对比

图9 理论模型能量与实际信号能量差值

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

表3 不同滤波方式下重构信号归一化方差

表4 不同滤波方式下速度误差

图10 一阶IMF

图11 待测光谱和其自适应阈值滤波重构

2.3 仿真结果及分析

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

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

图12 不同噪声下测速精度

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

表5 2018年不同波段速度精度对比(m·s-1)

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

表6 2015年不同波段速度精度对比(m·s-1)

2.4 EMD-NDFT算法的复杂度分析

为了验证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 程序运行时间对比

3 结 论

自适应EMD-NDFT太阳光谱多普勒红移测速方法是一种面向天文测速导航的信号处理方法。 该方法利用EMD算法将待测光谱进行了分解、 自适应滤波、 重构, 可消除航天器接收太阳光谱时的噪声; 通过NDFT算法将非均匀采样的太阳光谱转到频域, 并采用低频特征谱线进行泰勒匹配, 可计算出待测太阳光谱和标准太阳光谱之间的多普勒红移速度; 该方法将EMD的时域自适应滤波和NDFT的频域特征分解相结合, 在提高太阳光谱红移测速精度的同时, 降低了计算复杂度。 该方法针对不同年份和不同波长的光谱测速都有效果, 能适用于各类太阳光谱信号, 这对太阳光谱多普勒红移测速在天文自主测速导航中的应用具有探索意义。