高 凱，潘金波，賈世偉，張宇琛，沈俊逸

(上海機電工程研究所·上?！?01109)

0 引 言

彈道導彈飛行的中段將釋放彈頭、干擾機和誘餌；在再入段，彈頭和誘餌將高速再入。雷達探測彈頭時，經(jīng)過多目標關(guān)聯(lián)識別和干擾機對抗，真實彈頭回波依然可能會淹沒在噪聲下。這時需要對弱信號進行長時間積累，以期達到早發(fā)現(xiàn)的目的。弱信號的長時間積累依賴精確的目標動力學模型及其運動狀態(tài)外推。動力學模型通常用非線性微分方程來表示。而目標運動狀態(tài)外推由于動力學模型的復雜性一般都沒有解析解，故外推通常采用數(shù)值積分方法。隨著現(xiàn)代計算機的發(fā)展，數(shù)值積分方法精確求解高階非線性微分方程變得工程可行。

數(shù)值積分方法通常分為單步方法和多步方法。單步方法利用初始狀態(tài)求解不同時刻不同運動狀態(tài)下微分方程的值，組合這些值即可獲得單步目標軌道傳播，典型的方法為龍格-庫塔(Runge-Kutta，RK)方法。多步方法通過組合目標的多個后向狀態(tài)(即除了初始狀態(tài)，還需要初始狀態(tài)時刻之前時刻的狀態(tài))實現(xiàn)目標軌道傳播，常用方法有亞當斯-巴什福斯-莫爾頓(Adams-Bashforth-Moulton) 方法和尚平-戈登(Shampine-Gordon)方法。單步方法相比多步方法不需要目標狀態(tài)的后向值。多步方法通常采用單步方法來起步，但要注意階次的匹配。很多情況下，四階龍格-庫塔方法(RK4)的精度是足夠的，但在需要高階多攝動模型時，則需要更高階的數(shù)值積分方法，以保證積分的精度、削弱步進選擇帶來的精度和計算時間復雜度的影響。高階龍格-庫塔方法需要求解復雜的條件方程以獲得參數(shù)表，使其在需要高階場景下的應用變得復雜。

微分代數(shù)(Differential Algebra，DA)提供了一種使用現(xiàn)代計算機求解函數(shù)導數(shù)的方法。最早涉足微分代數(shù)方法的是Joseph Liouville，他將函數(shù)積分和有限項微分方程組聯(lián)系在一起。之后Joseph Fels Ritt給出了求微分方程組解的完整代數(shù)理論，使微分代數(shù)方法獲得了顯著的進步。Ellis Robert Kolchin和Robert Henry Risch在算法方面推進了微分代數(shù)方法的發(fā)展。當前另外兩種計算機求導方法分別是有限差分方法和符號計算方法。有限差分方法的缺點主要是差分間隔很難取得合適和存在無法準確獲得導數(shù)值的固有缺陷。在高維空間和高階偏導求解時，有限差分方法的時間復雜度非常高。符號計算求導方法可以和微分代數(shù)方法一樣獲得準確的導數(shù)值，但是該方法在處理復雜表達式求導時會難以化簡表達式，導致計算機內(nèi)存緊張，并且時間復雜度變得極高。這兩種方法的固有特點使得它們難以在大規(guī)模復雜表達式求導過程中應用。本文將重點關(guān)注微分代數(shù)技術(shù)在求解微分方程和偏微分方程中的應用，特別是初始條件已知時流微分方程的泰勒展開。

本文給出了一種基于微分代數(shù)的任意階空間目標軌道傳播方法。該方法不需要改變外推算法的計算過程，并避免了求解復雜條件方程。文中仿真分析采用空間目標的二體運動模型(該模型可求解析解，方便分析比較)，對軌道傳播方法進行分析。本文分析了步進時間和階次對空間軌道傳播精度的影響，并對高階微分代數(shù)方法和龍格-庫塔方法的精度進行了比較。

1 理論基礎(chǔ)

微分代數(shù)通過將實數(shù)代數(shù)變量替換為微分代數(shù)變量，使函數(shù)最終計算獲得的微分代數(shù)變量值包含函數(shù)對自變量的各階微分。微分代數(shù)支持所有實數(shù)代數(shù)的運算。對如下函數(shù)關(guān)系

()=((),)

(1)

(2)

(,)+(,)=(+,+)

(3)

×(,)=(×,×)

(4)

(,)·(,)=(·,·+·)

(5)

(6)

式(6)通過微分代數(shù)完成了導數(shù)的求解。當在一開始代入值計算時，可避免符號計算方法中復雜符號式子的存儲、化簡和計算。高階和多變量微分代數(shù)可通過擴展上述一元一階截斷泰勒展開表達式微分代數(shù)(式(3)～式(6))來獲得。

微分代數(shù)的詳細計算機代碼實現(xiàn)主要有源代碼轉(zhuǎn)換方法、操作符重載方法和表達式模板方法。源代碼轉(zhuǎn)換方法首先使用計算機代碼實現(xiàn)目標函數(shù)，然后用預處理器按照微分法則生成新的求函數(shù)導數(shù)的代碼。源代碼轉(zhuǎn)換方法僅能利用在代碼編譯時已知的信息，難以處理像循環(huán)、C++ 模板和其他面向?qū)ο筇匦浴２僮鞣剌d的核心思想是引入新的類對象，表達微分鏈式法則上的微分中間量，其缺點是難以完整地表示中間變量。而類模板方法則可以解決這一問題。類模板可以返回表達式模板，在編譯階段直接通過編譯器獲得計算導數(shù)鏈式法則層次結(jié)構(gòu)，并對中間變量自動形成合適的類對象。

2 空間目標軌道任意階導數(shù)求解

2.1 空間目標動力學模型

空間目標軌道傳播數(shù)值積分方法通常采用空間目標動力學方程的考埃爾形式，這是由于該形式可以很容易地添加任意攝動加速度模型?？及栃问饺缦?/p>

(7)

2.2 任意階導數(shù)求解

(8)

(9)

(10)

(11)

(12)

(13)

(14)

3 仿真分析

采用數(shù)值積分方法和二體運動解析解比較的方法進行仿真。將二體運動解析解作為空間目標二體運動的精確解，便于向數(shù)值積分方法提供精度參考。仿真采用的空間目標在地心慣性系(Earth-Centered Inertial，ECI)中的初始狀態(tài)如表 1所示。采用精度衰減因子(Position Dilution of Precision，PDOP)來量化和比較數(shù)值積分方法的解和二體運動解析解之間的差異。

(15)

式中，(=,,)為時刻數(shù)值積分方法獲得位置矢量與二體運動解析解位置矢量在方向上的差值。

表1 空間目標在地心慣性系下的狀態(tài)

首先,運用微分代數(shù)方法對空間目標初始狀態(tài)進行外推，外推時間為100.0s。仿真時依次獲得更高階導數(shù)，求積分后獲得外推狀態(tài)，之后與二體運動解析解進行比較。圖1給出了隨微分代數(shù)階次的變化情況，可以看出,隨微分代數(shù)階次近似以指數(shù)形式衰減，這符合從泰勒展開式中獲得的預期結(jié)果。圖1中，以分貝形式給出(=10×log())。需要指出的是，在求解高階導數(shù)時,需要采用高精度數(shù)據(jù)類型,以應對函數(shù)導數(shù)大的動態(tài)范圍。

圖1 σPDOP隨微分代數(shù)階次的變化Fig.1 σPDOP changes with the order of differential algebra

圖2給出了滿足<1時的最大軌道外推時間隨微分代數(shù)階次的變化情況。首先，通過仿真求解最大軌道外推時間對應階次的導數(shù)，然后將上一階次求得的最大軌道外推時間乘以2(為迭代次數(shù))，直至≥1，記錄此時的和-1對應的時刻，在和-1對應的時刻間隔內(nèi)用二分法查找=1的值。由圖2可以看出，滿足<1的最大軌道外推時間隨微分代數(shù)階次近似以指數(shù)形式增長，這符合從泰勒展開式中獲得的預期結(jié)果。圖3給出了分別用龍格-庫塔方法和本文所提的微分代數(shù)方法對二體運動數(shù)值積分與其解析解的比較曲線，可以看出，高階方法在大時間步進時可有效提高外推精度。

圖2 滿足σPDOP<1的最大軌道外推時間隨微分代數(shù)階次的變化Fig.2 Changes of the maximum orbit propagation time that satisfies σPDOP<1 with the order of differential algebra

圖3 微分代數(shù)方法與龍格-庫塔方法的比較Fig.3 Comparison of differential algebraic method and Runge-Kutta method

4 結(jié) 論

本文給出了一種基于微分代數(shù)的解決軌道傳播問題的單步數(shù)值積分方法。該方法可在攝動力模型解析區(qū)間內(nèi)給出空間目標狀態(tài)對時間的任意階導數(shù)，進而利用泰勒展開式進行軌道傳播。本文在二體模型下進行仿真分析，給出了隨微分代數(shù)階次的變化情況和滿足<1時的最大軌道外推時間隨微分代數(shù)階次的變化情況，并驗證了方法的有效性。本文還給出了龍格-庫塔方法和微分代數(shù)方法的外推精度比較，初步討論了算法性能。