基于Neyman-Pearson準則的空間目標軌道機動檢測方法

王慶瑞，鄒江威，吳文振，陳健

國防科技大學 自動目標識別重點實驗室，長沙 410073

隨著世界各國對太空探索的逐步深入，航天任務日益頻繁，空間目標種類與數量不斷增加，太空環境變得愈加復雜[1-2]，快速且準確的軌道機動檢測成為保障空間目標運行安全的重要需求[3-4]。

國內外學者從不同角度對空間目標機動檢測進行了研究。Kelecy等[5]利用最小二乘法與卡爾曼濾波法研究了脈沖小推力作用下的機動檢測，并分析比較了兩種方法的性能，但缺少可操作的檢測流程。Lubey等[6]利用最優控制估計方法，通過計算控制距離度量實現脈沖推力與連續推力作用下的機動檢測，其計算過程中的動態不確定性參數需要手動調整。國內董云峰等[7]率先提出利用小波分析工具實現軌道機動檢測，根據小波系數曲線隨小波分析尺度的變化趨勢判定是否存在軌道機動。張振軍等[8]通過比較小波分解后信號的相關系數閾值來判斷軌道機動。蘇建敏等[9]通過將小波分解后信號的平均值與給定閾值進行對比，檢測目標衛星軌道機動。上述利用小波分析工具檢測軌道機動，對于如何獲取合適的判斷閾值，缺乏定量的推導與分析，實用性尚有欠缺。崔紅正等[10]利用觀測值殘差比對不同推力下的機動進行了研究，對工程應用有較好的參考價值。

針對脈沖小推力下的軌道機動檢測問題，本文提出了一種基于Neyman-Pearson (NP)[11-12]概率判決模型的檢測方法。該方法能夠根據某段測量數據的概率密度函數與設定的虛警概率，自適應地生成判決門限，實現衛星變軌的有效檢測。采用滑窗策略，該方法能夠對整段測量數據實現持續性檢測，得到變軌發生的大致時間。文中首先介紹了變軌檢測的原理，然后給出了相應的變軌檢測步驟，最后通過不同機動幅度下目標衛星變軌檢測的仿真實驗，驗證了該方法的有效性。

1 特征量建模

利用星載平臺對目標衛星進行觀測的場景如圖 1所示，伴飛衛星與目標衛星相隔一定距離，同方向飛行，利用伴飛衛星上搭載的激光測距雷達[13]測量與目標衛星的相對距離，記為r0。

圖1 激光雷達測距示意圖Fig.1 Schematic diagram of lidar ranging

考慮到實際中星載激光雷達存在一定的測量誤差，測量得到的目標衛星相對距離為:

r=r0+w

(1)

式中：w為隨機測量誤差，設w滿足高斯同分布，標準差為σ。目前工程上使用的星載激光雷達測距精度已經達到毫米量級[14]，但考慮到實際中跟蹤誤差等不利因素的影響，本文設定測量噪聲的標準差為σ=0.05 m。

對相對距離r進行差分處理得到距離變化率V為：

(2)

當目標衛星未發生軌道機動時，相對距離在一定時間內是連續變化的，距離變化率V趨近平滑。當目標衛星發生軌道機動時，由于目標衛星受到脈沖推力作用，其速度發生突變[10]，即目標衛星相對伴飛衛星的距離變化率發生突變，此時距離變化率具有階躍信號特征[8]，因此可以將距離變化率作為特征量用于變軌檢測。

2 機動檢測方法原理

目標衛星軌道機動的檢測問題可抽象為對信號V的階躍突變檢測問題。但由于測量噪聲的存在，使得信號V的階躍特征被淹沒在噪聲中，無法直接檢測出來。本節以NP檢驗模型為基礎，構建機動檢測方法，檢測包含測量噪聲的信號V是否發生階躍突變。

2.1 NP檢驗模型

NP準則可用于檢驗某信號序列是否發生整體性偏移。設該信號序列為x，序列中各信號的概率密度分布函數已知，不妨設他們滿足等方差的正態分布。基于NP準則對序列x構造假設性檢驗如下：

(3)

式中：H0是零假設;H1是備選假設;幅度A表示序列的整體偏移，為未知常數；序列v服從期望為μ[n](n=1,2,…,N)、等方差σ2的獨立正態分布。序列x在H1、H0條件下的概率密度函數分別記為p(x;A,H1)與p(x;H0)。定義似然比：

(4)

似然比描述了對于序列x，H1的可能性與H0的可能性的比值。如果

(5)

那么檢測器判決為H1,其中γ為門限。將式(4)代入式(5)，兩邊取對數得：

假定A>0，化簡得:

將上式乘以系數1/N，得到序列x的均值，記為檢驗統計量T(x)：

(6)

記

(7)

其中γ′為新門限，得到：

(8)

式(7)說明了序列偏移量A與門限γ′的定量關系，但由于A為未知常數，無法通過式(7)解算出γ′。從另外一個方面，可以根據T(x)的概率密度函數，結合虛警概率PFA約束求解門限γ′。由于序列x為等方差的獨立正態分布，其均值T(x)也服從正態分布，T(x)的期望與方差分別為：

根據T(x)的正態分布特征，得到虛警概率PFA、檢測概率PD和門限γ′的關系為：

(9)

(10)

根據式(9)，在虛警概率約束下求解門限γ′：

(11)

當T(x)>γ′，判決為H1，此時為正確檢測；判決為H0，此時為虛警。將式(11)代入式(10)可得：

(12)

式中:Q為右尾函數，單調遞減；Q-1為其逆函數，單調遞增。上式表明，對于給定的PFA，檢測性能隨著NA2/σ2單調遞增，這里將NA2/σ2定義為階躍能量噪聲比，用來描述階躍能量和信號方差的大小關系。

同理，當A<0時，給定虛警概率下的檢測門限和檢測概率分別為：

(13)

PD=Pr {T(x)<γ″;H1}

(14)

值得注意的是，式(11)與式(13)的第一項互為相反數，第二項相同。這是因為T(x)服從正態分布，其概率密度函數是左右對稱的，對稱中心為T(x)的期望。因此A>0與A<0時的檢測問題是等效的，兩種情況下的檢測概率表達式相同。

2.2 衛星軌道機動檢測方法

依據上文的NP模型，采用距離變化率序列{V(n),n=1,2...,S}對衛星的脈沖小推力軌道機動進行檢測。從整段數據序列中抽取一段作為待檢測數據，將待檢測的某段數據分為前后兩段，假定前半段數據中目標衛星未發生變軌，作為參考窗口；后半段數據中目標衛星可能發生了變軌，作為檢測窗口。通過對參考窗口數據進行多項式擬合、外推，然后與檢測窗口軌道數據特征進行對比，實現目標衛星軌道機動檢測。

首先在序列V(n)確定位置相鄰的參考窗口和檢測窗口，長度分別為M和N，它們的關系設定為：

參考窗口：R(n)=V(n),n=n0,n0+1,n0+2,n0+3, …,n0+M-1。

檢測窗口：D(n)=V(n),n=n0+M,n0+M+1,n0+M+2, …,n0+M+N-1。

設參考窗口中目標衛星未發生變軌，根據距離變化率平滑變化這一先驗信息，估計數據的測量方差。利用最小二乘法對R(n)進行多項式曲線擬合得到距離變化率期望的估計值r(n)，r(n)和R(n)之間的差值為測量誤差。在獨立同分布的高斯白噪聲模型下，測量誤差的方差估計值為：

若檢測窗口中目標衛星未發生變軌，則檢測窗口中的數據和參考窗口中的數據具有相同的變化趨勢，兩部分的期望值滿足相同的多項式曲線參數。因此將參考窗口內的多項式曲線進行外推，得到檢測窗口中數據檢驗統計量的估計值為：

(15)

(16)

設定虛警概率為PFA，根據(11)、(13)得到A>0與A<0兩種情況下的變軌判決門限分別為：

(17)

(18)

檢測窗口內檢驗統計量的實際值為：

(19)

將T(D)和變軌判決門限進行對比，若

則認為目標在檢測窗口內已經發生了變軌。

上述討論的是針對固定參考窗口R(n)和檢測窗口D(n)的變軌檢測問題。將參考窗口和檢測窗口的位置在序列觀測數據V(n)上進行滑動，可實現對目標衛星變軌的連續檢測，檢測流程如圖 2所示。

圖2 檢測流程Fig.2 Flow chart of detection

3 衛星軌道機動檢測仿真

3.1 構造檢測信號

本文采用STK軟件仿真衛星運行軌道[15-17]，并獲取目標衛星和伴飛衛星之間的相對距離，伴飛衛星與目標衛星的初始軌道參數設置如表 1所示，其中a,e,i,Ω,ω,f為經典軌道六要素的半長軸、偏心率、傾角、升交點赤經、近地點俯角和真近點角[18-19]。仿真衛星軌道運行時間為300 s，激光測距雷達每隔0.1 s采集一次相對距離數據。

表1 衛星初始軌道參數

當目標衛星運行到200 s時，在目標衛星運行速度方向施加脈沖小推力，使該方向速度增量為I=0.3 m/s。目標衛星和伴飛衛星之間的相對距離如圖 3所示。可見伴飛衛星逐漸靠近目標衛星，相對距離近似為線性減小；但由于衛星變軌的幅度較小，不能從相對距離曲線中觀察出目標衛星是否發生軌道機動。

圖3 相對距離示意Fig.3 Diagram of relative distance

對相對距離量進行差分處理，得到目標衛星相對伴飛衛星的距離變化率如圖4所示。可見當目標衛星發生變軌時，距離變化率有一個明顯的階躍突變特征。對比相對距離量，目標衛星軌道機動特征在距離變化率序列中的表現更加明顯，因此將后者作為特征量用于軌道機動檢測是合理的。

圖4 相對距離變化率示意Fig.4 Diagram of relative distance change rate

包含測量噪聲的距離變化率，如圖5、圖6所示。從圖中可知，由于測量噪聲的存在，使得距離變化率的階躍突變特征被淹沒在噪聲中，無法直觀地判斷目標衛星是否發生軌道機動。

圖5 未變軌時包含測量噪聲的距離變化率Fig.5 The distance change rate when the measured noise is included without orbit change

圖6 變軌時包含測量噪聲的距離變化率Fig.6 The distance change rate when the measured noise is included with orbit change

3.2 基于NP準則的軌道機動檢測

設定參考窗口的長度M=1 300，檢測窗口長度N=100。選取圖 6中一段數據進行變軌檢測，利用最小二乘法對參考窗口內的數據進行多項式曲線擬合，多項式次數設定為2，而后將曲線外推至檢測窗口，如圖7所示。

圖7 數據多項式擬合、外推結果Fig.7 The results of data polynomial fitting and extrapolation

設定虛警概率PFA=10-2，針對圖6中的數據，其檢驗統計量T(D)與判決門限的關系，如圖8所示。

圖8 檢驗統計量與門限值的對比(N=100,I=0.3 m/s)Fig.8 The comparision between test statistics and threshold values(N=100,I=0.3 m/s)

將檢測窗口長度修改為N=200，其他參數保持不變，仿真得到檢驗統計量T(D)與變軌判決門限，如9所示。

圖9 檢驗統計量與門限值的對比(N=200,I=0.3 m/s)Fig.9 The comparision between test statistics and threshold values(N=200,I=0.3 m/s)

設速度增量I=-0.3 m/s，其他參數保持不變，仿真得到檢驗統計量T(D)與變軌判決門限，如10所示。

圖10 檢驗統計量與門限值的對比(N=100,I=-0.3 m/s)Fig.10 The comparision between test statistics and threshold values(N=100,I=-0.3 m/s)

根據式(8)，在滑窗策略下將距離變化率轉化為檢驗統計量，實際上是一個濾波過程。對比濾波前后的距離變化率曲線可見，濾波有效抑制了噪聲干擾，濾波后的階躍突變特征十分明顯。圖8與圖9中設定的脈沖速度I相同，對比兩者的檢測效果可得，檢測窗口寬度較大時，距離變化率的階躍特征更容易被檢測出來,但大檢測窗口導致信號階躍突變過程變長，檢測實時性降低。圖8與圖10的變軌檢測結果表明，在加速和減速的情況下，該檢測方法均能有效檢測出目標衛星軌道機動。

3.3 機動檢測效果影響因素

不同窗口寬度參數下，階躍能量噪聲比隨著階躍幅度|A|變化的曲線如圖11所示。從圖中可知，相同的變軌幅度|A|情況下，檢測窗口寬度越大，階躍能量噪聲比越大。

圖11 幅度|A|與階躍能量噪聲比關系Fig.11 The relationship between amplitude|A| and the of step energy to noise of ratio

根據式(12)計算得到階躍能量噪聲比、虛警概率PFA與檢測概率PD的關系，如圖 12所示。

圖12 階躍能量噪聲比與檢測概率關系Fig.12 Relationship between the step energy-to-noise ratio and detection probability

圖12表明，當虛警概率PFA一定時，階躍能量噪聲比與檢測概率PD成正相關關系，即目標衛星的變軌幅度越大，測量噪聲越小，其軌道機動越容易被檢測出來。選取虛警概率PFA=10-2，當階躍能量噪聲比≥10 dB時，檢測概率PD達到0.80以上。對應在窗口寬度N=100與N=200情況下，當幅度|A|分別大于或等于0.3與0.2時，能有效檢測出目標衛星軌道機動。

4 結束語

針對脈沖小推力作用下的衛星變軌檢測問題，本文提出了基于NP概率判決模型的檢測方法，有效抑制了目標衛星測距誤差的不利影響。該方法根據伴飛衛星和目標衛星之間相對速率數據的波動特征，與設定的虛警概率，能夠自適應生成變軌判決門限，實現軌道機動檢測。通過對檢測數據作滑窗處理，可以實現持續性檢測，并估計出變軌發生的時間。本文方法根據測量數據本身的特征自適應地生成變軌判決門限，相比于之前的方法具備更強的實用性。

文中通過公式推導與仿真模擬，獲得了虛警概率、檢測概率、機動程度與測量噪聲之間的定量關系。研究表明，在給定虛警概率條件下，當目標衛星機動幅度一定時，測量誤差越小，檢測窗口寬度越大，階躍能量噪聲比會越高，對目標衛星軌道機動的檢測概率越大。變軌檢測是獲取目標衛星軌道機動特征的第一步，我們后續的研究將在檢測目標軌道機動的基礎上，進一步估計軌道機動大小、空間方向等參數。