一种基于维纳滤波建立综合脉冲星时的优化方法 |
|||||||
| 申请号 | CN202310267715.2 | 申请日 | 2023-03-19 | 公开(公告)号 | CN117290661A | 公开(公告)日 | 2023-12-26 |
| 申请人 | 中国科学院国家授时中心; | 发明人 | 张哲浩; 童明雷; 高玉平; 赵成仕; 朱幸芝; | ||||
| 摘要 | 本 发明 提供了一种基于维纳滤波建立综合脉冲星时的优化方法，参考时为TT(BIPM)，利用已有数据拟合模型得到精确模型参数；更改参考时为被监测参考时，得到拟合前计时残差并在设定时长内做平均，得到若干残差点；根据互 功率谱 上界计算并调整所有互功率谱，对确定互功率谱上界的两组残差使用其所计算出的互功率谱进行滤波，而对其余残差使用所有互功率谱的算术平均进行滤波；计算各组残差加入微小 信号 前后的相关系数，并以非相关系数的归一化值作为每颗星的权重；对每颗星的滤波结果进行加权，最终得到综合脉冲星时。本发明使得综合脉冲星时相对参考时间的 波动 量级处在正确的范围，对参考时间的改变具有良好的敏感度。 | ||||||
| 权利要求 | 1.一种基于维纳滤波建立综合脉冲星时的优化方法，其特征在于，包括以下步骤： |
||||||
| 说明书全文 | 一种基于维纳滤波建立综合脉冲星时的优化方法技术领域[0001] 本发明属于脉冲星时间尺度建立领域，涉及脉冲星计时模型建立，尤其是建立综合脉冲星时间尺度的方法。 背景技术[0002] 时间是当前测量精度最高的基本物理量。精确的时间计量，是许多科学与工程领域的基础，高精度时间基准目前主要依赖于国际原子时(TAI)。TAI是由国际计量局BIPM根据世界上守时实验室的数据计算得出，每月发布一次。而由TAI导出的地球时TT(BIPM)是目前基于原子频标实现的最精确时间尺度，国际计量局一般于次年初发布其与上一年TAI的钟差数据。对于有高精度时间需求的部门，数据依赖于国外发布使得高精度时间服务的可靠性降低，因此，有必要寻找独立的高精度时间保持方法。 [0003] 脉冲星是恒星演化末期形成的致密星体，具有自转与辐射特性，其中毫秒脉冲星是脉冲星自转加速后的演化产物，具有毫秒量级的自转周期与很高的自转稳定性。研究表明，一些毫秒脉冲星的自转稳定性在10年尺度上优于原子时。但直接观测得到的脉冲到达信号并不能作为时间信号使用，还需要依靠计时模型，即一系列描述脉冲星信号发射、传播时延以及时间改正的数学方程，以扣除各种物理过程对信号传播的影响。为了获得精确的模型参数，一般需要数年乃至10年以上的观测，而精确的计时模型一旦得到，就可以建立脉冲星时间。 [0004] 作为一种天文信号，脉冲星辐射的脉冲在传播中受到多种因素的影响，导致以此确定的脉冲到达时存在波动，并使得建立起的脉冲星时间存在噪声。脉冲星时的噪声以白噪声与低频噪声为主，其中白噪声水平目前要高于原子钟信号及其导出时间尺度。为了提高脉冲星时的稳定度，需要使用某些方法降低噪声水平，这一般是通过对多颗星的观测数据进行综合来实现。综合脉冲星时的基本原理是基于不同脉冲星噪声生成过程的不相关性，使用某些算法消去不相关噪声，从而得到参考某种时间尺度下的共性信号，再通过加权平均，进一步降低噪声水平。根据已有研究，降低了噪声的综合脉冲星时能够成功检测出TAI的波动，体现出其重要的应用价值。 [0005] 当前，建立综合脉冲星时的基本方法有维纳滤波法、贝叶斯法等。维纳滤波法的计算相对简单，易于快速实现，其原理是通过对不同脉冲星计时残差做相关运算，得到被检测信号统计特性的估计，再利用其对每颗脉冲星的残差进行滤波，最后对滤波结果做加权平均。维纳滤波法存在的问题是当多颗脉冲星中有几颗具有较高噪声水平时，建立的综合脉冲星时亦会存在很高的噪声，导致直接淹没被检测信号。改进信号能量估计方法，并合理取权，是利用维纳滤波法构建综合脉冲星时的关键。 发明内容[0006] 为了克服现有技术的不足，本发明提供一种基于维纳滤波建立综合脉冲星时的优化方法，根据低噪声脉冲星残差间的互相关对信号能量的估计，调整高噪声脉冲星残差互相关(互功率谱)量级，保证每组残差滤波结果的波动水平均不高于参考时信号的波动水平；同时，根据不同脉冲星残差对微小信号改变的敏感度确定权值，使综合脉冲星时能正确跟随参考时间的变化而变化。 [0007] 本发明解决其技术问题所采用的技术方案包括以下步骤： [0008] 1)利用初始模型，将参考时改为TT(BIPM)，利用已有数据拟合模型直至拟合前、后残差的均方根相等，得到精确的模型参数；更改参考时为被监测参考时，得到拟合前计时残差； [0009] 2)对第1)步得到的拟合前计时残差在连续的设定时长内做平均，得到若干个残差点； [0010] 3)对第2)步得到的残差点拟合并去除2次项； [0011] 4)估计互功率谱上界；根据互功率谱上界计算并调整所有互功率谱，使所有互功率谱在指定的低频区间内不超过所述上界；给出所有的互功率谱估计后，对确定互功率谱上界的两组残差，使用其所计算出的互功率谱进行滤波，而对其余残差，使用所有互功率谱的算术平均进行滤波； [0012] 5)向每组残差中加入一个微小信号，信号的幅度取为参考时信号波动的量级估计，计算各组残差加入微小信号前后的相关系数，并以非相关系数的归一化值作为每颗星的权重； [0013] 6)利用得到的权重，对每颗星的滤波结果进行加权，最终得到综合脉冲星时。 [0014] 所述的第2)步平均残差点位于每个区间的中点，得到等间隔分布的若干个残差点。 [0015] 所述的第4)步使用AR模型或使用周期图法估计互功率谱上界。 [0016] 所述的第4)步中如果两种方法给出的互功率谱上界一致，或两种方法给出的互功率谱上界不一致但残差中无明显低频周期成分，即在低频区间内的谱形未均呈尖峰型，则使用AR模型法；如果两种方法给出的估计不一致，且残差中存在明显低频周期成分，即在低频区间内的谱形均呈尖峰型，则使用周期图法。 [0017] 所述的第4)步中低频区间根据互功率谱上界与计算出的各互功率谱的谱形确定，取为0频率点至转折频率，转折频率指超过该频率后互功率谱上界为白噪声谱，或取为互功率谱上界主峰的半高频率区间，使得最终计算出的互功率谱在低频区间内与上界的量值之差小于设定值。 [0018] 所述的第4)步中滤波公式为 其中modulus{sk}表示从残差xk恢复出的信号sk的模值估计， 是一个单位复数，对每颗星的残差是一个常数。 [0019] 所述的第5)步向每组残差中加入一个单位周期正弦波。 [0021] 图1是不同能量噪声对残差相位的影响示意图； [0022] 图2是改进的维纳滤波算法建立综合脉冲星时流程图； [0023] 图3是使用模拟数据得到的滤波结果。 具体实施方式[0024] 本发明提出一种基于已有的频域维纳滤波法建立综合脉冲星时的方法。该方法的核心思想是根据低噪声脉冲星残差间的互相关对信号能量的估计，调整高噪声脉冲星残差互相关(互功率谱)量级，保证每组残差滤波结果的波动水平均不高于参考时信号的波动水平；同时，根据不同脉冲星残差对微小信号改变的敏感度确定权值，使综合脉冲星时能正确跟随参考时间的变化而变化。 [0025] 频域内维纳滤波的基本待求量是频率响应，它被表示为 [0026] [0027] 其中x为观测信号，s为期望信号，S代表功率谱。建立综合脉冲星时的维纳滤波法被描述为 [0028] [0029] 其中角标k代表第k颗星，xk是第k颗星的残差，sk是第k颗星的滤波结果，F是傅里叶变换算子， 是所有残差间互功率谱的算术平均。最终的综合脉冲星时被表示为[0030] [0031] μk为第k颗星残差滤波结果在加权中的权重，一般取为滤波结果均方根倒数平方的归一化值。 [0032] 理论分析与实际计算表明，通过对互功率谱取算术平均而做出的参考时信号能量估计会因为高噪声脉冲星的引入而偏大。具体而言，在传统的维纳滤波法中，对参考时信号功率谱的估计一般被表示为 [0033] [0034] N为残差数。但 [0035] [0036] [0037] 根据互谱与自谱的关系 [0038] [0039] 可知(4)式给出的互功率谱估计其实是互功率谱模值的上界，在噪声较高时会严重高于真实互谱之模，从而导致(2)式中 偏大。对于由噪声造成的互功率谱估计偏高，必需进行调整。 [0040] 另一方面，根据残差频谱相位与信号频谱相位的关系(见图1)，可知当噪声ε较低时残差x的相位能大致反映出信号的原始相位；而当噪声很高时，残差相位是信号相位与某个随机改变的叠加，此时即便功率谱准确，从残差中恢复出的信号也不是所求的参考时波动，而是加入随机扰动后的结果。 [0041] 考虑到以上两点，将(2)式重写为 [0042] [0043] 其中modulus{sk}表示从残差xk恢复出的信号sk的模值估计， 是一个单位复数，对每颗星的残差是一个常数。下面着重考虑modulus{sk}的估计方法。 [0044] 为了准确地给出钟差信号的功率谱估计，基于互相关法，必须做出如下假设：至少有两颗脉冲星数据的噪声水平不高于参考时信号。实际上，如果参考时信号与残差噪声相比非常小，则(4)式一定会大于参考时信号的功率谱。在上述假设成立的条件下，[0045] 的实部会收敛于参考时信号的功率谱，但实际计算中其可能会是负数，因此通常会使用(4)式代替，这是一种偏大的保守估计，代表了参考时信号功率谱的上界。对于其余高噪声脉冲星，(4)式给出的估计都会高于上述两颗低噪声星给出的估计，因此需要将计算出的高量级功率谱压低，至少使其不高于前述两颗低噪声脉冲星给出的(4)式之值。这样一来，所有的互功率谱都被限制在参考时信号功率谱的上界之内。由功率谱计算出傅氏变换之模，再用(8)式滤波，恢复出的信号的波动就都被限制在参考时信号波动量级之内。 [0046] 上面所说的限制互功率谱不超过上界，是对一定的频率范围而言的。具体地说，由于所关心的钟差波动都是低频信号，因此对功率谱的限制也是对某个低频区间而言。这个低频区间，可以取为从0频率点至某个转折频率，超过该转折频率的功率谱可以被认为是白噪声谱；也可以取为功率谱低频峰值的半高频率区间。在实际计算中，具体取哪种频率区间，均需要确保由确定互功率谱上界的两组残差计算出的参考时功率谱在该区间内的量值与上界接近，以保证由这两组残差恢复出的信号具有正确的波动量级。具体问题中，不排除参考时功率谱估计结果为白噪声谱的情形，为防止恢复出的信号高频能量太大，可以将功率谱转折频率或半高频率中高频点以上的功率谱置零，再求信号模值。 [0047] 对于功率谱估计，除了上面提到的传统周期图法外，也可以使用AR模型法。AR模型功率谱估计是一种分辨率很高的估计法，能够克服周期图法方差大的缺陷，在功率谱调整上具有优势。AR模型法功率谱估计依赖于模型的选取，对于残差中含有周期信号的情形，其所估计的功率谱可能出现谱峰偏移，谱峰不准等问题。因此，在实际应用中，应使用两种方法估计互功率谱上界并进行比较，如果两种方法给出的互功率谱上界基本一致，或虽有差别，但体现出残差中没有明显周期成分，就可以使用AR模型估计功率谱。 [0048] 使用AR模型估计互功率谱时，根据模型/信号的频率特征，存在两种估计方式。 [0049] 具体来说，AR模型 [0050] [0051] 所描述的随机序列的功率谱为 [0052] [0053] 其中A(z)是复变量z的多项式，p为模型阶数，aj为模型系数(实数)，Xt代表描述残2 差的时间序列，σ是白噪声过程εt的方差，i为虚数单位，ω为圆频率(取值范围为[0,π])。模 2 型参数p，aj与σ可利用AR模型理论中的多种方法进行估计，参数估计的合理性可通过白噪声检验验证。使用AR模型估计互功率谱的第一种方法与使用周期图法类似，即基于(7)式，利用(11)式给出互功率谱。除此之外，也可以使用AR过程平稳部分的Wold系数表示，直接计算互相关函数，并求得互功率谱。具体采用哪种方法的原则与使用周期图法是类似的，即需要确保由确定互功率谱上界的两组残差计算出的互功率谱在该区间内的量值与上界接近。 使用第二种方法计算的互功率谱在滤波时最终要将可能出现的复值取为实数，具体做法为直接对复值取模或取复值的实部，并确保功率谱的非负性。 [0054] 利用(8)式进行滤波并建立综合脉冲星时，从本质上讲，不同脉冲星残差所起的作用并不相同。对于高噪声脉冲星，通过互相关给出的互功率谱估计体现的是残差中的噪声特性，而非参考时本身。只有决定所有互功率谱上界的两颗低噪声脉冲星给出的互功率谱才在最大程度上代表着参考时信号的特性，这两颗星的残差相位也与参考时信号最为接近。因此，(8)式中的modulus{sk}对决定互功率谱上界的两颗脉冲星而言，最佳取法就是这两颗星给出的互功率谱估计；对于其余高噪声脉冲星而言，计算得到的互功率谱反映的都是噪声特性，残差相位也是信号相位与随机影响叠加下的随机变量，因此，可以简单使用所有互功率谱的均值进行滤波，以得到一种平均功率谱下的滤波结果。从概率上讲，当参与建立综合脉冲星时的脉冲星数无限增加时，不同脉冲星残差恢复出的信号中的噪声成分会互相抵消，从而使真实的参考时信号得以显现。 [0055] 通过上述方式滤波后，对得到的每颗星的滤波结果进行加权时不能再使用(3)式中提到的取权方法，因为此时滤波结果的波动水平已不能反映原始残差的噪声水平。本发明提出一种新的取权方法：根据每颗星对被检信号的敏感度取权。具体而言，对每颗星的残差加入一微小信号(不失一般性，微小信号可以是某振幅的单位周期正弦波，正弦波的振幅取为参考时信号波动的量级)，制造出对比残差，计算每颗星原始残差与对比残差的相关系数，对于相关系数大的脉冲星自然应取较低的权，因为相关性高说明该星的噪声水平高，对参考时的改变不敏感。据此，定义“非相关系数” [0056] w＝1‑ρ， (12) [0057] 其中ρ是自相关系数 [0058] [0059] 归一化的w即为(3)式中的权值μ。 [0060] 所述AR模型理论下利用Wold系数计算互功率谱的方法如下： [0061] 根据AR模型理论，(10)式所描述的随机过程Xt的平稳部分可被表示为[0062] [0063] 式中 为Xt的Wold系数，可以通过下式计算 [0064] [0065] 且具有性质 [0066] [0067] 因此 绝对可和。 [0068] 设两个随机过程Xt、Yt的Wold系数与白噪声分别为 {ψj}和{εs},{εt′}，可以证明 [0069] [0070] 在AR模型理论中，两个随机过程的白噪声在s‑i≠t‑j时必然是不相关的。实际应用中，当两个随机过程的相关函数不为0时，白噪声在s‑i＝t‑j时必然应相关，否则会导致E(XsYt)＝0。定义 [0071] σc2＝E(εtεt′)， (18) [0072] 记t‑s＝τ，得到 [0073] [0074] 由于Wold系数以负指数阶收敛到0，在Xt、Yt平稳的条件下(19)式只需进行少量有限项的求和运算就可得到很高的精度。再根据维纳——辛钦定理，(19)式的傅里叶变换就是互功率谱。 [0075] (19)式中σc2是一个未知量，实际计算时只需要根据互功率谱上界调整其大小，即可得到互功率谱的估计。 [0076] 综上，本发明的技术方案包括以下步骤： [0077] 1)利用初始模型，将参考时改为TT(BIPM)，利用已有数据拟合模型直至拟合前，后残差的均方根(rms)相等，得到精确的模型参数；更改参考时为被监测参考时，得到拟合前计时残差； [0078] 2)对第1)步得到的拟合前计时残差在连续的设定时长内做平均，平均残差点位于每个区间的中点，得到等间隔分布的若干个残差点； [0079] 3)对第2)步得到的残差点拟合并去除2次项； [0080] 4)使用AR模型或使用周期图法估计互功率谱上界，互功率谱上界由噪声水平最低的两组残差确定；根据互功率谱上界计算并调整所有互功率谱，使所有互功率谱在特定的低频区间内不超过上述上界；低频区间根据互功率谱上界与计算出的各互功率谱的谱形确定，可以取为0频率点至转折频率(转折频率指超过该频率后互功率谱上界可以视为白噪声谱)，或取为互功率谱上界主峰的半高频率区间，具体问题中要使得最终计算出的互功率谱(尤其是确定互功率谱上界的两组残差计算出的互功率谱)在低频区间内与上界的量值类似；应用中具体使用AR模型或周期图法估计功率谱的判断标准是，如果两种方法给出的互功率谱上界一致，或两种方法给出的互功率谱上界不一致但残差中无明显周期成分，则使用AR模型法，如果两种方法给出的估计不一致，且残差中存在明显低频周期成分(在低频区间内的谱形均呈尖峰型)，则使用周期图法；给出所有的互功率谱估计后，对确定互功率谱上界的两组残差，使用其所计算出的互功率谱进行滤波，而对其余残差，使用所有互功率谱的算术平均进行滤波，滤波公式为式(8)； [0081] 5)向每组残差中加入一个微小信号(如单位周期正弦波)，信号的幅度取为参考时信号波动的量级估计，计算各组残差加入微小信号前后的相关系数(式(13))，并以“非相关系数”(式(12))的归一化值作为每颗星的权重； [0082] 6)利用得到的权重，对每颗星的滤波结果进行加权，最终得到综合脉冲星时。 [0083] 下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。 [0084] 1)使用脉冲星计时分析软件tempo2，参考TT(BIPM15)模拟位于MJD50188～55604之间，观测间隔为10天的6组数据，白噪声水平分别为100ns(两组)，200ns，300ns，1μs(两‑1组)。其中白噪声为100ns的第一组数据加入转折频率fc＝0.067yr ，谱指数α＝2，量级amp‑40 3 ‑1 ＝1.86×10 yr的低频噪声；白噪声为200ns的数据加入fc＝0.07yr ，α＝4，amp＝3×10‑25 3 ‑1 ‑25 3 yr的低频噪声；白噪声为300ns的数据加入fc＝0.06yr ，α＝2，amp＝5×10 yr的低频噪声。更改参考时为TT(TAI)，得到并输出参考TT(TAI)的拟合前计时残差。 [0085] 2)对第1)步得到的残差在连续的30天区间内做平均，平均残差点取在每个区间的中点，得到等间隔分布的180个残差点，位于MJD50206～55576。 [0086] 3)对第2)步得到的残差拟合并去除2次项，这等同于得到拟合自转频率及其一阶导数的拟合后残差。这去除了残差中主要的非平稳成分。 [0087] 4)使用周期图法与AR模型分别计算白噪声量级为100ns的残差的互功率谱上界，发现两种方法给出的互功率谱上界体现出残差中含有明显的低频周期成分，且存在量级与峰值漂移的差异，因此使用周期图法估计互功率谱上界并进行互功率谱计算与调整。互功率谱上界由噪声水平最低的两组残差(白噪声水平为100ns的两组残差)确定，并取该谱峰的半高频率区间为功率谱限制的频率区间。计算并调整互功率谱后，根据本发明所提出的方法对每组残差滤波。 [0088] 5)向每组模拟残差中加入振幅为200ns的单位周期正弦波，计算6组残差加入正弦波前后的相关系数，并确定每颗星的权重，结果列于表1。 [0089] 表1模拟数据加入正弦波前后残差的相关系数与确定的权 [0090] [0091] 6)利用得到的权重，对每颗星的滤波结果进行加权，最终得到的综合脉冲星时与TT(BIPM15)‑TT(TAI)的对比见图3。 |
||||||
QQ群二维码
意见反馈
意见反馈