一种天文时序信号的噪声估计方法与流程

本发明涉及一种天文时序信号的噪声估计方法，属于天文技术和信号处理技术领域。

背景技术：

几十年来，在天文观测时序信号中存在着各种的准周期振荡频率分量。研究和发现这些准周期振荡频率分量对于研究宇宙星系的结构演化有着至关重要的科学意义和价值。然而，由于观测信号中噪声的影响使得这些准周期振荡信号难于探测和提取。观测信号中的噪声主要以红噪声和白噪声为主。红噪声通常显示为不稳定的、非周期性的亮度波动。红噪声的功率谱服从幂律分布，即随着频率的增加而功率呈现衰减状态，而白噪声功率谱服从均匀分布，即功率谱是恒定的。这些噪声是天文观测信号中的固有性质，与观测设备和测量误差无关。因此，为了保证准周期振荡频率分量表征和提取的准确性必须对天文观测信号中的噪声水平进行估计。

通常的方式就是对时序信号的功率谱中的噪声进行建模，并获取模型中参数的数值，对信号的功率谱进行最优拟合，并从统计学的角度并给出信号功率谱噪声的置信水平。目前在对天文观测信号进行噪声估计时主要用的技术马尔科夫链蒙特卡罗(mcmc)方法；mcmc方法所采用的算法是metropolis–hastings(mh)算法以及改进的算法。然而，该传统算法在进行参数估计时都存在算法收敛速率慢和噪声参数估计准确性差的问题。

技术实现要素：

本发明提供了一种天文时序信号的噪声估计方法，以用于对天文时序信号的噪声进行估计。

本发明的技术方案是：一种天文时序信号的噪声估计方法，首先获取天文时序信号，对该时序信号进行归一化处理；接着，将归一化后的信号和窗函数相乘并运用周期图法得到该信号功率谱d；然后根据功率谱d定义该信号功率谱的噪声模型m和模型中的参数θ；定义该功率谱d在模型m和参数θ下的似然函数，以及参数θ的先验分布函数，根据贝叶斯公式定义参数θ在功率谱d和模型m下的后验概率分布函数；运用汉密尔顿蒙特卡洛算法来对后验概率分布函数进行采样，来获取模型m中参数θ的值，获得对该功率谱d的噪声模型m的最佳拟合曲线方程；对最佳拟合曲线方程进行卡方检验，计算置信水平。

所述方法具体步骤如下：

步骤1：获取天文时序信号，并利用该信号的均值和标准差对该时序信号进行归一化处理；

步骤2：将归一化后的信号和海明窗函数相乘并运用周期图法得到该信号功率谱d；

步骤3：首先根据功率谱d设计该功率谱的噪声模型m和模型参数θ；定义该功率谱d在模型m和参数θ下的似然函数p(d|θ,m)，以及参数θ的先验分布函数p(θ|m)，根据贝叶斯公式定义参数θ在功率谱d和模型m下的后验概率分布函数p(θ|d,m)＝p(d|θ,m)p(θ|m)；

步骤4：根据后验概率分布函数p(θ|d,m)，利用汉密尔顿蒙特卡洛算法对后验概率分布进行采样，来获取模型m中参数θ的值，获得对该功率谱d的噪声模型m的最佳拟合曲线方程，即噪声功率谱；

步骤5：对最佳拟合曲线方程进行卡方检验，计算99％的置信水平：将位于99％置信水平之上的频率成分认为是真正的振荡信号，在99％置信水平之下的认为是噪声。

本发明的有益效果是：本发明用均值和标准差来对信号进行归一化，可以移除信号直流分量对信号功率谱的影响；选择海明(hamming)窗来有效的抑制在频率泄露并很准确的表征信号的功率谱，可以避免由于天文观测信号中的振荡频率分量都是非平稳的不足；同时由于对信号进行了归一化进而保证了在定义模型参数的先验分布均值和方差的唯一性；而基于前述所有的结合，在结合使用汉密尔顿蒙特卡洛算法后，从而提高了对天文时序信号噪声估计的准确性和收敛速率；最后通过用99％的置信水平来保证噪声估计和振荡频率分量表征的准确性。

附图说明

图1是本发明天文观测信号中噪声估计方法的总体流程图；

图2是本发明中合成信号的归一化图；

图3是本发明中合成信号的功率谱图；

图4是本发明中mh算法的最佳拟合图和参数估计结果；

图5是本发明中hmc算法的最佳拟合图和参数估计结果；

图6是本发明中mh算法获取最佳拟合所需要收敛率变化图；

图7是本发明中hmc算法获取最佳拟合所需要收敛率变化图；

图8是本发明中hmc算法所获得的最佳拟合的95％置信水平图；

图9是本发明中hmc算法所获得的最佳拟合的99％置信水平图；

图10是本发明中天文观测中太阳耀斑的归一化时序信号；

图11是本发明中太阳耀斑的归一化时序信号的功率谱图；

图12是本发明中天文观测中太阳耀斑信号功率谱最佳拟合曲线以及它的99％置信水平。

具体实施方式

下面结合附图和实施例，对本发明作进一步说明，但本发明的内容并不限于所述范围。

实施例1：如图1所示，一种天文时序信号的噪声估计方法，首先获取天文时序信号，对该时序信号进行归一化处理；接着，将归一化后的信号和窗函数相乘并运用周期图法得到该信号功率谱d；然后根据功率谱d定义该信号功率谱的噪声模型m和模型中的参数θ；定义该功率谱d在模型m和参数θ下的似然函数，以及参数θ的先验分布函数，根据贝叶斯公式定义参数θ在功率谱d和模型m下的后验概率分布函数；运用汉密尔顿蒙特卡洛算法来对后验概率分布函数进行采样，来获取模型m中参数θ的值，获得对该功率谱d的噪声模型m的最佳拟合曲线方程；对最佳拟合曲线方程进行卡方检验，计算置信水平。

进一步地，所述方法的具体步骤如下：

步骤1：时序信号归一化：由于很难从天文实际的观测数据来评价本发明的方法和传统方法在噪声参数拟合精度的优劣，所以我们使用合成信号来展示和证明本发明的优越性。我们产生了一个功率谱的指数为2的红噪声(即红噪声功率谱为prednoise(f)＝af^-α，其中α＝2)合成信号，并利用信号的均值和标准差对其进行归一化：

其中，i为合成信号的幅度值，imean是均值，istd是标准差，inor是归一化之后的信号。结果如图2所示。此步骤可以有效的抑制信号的直流分量对信号噪声模型参数拟合精度的影响，同时通过此归一化步骤能够在步骤3时获取统一的噪声模型参数分布的均值和方差；

步骤2：获取信号的功率谱：为了抑制在计算信号功率谱时频率泄露，我们首先选择海明窗和归一化后的时序信号进行相乘，然后运用周期图法来获取信号的功率谱d；结果如图3所示。

步骤3：定义概率分布函数p(θ|d,m)：首先根据功率谱d设计该功率谱的噪声模型m和模型参数θ；定义该功率谱d在模型m和参数θ下的似然函数p(d|θ,m)，以及参数θ的先验分布p(θ|m)，根据贝叶斯公式来计算参数θ的后验概率分布p(θ|d,m)＝p(d|θ,m)p(θ|m)。根据图3所示的功率谱，表明该时序信号的功率谱存在红噪声。由于红噪声服从幂律分布，即p1(f)＝af^-α，则该信号功率谱的噪声模型m为：

m(f)＝p1(f)＝af^-α

其中m为模型，θ为模型参数，它包含了2个参数即a，α。由于在步骤1中已经对信号进行了归一化，因此定义参数α的先验分布为的正态分布，而参数a的先验分布为独立的对数正态分布：

任何长度为n的时序的功率谱似然函数p(d|θ,m)定义为：

d(fj)和m(fj)是在频率fj处的信号的功率谱和噪声模型m的功率谱值。则参数的后验概率分布为：

p(θ|d,m)＝p(d|θ,m)p(θ|m)

步骤4：汉密尔顿蒙特卡洛算法采样：根据后验概率分布函数p(θ|d,m)，利用汉密尔顿蒙特卡洛算法对后验概率分布进行采样，来获取模型m中参数θ的值，获得对该功率谱d的噪声模型m的最佳拟合曲线方程。

hmc算法区别于传统mcmc采样方法之处在于，它是运用汉密尔顿动力学理论从目标分布中的抽样，其中引入了虚拟的动量变量，在汉密尔顿系统下利用蛙跳(leapfrog)技术得到马氏链下一时刻的备选状态，然后依据metropolis准则判断是否接受汉密尔顿动力系统下得到的此备选状态，从而完成马尔科夫链的更新迭代。它减少了传统mcmc算法的随机游走行为，改进了马尔科夫链的有效性。

汉密尔顿系统中的每一个状态都由一个位置矢量x和一个随机动量矢量z组成，它们的联合密度为p(x,z)＝f(x)g(z)。目标就是从这个联合概率密度里面抽取样本，其中能量函数被定义为：

其中u(x)＝-logf(x)代表势能，代表着动能。联合概率密度p(x,z)∝exp{-e(x,z)}。汉密尔顿蒙特卡洛算法的采样过程如下：首先根据一个给定的随机位置变量x中生成动量分量z，然后根据汉密尔顿方程：

是否将状态(x,z)更新为下一个状态(x′,z′)的接受概率为：

pacc＝min(1,exp{e(x,z)-e(x′,z′)})

接受x′为下一状态值的概率为pacc，拒绝x′而接受x为下一状态值的概率为1-pacc，重复此过程，直到产生足够的样本数。

同时汉密尔顿蒙特卡洛算法利用蛙跳(leapfrog)技术，使得新状态(x′,z′)即使距离前序状态(x,z)很远也拥有很好的采样接受概率。蛙跳技术是一种被用作就解决高维状态下求解后验分布函数的积分的一种方法。蛙跳技术的运算迭代步骤是：首先，根据初始状态迭代得到1/2步长的位置变量其次，根据新的位置变量计算整个步长ε处的动量变量zi(t+ε)；最后，依据得到新的动量变量zi(t+ε)以及位置变量可以更新迭代剩下的1/2步长的动量变量，从而得到xi(t+ε)。具体公式如下：

图4是采用传统的mcmc方法对合成信号的拟合曲线以及对参数α的估计结果，图5则是本发明所采用的方法的拟合曲线和对参数α的估计结果。合成信号功率谱中α的真值等于2，从图可以看到目前传统方法mcmc所采用的mh算法的最优拟合和参数α的拟合结果(即图4)为1.9，而本发明的方法的最优拟合和参数α的拟合结果(图5)为2.01。从精度来说，本发明的方法的参数估计结果优于传统方法。其中，传统的方法获得的功率谱采用最接近本发明步骤的方式处理获得(采用利用均值的归一化及汉宁窗)，当功率谱采用和本发明步骤的方式处理获得时参数α的拟合结果优于1.9，但是本发明的结果2.01仍然更优。

图6、7显示了传统方法和本发明方法收敛速率变化的图。横轴表示参数估计方法达到收敛所需要的执行次数，纵轴表示收敛因子，当收敛因子接近于1时表示采样过程达到收敛。图6和图7对比来看，传统的mcmc方法大概需要3万次的迭代才能达到收敛，即收敛因子接近于1，而本发明的方法仅需要3000次左右即达到收敛。收敛的执行次数比传统方法快了一个量级。通过图4-图7的对比实验证明本发明比传统方法在对噪声进行估计时，可以有效的提高参数拟合精度和收敛速率。

步骤5：噪声估计：对图5得到的最佳拟合曲线进行卡方检验，并计算出99％的置信水平。由于该功率谱服从卡方分布因此用卡方检验来计算噪声模型m(f)的置信水平。定义(1-∈)的置信界限为γ∈，其中∈＝0.01，然后通过卡方概率密度来计算γ∈：

图8展示了95％置信水平、图9展示了99％置信水平。高于置信水平的功率谱的峰则被认定为信号中的振荡频率分量。从图8中的箭头所示的位置，可以看出有峰高于95％的置信水平。如果用95％置信水平来判断信号中的频率分量，那么箭头所指向的位置存在频率信号。但该信号是红噪声并不存在任何的频率分量，因而尽管采用了前面的步骤1-4提高参数估计的准确性，但由于选择了较小的置信水平而导致噪声被低估。我们从图9中可以看到并没有任何峰高于99％的置信水平，这意味着该信号中不含有频率分量，是噪声。图8、图9进一步证明了本发明通过创造性劳动获得的99％置信水平能有效提高时序信号噪声估计的准确性，存在显著的进步。

实施例2：如图1所示，一种天文时序信号的噪声估计方法，首先获取天文时序信号，对该时序信号进行归一化处理；接着，将归一化后的信号和窗函数相乘并运用周期图法得到该信号功率谱d；然后根据功率谱d定义该信号功率谱的噪声模型m和模型中的参数θ；定义该功率谱d在模型m和参数θ下的似然函数，以及参数θ的先验分布函数，根据贝叶斯公式定义参数θ在功率谱d和模型m下的后验概率分布函数；运用汉密尔顿蒙特卡洛算法来对后验概率分布函数进行采样，来获取模型m中参数θ的值，获得对该功率谱d的噪声模型m的最佳拟合曲线方程；对最佳拟合曲线方程进行卡方检验，计算置信水平。

进一步地，所述方法的具体步骤如下：

步骤1：时序信号归一化：获取天文时序信号，并利用该信号的均值和标准差对该时序信号进行归一化处理，如图10所示；

步骤2：获取信号的功率谱：为了抑制在计算信号功率谱时频率泄露，我们首先选择海明窗和归一化后的时序信号进行相乘，然后运用周期图法来获取信号的功率谱d；功率谱图如11所示。

步骤3：定义概率分布函数p(θ|d,m)：首先根据功率谱d设计该功率谱的噪声模型m和模型参数θ；定义该功率谱d在模型m和参数θ下的似然函数p(d|θ,m)，以及参数θ的先验分布p(θ|m)，根据贝叶斯公式来计算参数θ的后验概率分布p(θ|d,m)＝p(d|θ,m)p(θ|m)。根据图11所示的功率谱，该时序信号的对数功率谱在低频部分存在红噪声，在高频部分则存在白噪声。由于红噪声服从幂律分布，即p1(f)＝af^-α，而白噪声服从均匀分布，即p2(f)＝c，则该信号的噪声模型为：

m(f)＝p1(f)+p2(f)＝af^-α+c

其中a为大于0的常数，α为大于0的指数，c为尾部高频白噪声，f为频率，p为噪声的功率谱。这个模型的参数θ有三个：a，α，c。由于在步骤1中已经对信号进行了归一化，因此定义参数α的先验分布为的正态分布，而参数a和参数c的先验分布为独立的对数正态分布：和关于参数a，α的先验分布定义同实例1的定义。

步骤4：汉密尔顿蒙特卡洛算法采样：根据后验概率分布函数p(θ|d,m)，利用汉密尔顿蒙特卡洛算法对后验概率分布进行采样，来获取模型m中参数θ的值，获得对该功率谱d的噪声模型m的最佳拟合曲线方程(即噪声功率谱)。功率谱的最佳拟合曲线如图12虚线所示；

步骤5：噪声功率谱的水平估计：对最佳拟合曲线方程进行卡方检验，计算99％的置信水平。在99％置信水平之下的功率谱认为是信号噪声(此耀斑信号中包含红噪声和白噪声)的功率谱，而高于99％的功率谱则是该信号所蕴含的频率分量。功率谱的99％置信水平如图12点划线所示；在图12中用箭头指示了该耀斑信号中存在的振荡频率分量的位置。

由于传统的噪声估计方法在选择置信水平时是选择95％的置信水平，但这个阈值会存在对噪声低估的问题，因而，在本发明中将此置信水平定义为99％，即高于99％的功率谱才是该信号所蕴含的频率分量。这样有效的避免在对信号振荡频率分量进行表征时的误判。

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

起点商标作为专业知识产权交易平台，可以帮助大家解决很多问题，如果大家想要了解更多知产交易信息请点击 【在线咨询】或添加微信 【19522093243】与客服一对一沟通，为大家解决相关问题。

此文章来源于网络,如有侵权,请联系删除