一种X射线脉冲星光子序列的频域加权比相方法

摘要

本发明公开了一种X射线脉冲星光子序列的频域加权比相方法，解决现有技术有用信号的细节信息损失，计算复杂度高和信噪比低的问题；实现步骤是：1.对两航天器接收到的两列光子到达时间序列进行等间隔采样得到对应的光子强度序列；2.将时域内的两光子强度序列经FFT变换到频域内；3.提取两光子强度序列的相位信息，进行做差运算得到两者之间的相位差；4.对出现分段的相位差进行累积，并提取相位差关于对应频点的比值；5.对频域内相位差关于频率的比值进行能量加权及比例乘法运算即可得到时域内的归一化延迟相位。本发明能够有效地提高相位估计精度，降低计算复杂度，可用于X射线脉冲星在相对导航定位系统中的相对运动状态估计。

说明书

技术领域

本发明属于相对导航技术领域，尤其涉及一种X射线脉冲星光子序列的频 域加权比相方法。

背景技术

相对导航是航天器间执行自主交会对接、协作通信以及编队飞行等空间任 务的一项关键技术。现阶段较为成熟的相对导航系统需要地面设施的辅助和支 撑，具有自主性弱、抗干扰能力差和不适用深空等缺点。而基于X射线脉冲星 的相对导航系统因其自主性高、可靠性好、适用性广等优点，与现有的导航系 统形成优势互补，是一种全新的自主导航方法，正逐渐受到航空航天领域的高 度重视和青睐。一般地，脉冲到达时间延迟量主要通过脉冲信号的延迟相位来 获得，其精度直接决定了两航天器相对定位的精度。因此，如何获得高精度的 比相估计便成了亟待解决的问题。经典的相位估计算法是时域互相关算法和频 域TaylorFFT算法。时域互相关算法主要采用光子的经验密度函数及光子的真 实密度函数间的互相关函数来得到时间延迟量，其测量的精度取决于脉冲星信 号的采样频率；频域TaylorFFT主要用于射电脉冲星计时，其基本原理是对傅 里叶变换后的经验速率函数和真值速率函数的差值进行最小化处理，该时间延 迟量的测量精度不受时域采样频率的限制，却依赖于累积脉冲轮廓的信噪比。 随后，Emadzadeh从概率统计学角度出发，直接采用测量到的光子到达时间进行 延迟估计，无需再获取速率数据，免去了复杂的累积过程，并采用克拉美罗下 界为其测量均方根误差的下界，分析了最大似然的测量误差渐变曲线。但总的 来说，上述这些对航天器间的延迟相位方法主要采用脉冲信号的累积轮廓或直 接测量得到的光子到达时间序列作为研究对象。而事实上，在对脉冲信号轮廓 进行累积的过程中，由于平滑作用极易造成信号的细节信息损失而降低信噪比； 如果选择直接测量得到的光子到达时间序列则会加大数据对象，给运算单元带 来超重负荷而降低运算效率。

发明内容

本发明的目的在于提供一种X射线脉冲星光子序列的频域加权比相方法， 旨在有效地保留有用信号的细节信息，降低计算复杂度，提高比相精度。

本发明是这样实现的，一种X射线脉冲星光子序列的频域加权比相方法， 其特征在于，所述X射线脉冲星光子序列的频域加权比相方法是基于同一颗X 射线脉冲星源，对两航天器测量到其发射的光子序列间的归一化延迟相位的测 量。其主要实现方案为：首先对两航天器测量到的光子到达时间序列进行等间 隔采样获取其对应的光子强度序列，经傅里叶变换将其变换到频域内；然后提 取频域内光子强度序列间的相位差信息，并对出现分段的相位差进行累积，获 取相位差关于对应频点的比值；最后，对频域内相位差与各频点的比值进行能 量加权及比例乘法运算即可得到脉冲信号在时域内的归一化延迟相位。该方法 主要基于具有非齐次泊松分布特点的脉冲信号开展的实验研究，可同样适用于 与X射线脉冲星信号特性类似的周期性信号模型。又由于延迟相位与两航天器 沿脉冲星矢量方向上的距离增量之间存在的固有关系可实现对航天器间相对运 动状态的估计，因此该方法也可用于对位置估计精度要求十分严格的航天器交 互对接、编队飞行及深空探测等领域。

进一步，所述X射线脉冲星光子序列的频域加权比相方法包括以下步骤：

步骤一，已知相同观测时间内两光子到达时间序列T₁(t)、T₂(t)，进行单位 周期内等间隔采样得到对应的光子强度序列λ₁(n)、λ₂(n)(n＝1,2.,...,N)；

步骤二，对比处理器内存单元与光子强度序列长度，若N＜N₀，N₀为已知 固定的阀值，则继续步骤三、步骤四、步骤五；反之，将N分成L段时间区间， 每段包含M个周期，对每个时间区间内的光子强度序列依次执行步骤三、步骤 四、步骤五；

步骤三，对λ₁(n)、λ₂(n)进行FFT变换得到X₁(k)、X₂(k)(k＝1,2,...,N)，λ₁(n)与 λ₂(n)在时域内的超前或滞后的n₀个时间片段会反映在频域内的比相信息上；

步骤四，获取X₁(k)、X₂(k)的相位并进行做差运算得到相位差

步骤五，对相位差进行累积运算，得到累积后相位差与对应频率点的 比值s(k)；

步骤六，赋予各频点处s(k)一个对应的能量权重ω(k)得到平均比值s₀，s₀乘 以比例因子1/2π即可得到两光子强度序列在时域内的归一化延迟相位。

进一步，所述步骤三的将时域内的光子强度信息经FFT变换得到频域内的 X₁(k)、X₂(k)，按照如下步骤进行：

对航天器SC1测量得到的光子强度信息进行FFT变换得到：

$X_1\left(k\right)=FT\left[{\lambda }_1\left(n\right)\right]=\underset{n=0}{\overset{N}{\Sigma }}{\lambda }_1\left(n\right)e^{-j\frac{2\pi }{N}kn},n=1,2,...,N;k=1,2,...,N;$

对航天器SC2测量得到的光子强度信息进行FFT变换得到：

$X_2\left(k\right)=FT\left[{\lambda }_2\left(n\right)\right]\approx FT\left[{\lambda }_1(n-n_0)\right]=e^{-j\frac{2\pi }{N}{\mathrm{kn}}_0}{X}_1\left(k\right);$

式中，n表示第n个采样时间片段，k表示各时间片段对应的频率，且k和n 均为整数值；N表示采样时间片段的总数；式中采用λ₁(n-n₀)约等于λ₂(n)，是由 于探测器接收到的两光子信号来自于同一颗X射线脉冲星。

进一步，所述步骤四中分别计算频域内两光子强度的相位并提取两者的相 位差，通过如下步骤计算：

采用反正切函数分别提取出两光子强度序列在频域内各频点处的相位值 为：

对两光子强度在频域内对应频点处的相位进行做差运算，得到两者之间的 相位差值：

进一步，所述步骤五中累积相位差计算累积后相位差与对应频率点 的比值s(k)，按如下步骤进行：

第一步，分析相位差的在[-2π,2π]之间变化，即频域内的是关于 频率k的不连续函数：

由上式可见，k在频率域内实际上被分成了n₀(N/(N/n₀))段，每段内的 频率变化范围均为[0,N/n₀]；

第二步，对[0,2π]范围内的相位差进行累积，消除分段现象，得到频域区间 内连续的相位差变化：

式中，设定的个数为p，且其在频域内对应的频率为p(i)，i＝1,2,...p； 且每段频率区间内相位差的个数均为n_d＝N/n₀；

第三步，对[-2π,0]范围内的相位差进行累积，消除分段现象，得到频域区 间内连续的相位差变化：

式中，设定的个数为q，且其在频域内对应的频率为q(j)，j＝1,2,...q； 且每段频率区间内相位差的个数均为n_d＝N/n₀；

第四步，设[0,2π]及[-2π,0]区间内的相位差与对应频率的比值分别为s_p(i)和 s_n(j)，则有：

令s＝{s_p,s_q}，记得到了整个频率区间内相位差与其对应频点之间的比值。

进一步，所述步骤六中对s(k)进行能量加权、求和及比例乘法运算得到时域 内光子强度序列的比相值，按如下步骤进行：

第一步，基于各频率点处的两列光子强度的幅度信息对相应的s(k)进行能量 加权，得到一个相对的稳定的s₀如下：

$s_0=\underset{i=1}{\overset{p}{\Sigma }}{\omega }_p\left(i\right)s_p\left(i\right)+\underset{j=1}{\overset{q}{\Sigma }}{\omega }_q\left(j\right)s_q\left(j\right),i=1,2,...,p,j=1,2,...,q;$

其中：

${\omega }_p\left(i\right)=\frac{\sqrt{{X_1}^2\left(p\right(i))+{X_2}^2\left(p\right(i))}}{\underset{i=1}{\overset{p}{\Sigma }}\sqrt{{X_1}^2\left(p\right(i))+{X_2}^2\left(p\right(i))}},{\omega }_q\left(j\right)=\frac{\sqrt{{X_1}^2\left(q\right(j))+{X_2}^2\left(q\right(j))}}{\underset{j=1}{\overset{q}{\Sigma }}\sqrt{{X_1}^2\left(q\right(j))+{X_2}^2\left(q\right(j))}};$

第二步，s₀乘以比例因子1/2π即得到两光子强度序列在时域内的延迟相位:

$\frac{{\hat{n}}_0}{N}=s_0/2\pi .$

本发明的另一目的在于提供一种使用所述X射线脉冲星光子序列的频域加 权比相方法的航天器交互对接系统。

本发明的另一目的在于提供一种使用所述X射线脉冲星光子序列的频域加 权比相方法的编队飞行控制系统。

本发明的另一目的在于提供一种使用所述X射线脉冲星光子序列的频域加 权比相方法的深空探测控制系统。

本发明提供的X射线脉冲星光子序列的频域加权比相方法，与现有技术相 比，具有以下优势：

1.本发明对与X射线脉冲星信号具有类似的非齐次泊松分布特性的周期性 信号具有普适性；

2.本发明选用整个观测时间内的光子强度序列作为研究对象，无须进行轮 廓累积过程，很好地保留了脉冲信号的有效性和完整性，提高了信噪比；

3.本发明对直接观测到的光子到达时间进行了采样处理，减少了数据运算， 相比于最大似然比相估计器节省了近10倍左右的计算时间；

4.本发明在比相过程中采用了频域内能量加权平均的计算方法，削弱或防 止了噪声的干扰，对有效信号赋予了较高的权值，对比相精度的提高具有一定 的贡献，该发明相比于非线性最小方差法及最大似然法得到的比相精度提高了 近10²的数量级。

附图说明

图1是本发明实施例提供的X射线脉冲星光子序列的频域加权比相方法流 程图。

图2是本发明实施例提供的实施例的实现流程图。

图3是本发明实施例提供的相位差的分段情况示意图。

图4是本发明实施例提供的相位差累积的结果示意图。

图5是本发明实施例提供的基于加权FFT、NLS及ML在比相均方根误差 方面的对比示意图。

图6是本发明实施例提供的基于加权FFT、NLS及ML在计算复杂度方面 的对比示意图。

图7是本发明实施例提供的基于加权FFT及分段加权FFT在比相均方根误 差方面的对比示意图。

图8是本发明实施例提供的基于加权FFT及分段加权FFT在计算复杂度方 面的对比示意图。

图9是本发明实施例提供的在不同的观测时间及采样间隔数下加权FFT估 计的均方根误差变化示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例， 对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以 解释本发明，并不用于限定本发明。

下面结合附图对本发明的应用原理作详细的描述。

如图1所示，本发明实施例的

S101：对两航天器接收到的两列光子到达时间序列进行等间隔采样得到对 应的光子强度序列；

S102：将时域内的两光子强度序列经FFT变换到频域内；

S103：提取两光子强度序列的相位信息，进行做差运算得到两者之间的相 位差；

S104：对出现分段的相位差进行累积，并提取相位差关于对应频点的比值；

S105：对频域内相位差关于频率的比值进行能量加权及比例乘法运算即可 得到时域内的归一化延迟相位。

下面结合具体实施例对本发明应用原理作进一步的描述。

本发明基于X射线脉冲星的加权FFT相位估计方法，如图2所示，包括如 下步骤：

步骤1：等间隔采样光子到达时间序列得到光子强度序列。

在相同的观测时间内，记录两航天器各自接收到的光子到达时间序列T₁(t)、 T₂(t)，对其进行等间隔时间采样(设采样时间片段为T_b，则采样间隔个数 N_s＝P/T_b，且满足N_s＝2ⁿ)，得到强度序列λ₁(n)、λ₂(n)(n＝1,2.,...,N)，两光子强 度序列可反映出脉冲星辐射的光子的概率密度信息及不同航天器接收到的光子 到达时间延迟量的信息。

步骤2：判断光子强度序列的长度是否满足处理单元的内存需求。

由于在有限的内存单元中对大量的光子强度数据进行FFT处理容易由于过 载而发生溢出，根据经验设定一个序列长度阀值N₀，当采样得到的λ₁(n)、λ₂(n)的 长度N≤N₀时，继续进行下面的步骤，反之则将时间长度N分成L段时间区间， 每段区间内包含M个周期，对每段时间区间内的光子强度序列依次进行以下各 步骤的操作。

步骤3：将时域内光子强度信息经FFT变换到频域。

设采样得到的光子强度序列中λ₂(n)相比于λ₁(n)超前或滞后n₀的时间片段， 则将其归一到[0,1]内的延迟相位为n₀/N。对两光子强度序列进行FFT变换可以 得到它们在频域内的幅值及相位的变化特性。

$X_1\left(k\right)=FT\left[{\lambda }_1\left(n\right)\right]=\underset{n=0}{\overset{N}{\Sigma }}{\lambda }_{1}h\left(n\right)e^{-j\frac{2\pi }{N}kn},n=1,2,...,N;k=1,2,...,N$

$X_2\left(k\right)=FT\left[{\lambda }_2\left(n\right)\right]\approx FT\left[{\lambda }_1(n-n_0)\right]=e^{-j\frac{2\pi }{N}{\mathrm{kn}}_0}{X}_1\left(k\right)$

式中，n表示第n个采样时间片段，k表示各时间片段对应的频率，且k和n 均为整数值；N表示采样时间片段的总数。式中采用λ₁(n-n₀)约等于λ₂(n)，是由 于探测器接收到的两光子信号来自于同一颗X射线脉冲星，理想情况下两光子 信号具有相同的轮廓特征。同时，可以观察到时域内两光子强度序列间超前或 滞后的时间片段反映在频域内的相位信息上。

步骤4：提取频域内两光子强度的相位进行做差运算。

采用反正切函数分别提取出两光子强度序列在频域内各频点处的相位值 为：

对两光子强度在频域内对应频点处的相位进行做差运算，可以得到两者之 间的相位差值：

步骤5：累积相位差并得到其关于各频率的比值s(k)。

由FFT变换的时延特性，两光子强度序列在时域内的超前或滞后的时间片 段n₀会与步骤4所求得的相位差之间存在一定的等式关系。而实际上，又由于 通过反正切函数求解得到的相位值的范围为[-π,π]，相应得到的相位差的会在[-2π,2π]之间变化，即频域内的是关于频率k的不连续函数，如图3所 示。此时，有：

由上式可见，k在频率域内实际上被分成了n₀(N/(N/n₀))段，每段内的 频率变化范围均为[0,N/n₀]。

为了获得关于k的连续变化关系，分别对[0,2π]及[-2π,0]范围内的相位 差进行累积运算。相位累积是对在[0,2π]及[-2π,0]两个值域区间的相位差分别进 行的用于消除分段现象的一个叠加过程。假定的个数为p，的 个数为q，且它们在频域内对应的频率分别为p(i)、q(j)(p(i),q(j)∈[1,N])，其中 i＝1,2,...p，j＝1,2,...q。又由相位差在频域内的分段特性容易知道：每段频率区间 内相位差的个数均为n_d＝N/n₀，则对所有位于[0,2π]或[-2π,0]内的相位差值分别 进行累积即可得到累积后的连续的相位差变化，如图4所示。[0,2π]、[-2π,0]内 连续的相位差表示为：

设[0,2π]及[-2π,0]区间内的相位差与对应频率的比值分别为s_p(i)和s_n(j)，则 有：

令s＝{s_p,s_q}，即得到了整个频率区间内相位差与其对应频点之间的比值。

步骤6：对s(k)进行能量加权及比例乘法运算得到时域内光子强度序列的延 迟相位。

理论上，s(k)会与两光子强度序列间的延迟时间片段n₀呈线性关系。然而， 由于光子信号在发射、传输及接收过程中总是会受到外来噪声的干扰，该线性 关系不会十分显著。这里，基于各频率点处的两列光子强度的幅度信息对相应 的s(k)进行能量加权，得到一个相对的稳定的s₀如下：

$s_0=\underset{i=1}{\overset{p}{\Sigma }}\omega \left(i\right)s_p\left(i\right)+\underset{j=1}{\overset{q}{\Sigma }}\omega \left(j\right)s_q\left(j\right),i=1,2,...,p,j=1,2,...,q$

其中

$\omega \left(i\right)=\frac{\sqrt{{X_1}^2\left(p\right(i))+{X_2}^2\left(p\right(i))}}{\underset{i=1}{\overset{p}{\Sigma }}\sqrt{{X_1}^2\left(p\right(i))+{X_2}^2\left(p\right(i))}},\omega \left(j\right)=\frac{\sqrt{{X_1}^2\left(q\right(j))+{X_2}^2\left(q\right(j))}}{\underset{j=1}{\overset{q}{\Sigma }}\sqrt{{X_1}^2\left(q\right(j))+{X_2}^2\left(q\right(j))}}$

这里，采用能量进行加权的物理意义主要有两点：(1)使s(k)从一个不稳定 的值收敛到一个稳定值，便于延迟相位的求解；(2)在一定程度上能够有效地 提取出在整个观测时间内的有效光子信号，减弱边缘噪声对有效信号的影响。

最后，s₀乘以比例因子1/2π即可得到两光子强度序列在时域内的延迟相位。

$\frac{{\hat{n}}_0}{N}=s_0/2\pi$

此外，如果在步骤2采用的是分段加权FFT，则经过对多段区间内的光子 强度序列进行了步骤3～6的操作，则此时两光子强度序列在时域内的延迟相位 为：

$\left(\frac{{\bar{n}}_0}{N}\right)=\frac{1}{L}\left(\underset{i=1}{\overset{L}{\Sigma }}\right({s_0}^{\left(i\right)}/2\pi \left)\right)$

这是经过多次平滑后的结果，且由于其每段内数据对象均小于直接加权FFT 相位估计，因此必然存在关系：

$\frac{{\bar{n}}_0}{N}<\frac{{\hat{n}}_0}{N}$

虽然分段加权FFT进行相位估计的效果没有直接加权FFT的好，但却能很 好的解决加权FFT对内存的苛刻要求，是加权FFT很好的补充和替代。在工程 实践过程中，可以采用加权FFT和分段加权FFT进行交替使用来实现高精度的 比相估计。

下面结合仿真实验对本发明的应用效果作详细的说明。

本发明的应用效果可以通过以下对采用加权FFT、非线性均方差法(NLS) 及最大似然法(ML)估计比相进行对比仿真实验来进一步说明：

1.仿真实验条件

本发明的仿真在Window7专业版，CPU为Intel(R)Core(TM)i5-3470， 内存为4.00GB，软件平台为MatlabR2011a上实现。仿真实验选用蟹状脉冲星 (PSRB0531+21)作为X射线源，测试其发射的脉冲信号在两航天器间的比相 求解。其中，蟹状脉冲星发射脉冲信号的频率为29.86866613087979871466，航 天器上安装的探测器的表面积为10000平方厘米；X射线源发射的光子序列及 探测器的背景噪声主要采用由高斯和模型和泊松模型拟合生成，设置光子流量 密度为0.064092419551171光子/秒/平方厘米，背景噪声密度为 0.445211779062184光子/秒/平方厘米。

2.仿真内容与结果

仿真一，用加权FFT、NLS及ML算法分别估计出两航天器接收的光子信 号之间的相位值，求解估计相位与标准相位之间的均方根误差进行对比，如图5 所示。

从图5可以看出，随着观测时间的延长，NLS、ML及加权FFT估计的比相 的均方根误差均呈现逐渐减小的趋势，且加权FFT比前两者估计出的比相的均 方根误差要低10²的数量级左右，可知其比相精度要明显高于NLS和ML算法。

仿真二，统计加权FFT、NLS及ML相位估计器在进行相位估计的过程中 所占用的CPU时间，如图6所示。

从图6可以看出，随着观测时间的延长，加权FFT算法的计算代价介于NLS 与ML算法之间，但加权FFT仍然比ML占用CPU的时间要少10倍左右，要 明显比ML节约计算时间。

仿真三，用加权FFT及分段加权FFT分别估计出两航天器接收的光子信号 之间的相位值，求解估计相位与标准相位之间的均方根误差进行对比，如图6 所示。

从图7可以看出，初始时刻，分段加权FFT比相的均方根误差几乎与该时 间区间内进行加权FFT比相所得到的RMS值接近。但随着观测时间的延长，由 于测试数据的增多，比相均方根误差会被多次平均，其比相的RMS值会呈现缓 慢的下降趋势，虽然下降的速度明显要低于加权FFT比相，总体上来说仍然高 于ML和NLS算法。

仿真四，统计加权FFT分段加权FFT相位估计器在进行相位估计的过程中 所占用的CPU时间，如图8所示。

从图8可以看出，随着观测时间的延长，分段加权FFT的CPU时间仍然保 持着显著增长的趋势，但其在每个观测点处所使用的CPU时间明显要少于直接 加权FFT比相所使用的时间。

仿真五，采用控制变量法及蒙特卡罗法对不同的观测时间及采样间隔数求 解加权FFT估计的相位与标准相位之间的均方根误差，如图9所示。

从图9可以看出，在一定的采样间隔下，当观测时间逐渐延长时，加权FFT 的比相精度逐渐提高；在一定的观测时间下，当采样间隔逐渐增大时，加权FFT 的比相精度逐渐提高。

综上，加权FFT明显比NLS及ML相位估计器具有更高的比相精度，虽然 其所占CPU时间在NLS及ML之间，但仍然比ML要小10倍左右。分段加权 FFT能够很好地改善加权FFT对内存的苛刻要求，其得到的比相精度虽然要比 加权FFT低，但基本上比NLS及ML算法均高，且分段加权FFT的总CPU时 间比直接加权FFT的低。采样间隔数及观测时间的增大均对加权FFT估计的比 相精度有改善作用。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发 明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明 的保护范围之内。