基于X射线脉冲星的最大似然相位估计方法

摘要

本发明公开了一种基于X射线脉冲星的最大似然相位估计方法，包括：根据测量击打到X射线探测器探测材料上的光子到达时间的过程，构造X射线脉冲星信号泊松模型构造单元；利用多个高斯分布概率函数对X射线脉冲星轮廓进行拟合构造X射线脉冲星信号模型构造单元；构造用来估计初始相位和周期的最大似然估计单元；构造提高新模型中代价函数的峰值搜索性能的相位并行计算单元。既降低了计算量，减少了计算时间，又保证了较高的运算精度，在X射线脉冲星导航的研究中具有重要的意义。

说明书

技术领域：

本发明属于导航技术领域，特别涉及X射线脉冲星产生信号的相 位估计，用于近地轨道及深空中飞行器的位置和速度的确定。

背景技术：

早在1974年，由于能产生频率稳定的脉冲信号，人们考虑脉冲 星做定位使用。近年来，许多研究证明了基于脉冲星到达时间的导航 是可行的。X射线脉冲星其辐射出的X射线含有其辐射的绝大部分 的能量，从而只需要小体积和轻重量的X射线探测器。现有的全球 导航卫星系统工作也可为近地轨道飞行器进行自主导航，但是由于它 们的导航信号深空中无法获取，这些系统在深空中工作受到限制。X 射线脉冲星除了分布于整个银河系外，还能通过测量时间或相位提供 精确的定位，因此它在深空中可以全面进行自主导航。

近年来，许多学者都在不同方面对X射线脉冲星导航技术进行 研究。研究表明，如何精确测量X射线脉冲的到达时间成为X射线 脉冲星导航技术最关键的一点。Hanson等人提出了一种基于高斯分 布相位估计的周期折叠方法。Emadzadeh等人完善了实验结果，并建 立了周期折叠方法的数学模型。除此之外，Emadzadeh等人基于脉冲 星光子到达时间的统计，提出了基于累积轮廓的非线性最小二乘法和 最大似然法来对脉冲星信号的相位进行估计

上述方法的前提条件是要求累积出信号轮廓，轮廓累积的过程实 际上是信号平均的过程，这一过程中会损失部分有用信息，从而导致 信号的测相精度下降，而且累积过程屏蔽了轨道运动所形成的多普勒 频率变化，降低了信号信息利用率。

发明内容：

本发明的目的在于克服上述现有技术的缺点，提供一种基于X射 线脉冲星的最大似然相位估计方法，包括：

根据测量击打到X射线探测器探测材料上的光子到达时间的过 程，构造X射线脉冲星信号泊松模型构造单元；

利用多个高斯分布概率函数对X射线脉冲星轮廓进行拟合构造X 射线脉冲星信号模型构造单元；

构造用来估计初始相位和周期的最大似然估计单元；

构造提高新模型中代价函数的峰值搜索性能的相位并行计算单 元。

所述的X射线脉冲星信号泊松模型构造单元包括基于泊松分布的 时序模型，其构造如下：脉冲星具有独特且稳定的与脉冲相位有关 的完整脉冲轮廓，假设是标准脉冲轮廓，λ_b和λ_s分别是有效的噪 声强度和流量强度；同时，探测器A的有效区域，整个观察期间K 的叠加时间以及探测器效率η也都与该模型相联系；因此，X射线脉 冲星的到达率模型可以表示成轮廓函数：

其中T_i是采样间隔，x是参考点，v是探测器速率，λ_s是X射线脉冲 星信号辐射强度比例因子；是探测器向辐射源运动形成的无直流 的标准脉冲星轮廓，为具有多普勒情况下的相位；。

由于T_i是非重叠性的时间间隔，落在时间间隔T_i中的光子k_n的周期 独立的泊松分布，如下：

其中k_n=1,2,3…，所描述的泊松分布的均值和方差为：

对于时间序列它的联合概率分布函数表示为：

所述X射线脉冲星信号模型构造单元包括多个高斯分布概率曲线 拟合模块，采用多个高斯分布概率函数来拟合X射线脉冲星轮廓； 通过一组关于标准平均脉冲轮廓的高斯函数来描述：

其中i代表第i分量，a_i，μ_i和δ_i分别为第i高斯分量的比例因子，均 值和方差；

所述；X射线脉冲星信号模型构造单元包括基于GFSAP方法的新 模型构造模块，其新模型构造模块如下：提出了一个假设，对于脉冲 星光子序列，如果n个光子在同一周期被捕获到，理论上是它就可以 看成是单个概率分布函数为标准脉冲星轮廓的光子的独立事件，在这 一周期发生了n次；

假设代表标准脉冲星轮廓，且如果在一个脉冲星周 期内仅有一个光子，它的到达时间同样遵循概率分布函数为的分 布；

令λ_b＝0，那么，脉冲星信号密度可以重新表示成 其中Δ是开始时间到光子到达时间的时间间隔； 令P为脉冲星周期，如果在周期P内有一个光子到达，在间隔Δ内这 个时间发生的概率为

此外，由于则有得出一个周 期内事件在时刻τ发生的概率如下

P(τ)＝h(τ)；

令τ_n代表第n个光子的小数部分，N代表循环计数；则有

${\tau }_n=\frac{t_i-\mathrm{NP}}{P}$

如果脉冲星频率由于受航天器速率影响是个未知的常量，那么周期要 修正为

${\tau }_n=\frac{t_i-{\mathrm{NP}}_v}{P_v}$

其中由上述知，v是航天器的速率；这样，在一个周期内 光子的到达时间的概率可以重新表达成

P(τ_n)＝h_g(τ_n)  τ_n∈[0,1)；

此式表明，第N个周期的光子到达时间的分布可以近似为用标准轮 廓作为它的分布函数时的概率；为了定义初始相位脉冲星轮廓的 高斯函数可重新写为

其中δμ_i是初始相位和脉冲星轮廓的第i分量之间的偏差； 可以设置为[0,1)之间的任意值。

所述最大似然估计单元包括最大似然估计模块，其中最大似然估 计模块如下：对于序列的概率分布函数通过 找关于参数的最大值进行最大似然估计；即将该函数的自然对数最 大化，如下

其中LLF表示对数似然函数；那么，初始相位可以通过解决如下优 化问题来估计

假如周期P是一个未知常量，它也可以通过新型最大似然估计方 法来估计；在这种情况下，周期P和初始相位是下述函数的两个变 量：

并且未知量可以通过下列公式求解：

其中Γ是P_v的搜索空间。

GFSAP模型必须使用多个高斯分量来拟合含有多个峰值的脉冲 星轮廓，利用牛顿-拉普森迭代搜索方法，它的迭代因子如下

${\overrightarrow{\theta }}_{k+1}={\overrightarrow{\theta }}_k+{\left[J\left(\overrightarrow{\theta }\right)\right]}^{-1}\mathrm{pLLF}\left(\overrightarrow{\theta }\right){|}_{\theta ={\theta }_k}$

其中，k=0,1,2,…，当||θ_k+1-θ_k||＜ξ成立时 迭代过程将会停止，其中ξ是收敛极限；然而，由于观察时间长，等 式（29）中的将会产生很大的计算量，因为对于每一个光子 的到达时间，都必须计算非线性高斯求和函数；由于所有光子的到达 时间都有相同的概率分布函数；因此，可以重新改写为

$\mathrm{pLLF}\left(\overrightarrow{\theta }\right)=\underset{l=1}{\overset{L}{\Sigma }}i=\underset{i=\mathrm{la}+1}{\overset{\mathrm{la}+a}{\Sigma }}\frac{\partial \ln \left(h_g({\tau }_n;\overrightarrow{\theta })\right)}{\partial \overrightarrow{\theta }}$

此式表明，计算过程可以分解为L个含有a个光子的独立子过程；据 此，计算的步骤如下：

步骤1、将光子序列划分为L个部分，每个部分含有[m/L]个光子, 其中[·]表示循环操作；

步骤2、对于步骤1中的每一部分，计算函数；

步骤3、将步骤2中的结果进行相加。

本发明的有益效果在于：既降低了计算量，减少了计算时间， 又保证了较高的运算精度，在X射线脉冲星导航的研究中具有重要 的意义。

附图说明：

图1是X射线脉冲星信号泊松模型构造单元流程图；

图2是新型X射线脉冲星信号模型构造单元流程图；

图3是新型最大似然估计单元流程图；

图4是最大似然估计与克拉美罗界对比图；

图5是本发明的的并行计算方法的流程图。

具体实施方式：

下面结合附图对本发明做进一步详细描述：

参照图1，该发明X射线脉冲星信号泊松模型构造单元中的基于 泊松分布的时序模型，其构造方法如下：由于辐射周期的高稳定性， X射线脉冲星辐射可以看成是一种周期平稳的过程。脉冲星的平稳让 我们能在任何时刻在太阳系中心（SSB）参考坐标中准确预测脉冲星 信号相位。在这个坐标系中，脉冲星信号的相位可以用下列表达式表 达：

${\Phi }_n^{\mathrm{SSB}}\left(t\right)={\Phi }_n^{\mathrm{SSB}}\left(t_0\right)+f_n·(t-t_0)+O\left(m\right)---\left(1\right)$

其中是在通用时刻t时第n颗星相位，f_n是第n颗脉冲星的脉 冲频率，O(m)是的高阶项，由于高阶项仅仅影响长期范围，我 们可以把O(m)看成是无关紧要的变量并且忽略它对相位的影响。为 了进一步简化表达式，我们用来替换

探测器的多普勒速率会产生多普勒频率f_d，它可以分解为两个不 同的部分：X射线辐射源的频率f₀和多普勒频移f₀·ν/c。即 f_v＝f₀·(1+ν/c)    （2）

其中ν是探测器的速率，c是光速。因此，假设是观察到的相位 且探测器的速度是常量，那么探测到的相位等于：

其中x₀的探测器的初始位置，是初始相位。

脉冲星具有独特且稳定的与脉冲相位有关的完整脉冲轮廓。假设 是标准脉冲轮廓，λ_b和λ_s分别是有效的噪声强度和流量强度。同 时，探测器A的有效区域，整个观察法K的叠加时间以及探测器效 率η也都与该模型相联系。X射线脉冲星的到达率的模型可以表示成 轮廓函数：

其中T_i是采样间隔，x是参考点，v是探测器速率，λ_s是X射线脉冲 星信号辐射强度比例因子。是探测器向辐射源运动形成的无直流 的标准脉冲星轮廓。

由于T_i是非重叠性的时间间隔，落在时间间隔T_i中的光子k_n遵循独 立的泊松分布，如下：

其中k_n=1,2,3…..等式（5）描述的泊松分布的均值和方差为：

对于时间序列它的联合概率分布函数表示为

参照图2，本发明的多个高斯分布概率曲线拟合模块，提出了一 种用多个高斯分布概率函数来拟合X射线脉冲星轮廓的方法。这种 方法能够很好地表示出脉冲星轮廓的结构和细节。因此，我们通过一 组关于标准平均脉冲轮廓的高斯函数来描述这种方法

在等式（8）中，

其中i代表第i个成分，a_i，μ_i和δ_i分别为第i个高斯分量的比例因子， 均值和方差。

参照图2，本发明基于GFSAP方法的新模型的构造方法如下： 假设记录数据从t₀开始，t_end代表结束时间，那么观测时间间隔为 t_obs＝t_end-t₀。此外，t_i代表第i个光子到达时间，光子序列可以表示成 所述的泊松模型的定义表明记录的序列的到达时 间是递增的，即

t₀＜t₁＜t₂＜t₃＜...＜t_m＜t_end    （10）

我们提出了一个假设，对于脉冲星光子序列，如果n个光子在同 一周期被捕获到，它就可以看成是一个概率分布函数为标准脉冲星轮 廓的光子的独立事件，理论上是在这一周期发生了n次。这种假设由 以下定理支持。

定理1、假设代表标准脉冲星轮廓，且如果在一个 脉冲星周期内仅有一个光子，它的到达时间同样遵循概率分布函数为 的分布。

论证：不失一般性，我们令λ_b＝0。那么，根据等式（4），脉冲星 信号密度可以重新表示成其中Δ是开始时间到光 子到达时间的时间间隔。令P为脉冲星周期。如果在周期P内有一个 光子到达，在间隔Δ内这个时间发生的概率为

此外，由于则有我们将很容 易得出一个周期内事件在时刻τ发生的概率如下

P(τ)＝h(τ)    （12）

因此，定理1成立。

根据定理1和等式（12），我们将提出一种研究脉冲星信号的新型 方法。令τ_n代表第n个光子的小数部分，N代表循环计数。则有

${\tau }_n=\frac{t_i-\mathrm{NP}}{P}---\left(13\right)$

如果脉冲星频率由于受航天器速率影响是个未知的常量，那么周期要 修正为

${\tau }_n=\frac{t_i-{\mathrm{NP}}_v}{P_v}---\left(14\right)$

其中由等式（2）知，v是航天器的速率。这样，根据定 理1和等式（8），在一个周期内光子的到达时间的概率可以重新表达 成

P(τ_n)＝h_g(τ_n)  τ_n∈[0,1)（15）

等式（15）表明，第N个循环周期的光子到达时间的分布可以近 似为用标准轮廓作为它的分布函数时的概率。为了定义初始相位等式（9）可重新写为

其中δμ_i是初始相位和脉冲星轮廓第i部分之间的偏差。 总而言之，可以设置为[0,1)之间的任意值。在下述部分，我们令 即δμ₀=0。

参照图3，本发明采用的新型最大似然估计方法，是一种用来估 计初始相位和时间的最大似然估计方法，其中初始相位的最大似然估 计方法如下：

对于序列的概率分布函数

我们可以通过找关于参数的最大值进行最大似然估计。即将该函数 的自然对数最大化，如下

其中LLF表示对数似然函数。那么，初始相位可以通过解决如下优 化问题来估计

周期最大似然估计方法如下：假如周期P是一个未知常量，它 也可以通过新型最大似然估计方法来估计。在这种情况下，周期P和 初始相位是下述函数的两个变量：

并且未知量可以通过下列公式求解：

其中Γ是P_v的搜索空间。

参照图4，本发明采用的新型最大似然估计方法具有较高的估计 性能，其性能分析如下：克拉美罗界是任一未知参数的无偏估计变量 的下界。它是一个在很小的误差范围内估计算法（如最大似然估计算 法）性能的有效且严谨的下界。令为观测到的向量τ_n的概率分 布函数，未知参数为则的费雪信息量可表示为

其中假设满足规律性条件

$\int \frac{\partial p(x;\overrightarrow{\theta })}{\partial \overrightarrow{\theta }}\mathrm{dx}=0---\left(22\right)$

对矢量无偏估计的克拉美罗界可表示为

$\mathrm{cov}\left(\hat{\theta }\right)\ge J^{-1}\left(\overrightarrow{\theta }\right)---\left(23\right)$

对于相位费雪信息量可以表示成

相似地，速率v的费雪信息量表示为

$J\left(v\right)=-\underset{i=1}{\overset{N}{\Sigma }}\frac{1}{\lambda (T_i;\overrightarrow{\theta })}{\left(\frac{\partial \lambda (T_i;\overrightarrow{\theta })}{\partial v}\right)}^2---\left(25\right)$

由于λ(·)是个周期函数，相位和速率的克拉美罗界经过一些推导 后分别表示为

和

$\mathrm{CRLB}\left(v\right)={\left(\frac{\mathrm{AK\eta }{\lambda }_s^2{f_0}^2}{3P{c}^2}{T_{\mathrm{obs}}}^3{T}_i\Omega \right)}^{-1}---\left(27\right)$

其中t_obs是观测时间，P是脉冲星周期，c是光速，且

其中C是排列组合。等式（26）表明正比于脉冲周期而反比 于有效面积A，探测器效率η，流量的平方观测时间T_obs，采样时 间间隔T_i，脉冲频率的平方f₀²和轮廓因子Ω。与（26）相比较，等式 （27）中大多数CRLB(v)的参数是相同的，除了CRLB(v)反比于T_obs³/3。 等式（28）表明，基于多个高斯概率分布曲线拟合方法，脉冲星轮廓 对克拉美罗界的贡献与方差的乘积δ_kδ_j和每个变量的相关关系 f_k(·)f_j(·)有关。它表明要获得更高的导航性能或者降低卡拉美罗界， 我们应该采用有较低δ_kδ_j或者较高f_k(·)f_j(·)的脉冲星轮廓，并且，含有 更多分量的脉冲星轮廓将可能获得更好的估计性能。

参照图5，本发明中新型最大似然估计单元中的数值搜索方法， 能够提高新模型中代价函数的峰值搜索性能。通常的，GFSAP模型 必须使用多个高斯分量来拟合含有多个峰值的脉冲星轮廓。因此，价 值函数一般来说不是凸起的即含有多个极小值。为了避免困于局部极 值，网格离散搜索方法将用在直接搜索最大值的过程中。这种方法可 以权衡计算复杂性和估计精确性。考虑到相位估计通常用在迭代定位 过程中，先验的相位信息可以容易的从先前测定的位置中获得。因此， 我们可以利用牛顿-拉普森迭代搜索方法，它的迭代因子如下

${\overrightarrow{\theta }}_{k+1}={\overrightarrow{\theta }}_k+{\left[J\left(\overrightarrow{\theta }\right)\right]}^{-1}\mathrm{pLLF}\left(\overrightarrow{\theta }\right){|}_{\theta ={\theta }_k}---\left(29\right)$

其中，k=0,1,2,…，当||θ_k+1-θ_k||＜ξ成立时 迭代过程将会停止，其中ξ是收敛极限。然而，由于观察时间长，等 式（29）中的将会产生很大的计算量，因为对于每一个光子 的到达时间，都必须计算非线性高斯求和函数。该模块提出的一种 的并行计算方法来降低计算的复杂性，由于第三节提出的模 型假设X射线脉冲星信号是周期平稳过程，所有光子的到达时间都 有相同的概率分布函数。因此，可以重新改写为

$\mathrm{pLLF}\left(\overrightarrow{\theta }\right)=\underset{l=1}{\overset{L}{\Sigma }}i=\underset{i=\mathrm{la}+1}{\overset{\mathrm{la}+a}{\Sigma }}\frac{\partial \ln \left(h_g({\tau }_n;\overrightarrow{\theta })\right)}{\partial \overrightarrow{\theta }}---\left(30\right)$

等式（30）表示，计算过程可以分解为L个含有a个光子的独立子过 程。根据式（30），计算通过如下步骤：

步骤1、将光子序列划分为L个部分。每个部分含有[m/L]个光 子,其中[·]表示循环操作。

步骤2、对于步骤1中的每一部分，计算函数。

步骤3、将步骤2中的结果进行相加。

本发明提出了X射线脉冲星的一种新型最大似然相位估计方法， 可以直接利用所测量的到达时间进行相位估计，并提出了一种并行最 大似然估计方法来更有效率的求解最大似然问题。本专业领域人员在 不背离本发明权利要求范围和主旨的前提下可以实现多种显而易见 的改进，本发明的权利要求范围并不限于以上论述。

本发明未做详细描述的内容属于本领域专业技术人员公知的现 有技术。