一種基于微分代數(shù)的任意階空間目標(biāo)軌道傳播方法及其分析*

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

(上海機(jī)電工程研究所·上海·201109)

0 引 言

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

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

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

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

1 理論基礎(chǔ)

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

()=((),)

(1)

(2)

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

(3)

×(,)=(×,×)

(4)

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

(5)

(6)

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

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

2 空間目標(biāo)軌道任意階導(dǎo)數(shù)求解

2.1 空間目標(biāo)動力學(xué)模型

空間目標(biāo)軌道傳播數(shù)值積分方法通常采用空間目標(biāo)動力學(xué)方程的考埃爾形式，這是由于該形式可以很容易地添加任意攝動加速度模型。考埃爾形式如下

(7)

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

(8)

(9)

(10)

(11)

(12)

(13)

(14)

3 仿真分析

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

(15)

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

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

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

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

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

圖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)給出空間目標(biāo)狀態(tài)對時間的任意階導(dǎo)數(shù)，進(jìn)而利用泰勒展開式進(jìn)行軌道傳播。本文在二體模型下進(jìn)行仿真分析，給出了隨微分代數(shù)階次的變化情況和滿足<1時的最大軌道外推時間隨微分代數(shù)階次的變化情況，并驗證了方法的有效性。本文還給出了龍格-庫塔方法和微分代數(shù)方法的外推精度比較，初步討論了算法性能。