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

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

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

0 引 言

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

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

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

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

1 理論基礎(chǔ)

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

()=((),)

(1)

(2)

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

(3)

×(,)=(×,×)

(4)

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

(5)

(6)

式(6)通過(guò)微分代數(shù)完成了導(dǎo)數(shù)的求解。當(dāng)在一開(kāi)始代入值計(jì)算時(shí)，可避免符號(hào)計(jì)算方法中復(fù)雜符號(hào)式子的存儲(chǔ)、化簡(jiǎn)和計(jì)算。高階和多變量微分代數(shù)可通過(guò)擴(kuò)展上述一元一階截?cái)嗵├照归_(kāi)表達(dá)式微分代數(shù)(式(3)～式(6))來(lái)獲得。

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

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

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

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

(7)

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

(8)

(9)

(10)

(11)

(12)

(13)

(14)

3 仿真分析

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

(15)

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

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

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

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

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

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

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

4 結(jié) 論

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