周歡， 丁智堅(jiān)， 鄭偉

(1.中國工程物理研究院 總體工程研究所， 四川 綿陽 621999； 2.中國空氣動力研究與發(fā)展中心， 四川 綿陽 621000；3.國防科技大學(xué) 航天科學(xué)與工程學(xué)院， 湖南 長沙 410073)

0 引言

地球外部空間擾動引力作用于飛行過程中的導(dǎo)彈，影響其導(dǎo)航解算和制導(dǎo)計算，進(jìn)而影響落點(diǎn)精度，故需要在發(fā)射前構(gòu)建引力模型?，F(xiàn)有點(diǎn)質(zhì)量法、球諧函數(shù)、Stokes方法等空間擾動引力計算方法難以滿足諸元計算及飛行過程中擾動引力賦值的存儲量和計算速度要求，由此需要建立快速逼近理論?？焖俦平碚摰暮诵氖墙⑿问胶唵?、計算效率高、存儲量小和適應(yīng)范圍廣的數(shù)值函數(shù)，既能精確地反映擾動引力空間特征，又能滿足彈道計算要求。其中涉及到空間結(jié)構(gòu)剖分、重構(gòu)函數(shù)確定、效能損失控制以及重構(gòu)精度與速度優(yōu)化抉擇等問題。

國外，點(diǎn)質(zhì)量法[1]曾應(yīng)用于民兵Ⅲ導(dǎo)彈的重力場計算中。為提高解算速度，一些學(xué)者探討了譜方法[2-5]和數(shù)值逼近方法在擾動引力求解中的應(yīng)用。Junkins[6]首次提出了重力位的有限元表達(dá)，并與Engels等[7]探討了其在慣性導(dǎo)航計算中的應(yīng)用。國內(nèi)，廣義延拓逼近[8]、內(nèi)插[9-10]、三次等距B樣條函數(shù)[11]等方法被用于求解主動段擾動引力，球諧函數(shù)換極法[12]被用于求解被動段擾動引力。然而，上述方法或難以滿足彈上存儲量和計算效率的要求，或因無法兼顧數(shù)據(jù)的空間趨勢信息而提高逼近精度。

針對以上問題，本文提出了一種基于網(wǎng)函數(shù)逼近理論的沿飛行彈道擾動引力快速重構(gòu)方法，推導(dǎo)了該逼近方法的余項(xiàng)誤差，分析了影響賦值精度的因素，驗(yàn)證了該方法在各種情況下的適應(yīng)性。

1 沿飛行彈道的空域剖分

基于有限元思想，按如下思路構(gòu)建擾動引力快速計算模型：1)確定不考慮擾動引力的參考彈道；2)形成以參考彈道為中心的飛行管道，進(jìn)行空域剖分并確定各節(jié)點(diǎn)位置；3)通過高精度引力計算方法進(jìn)行節(jié)點(diǎn)擾動引力賦值；4)判斷計算點(diǎn)所在單元及其與該單元各節(jié)點(diǎn)的相對位置關(guān)系；5)根據(jù)一定的逼近方法，由節(jié)點(diǎn)擾動引力計算當(dāng)前點(diǎn)擾動引力。擾動引力對落點(diǎn)精度的影響隨飛行時間累積，擬采用“漏斗形”管道以減少單元數(shù)和防止實(shí)際彈道超出管道。

以主動段飛行管道的構(gòu)建為例。利用主動段彈道相對平直的特性，在發(fā)射坐標(biāo)系內(nèi)生成管道，進(jìn)行空域剖分：1)在標(biāo)準(zhǔn)彈道上依次選定基準(zhǔn)點(diǎn)d0，d1，d2，d3，…，di，相鄰兩點(diǎn)間沿y軸方向的距離為δy1,δy2,δy3,…，δyi；2)以di為幾何中心，在平行于Oxz的平面內(nèi)生成邊長分別為δxi和δzi的長方形；3)令δxi<δxi+1，δyi<δyi+1，δzi<δzi+1，依次連接各長方形頂點(diǎn)ki，構(gòu)成主動段“漏斗形”飛行管道。采用類似方法在軌道坐標(biāo)系內(nèi)構(gòu)建被動段飛行管道。

2 擾動引力快速逼近方法

快速逼近方法是高精度擾動引力快速賦值的主要方法，擾動引力的快速逼近主要涉及兩個問題：一是判斷計算點(diǎn)所在單元；二是單元內(nèi)部的逼近計算。前者通過比較計算點(diǎn)與di的相對位置關(guān)系獲得，后者基于網(wǎng)函數(shù)逼近理論實(shí)現(xiàn)。

由Cook[13]提出的網(wǎng)函數(shù)逼近理論是一種多方向擬合一次誤差調(diào)整的多變元函數(shù)逼近方法，其本質(zhì)上是單變元函數(shù)的Lagrange插值算法的極自然推廣。其插值函數(shù)類是簡單的高階偏微分方程邊值問題的解，具有Coons型結(jié)構(gòu)以及明顯的統(tǒng)計特征，算法簡單。該方法的基本思想是通過n維歐式空間中網(wǎng)格邊界上值來近似計算網(wǎng)格中任意點(diǎn)的值，與單元內(nèi)部的擾動引力逼近問題契合。

令L(x)、L(y)、L(z)分別為關(guān)于x、y、z的一次Lagrange插值算符，其插值基函數(shù)為

(1)

記主動段某六面體單元的8個頂點(diǎn)為ki,1、ki,2、ki,3、ki,4、ki+1,1、ki+1,2、ki+1,3、ki+1,4，其擾動引力值分別為gi,1、gi,2、gi,3、gi,4、gi+1,1、gi+1,2、gi+1,3、gi+1,4；記12條棱為L0～L11，稱其為計算單元上的1-網(wǎng)，定義在1-網(wǎng)上的擾動引力值為fi(x,y,z)(i=0,1,…,11)。根據(jù)網(wǎng)函數(shù)逼近理論，已知1-網(wǎng)函數(shù)，可由三維1-網(wǎng)插值方法求取單元內(nèi)任意一點(diǎn)的值。令l(3)為三維1-網(wǎng)插值算符，則

l(3)=L(x)L(y)+L(y)L(z)+L(z)L(x)-
2L(x)L(y)L(z).

(2)

據(jù)此，可求得單元內(nèi)任意一點(diǎn)A(x,y,z)的值F(x,y,z)，

F(x,y,z)=F1(x,y,z)+F2(x,y,z)+
F3(x,y,z)-F4(x,y,z)，

(3)

式中：

(4)

節(jié)點(diǎn)擾動引力可由傳統(tǒng)擾動引力計算方法獲得，而fi(x,y,z)(i=0,1,…,7)則通過對每條棱所對應(yīng)的兩節(jié)點(diǎn)插值獲得。對于沿彈道方向的棱fi(x,y,z)(i=8,9,10,11)，為保證飛行管道的平滑，納入相鄰網(wǎng)格數(shù)據(jù)的空間趨勢信息；同時為避免因多點(diǎn)插值引起的龍格現(xiàn)象，采用加權(quán)3點(diǎn)插值的方式求取。以f9(x,y,z)的求解為例，

(5)

通過相同方法求取其他3條沿彈道方向的棱上的擾動引力值。至此，已完成了主動段單元內(nèi)部任意一點(diǎn)擾動引力的求解，通過類似方法可完成被動段擾動引力的快速逼近。

3 逼近精度分析

3.1 插值余項(xiàng)及誤差估計

稱R(3)[F]=F(x,y,z)-l(3)[fi(x,y,z)]為函數(shù)F(x,y,z)的l(3)插值余項(xiàng)，記γ(x,y,z)=R(3)[F]。當(dāng)F(x,y,z)有直到6階的連續(xù)偏導(dǎo)數(shù)，則可推導(dǎo)γ(x,y,z)的估計式，

(6)

式中：

(7)

(8)

η1、ζ1為Oyz平面上的給定點(diǎn)坐標(biāo)，ζ2、ξ2為Oxz平面上的給定點(diǎn)坐標(biāo)，ξ3、η3為Oxy平面上的給定點(diǎn)坐標(biāo)，ξ0、η0、ζ0為Oxyz空間坐標(biāo)系下的給定點(diǎn)坐標(biāo)。

由(6)式可知，γ(x,y,z)取決于計算單元大小以及表征F(x,y,z)變化率的高階偏導(dǎo)數(shù)。因此，單元越小，數(shù)據(jù)變化越平緩，計算點(diǎn)離節(jié)點(diǎn)越近，則逼近程度越高。

3.2 逼近精度影響因素分析

本節(jié)通過數(shù)值仿真，對3.1節(jié)的結(jié)論進(jìn)行驗(yàn)證。以分層點(diǎn)質(zhì)量法進(jìn)行節(jié)點(diǎn)擾動引力賦值。針對某一固定空域進(jìn)行等大小的網(wǎng)格剖分，以網(wǎng)函數(shù)方法和分層點(diǎn)質(zhì)量法分別計算空域內(nèi)1 000個隨機(jī)點(diǎn)的擾動引力，二者之差即為逼近誤差。

3.2.1 改變網(wǎng)格大小

計算原點(diǎn)坐標(biāo)(λ,φ,H)=(110°,40°,0 km)，其中λ為經(jīng)度，φ為緯度，H為高程；計算范圍(Δλ,Δφ,ΔH)=(±1.5°,±1.5°,300 km)。選取3種大小的網(wǎng)格：

1)尺寸1網(wǎng)格(dλ,dφ,dH)=(0.25°,0.25°,5 km)；

2)尺寸2網(wǎng)格(dλ,dφ,dH)=(0.5°,0.5°,10 km)；

3)尺寸3網(wǎng)格(dλ,dφ,dH)=(0.75°,0.75°,15 km)。

擾動引力三分量逼近誤差的統(tǒng)計結(jié)果見表1.

表1 不同網(wǎng)格大小的逼近誤差統(tǒng)計

3.2.2 改變網(wǎng)格形狀

計算原點(diǎn)坐標(biāo)與3.2.1節(jié)相同，選取體積大致相等、形狀不同的4種網(wǎng)格：

1)形狀1網(wǎng)格(dλ,dφ,dH)=(0.1°,0.1°,10 km)；

2)形狀2網(wǎng)格(dλ,dφ,dH)=(0.2°,0.2°,2.5 km)；

3)形狀3網(wǎng)格(dλ,dφ,dH)=(0.05°,0.05°,40 km)；

4)形狀4網(wǎng)格(dλ,dφ,dH)=(0.2°,0.1°,20 km)。

4種網(wǎng)格形狀的示意圖見圖1. 圖1中，R為沿地球矢徑方向，L為沿經(jīng)度方向，B為沿緯度方向。擾動引力三分量逼近誤差的統(tǒng)計結(jié)果見表2.

表2 不同網(wǎng)格形狀的賦值誤差統(tǒng)計

3.2.3 改變數(shù)據(jù)空間變化趨勢

基于同一地區(qū)擾動引力的變化隨高程增加而趨于平緩的結(jié)論，設(shè)計數(shù)值仿真算例。計算范圍及計算原點(diǎn)經(jīng)度和緯度同3.2.2節(jié)，計算原點(diǎn)高程分別為：1)H=0 km；2)H=20 km；3)H=40 km；4)H=60 km. 網(wǎng)格大小為(dλ,dφ,dH)=(0.5°,0.5°,10 km)。擾動引力三分量逼近誤差的統(tǒng)計結(jié)果見表3.

表3 數(shù)據(jù)變化趨勢對賦值誤差統(tǒng)計結(jié)果的影響

由表2的計算結(jié)果可以看出，網(wǎng)格形狀越趨近正方體，逼近誤差越小。由表3的計算結(jié)果可以看出，隨高程增加，擾動引力數(shù)據(jù)趨于平緩，逼近誤差逐漸減小。對比表1和表2的結(jié)果可知，形狀2～形狀4的網(wǎng)格遠(yuǎn)小于尺寸1～尺寸3的網(wǎng)格，但逼近精度并沒有得到明顯提高，甚至低于尺寸1的賦值精度。由此可知，網(wǎng)格形狀比網(wǎng)格大小對逼近精度具有更大的影響。由表1～表3數(shù)值仿真得到的結(jié)論與3.1節(jié)中理論推導(dǎo)結(jié)論一致。

4 快速逼近方法的適應(yīng)性分析

由于地球外部空間擾動引力場受地球表面形狀及地球內(nèi)部物質(zhì)分布的影響，通過將快速逼近方法應(yīng)用于不同發(fā)射區(qū)域、不同射向和不同射程的彈道導(dǎo)彈全程彈道計算中檢驗(yàn)方法的適應(yīng)性。適應(yīng)性判斷準(zhǔn)則如下：

1)擾動引力計算模型適應(yīng)于任意彈道，全程彈道積分不出現(xiàn)奇異；

2)實(shí)際飛行彈道不超出預(yù)先構(gòu)建的“漏斗形”彈道管道逼近空間；

3)根據(jù)彈道導(dǎo)彈射前引力模型構(gòu)建的一般要求，擾動引力單向分量的逼近誤差均值應(yīng)不大于1 mgal；

4)根據(jù)彈道導(dǎo)彈對落點(diǎn)精度的一般要求，由擾動引力二次逼近誤差導(dǎo)致的落點(diǎn)偏差應(yīng)不超過100 m.

4.1 單彈頭全程擾動引力快速賦值結(jié)果分析

設(shè)置主動段起始單元邊長為0.5 km，單元邊長增幅為2 km. 被動段起始單元邊長為δr1=δξ1=200 km，δf1=2°，單元邊長增幅為20 km和0.5°. 節(jié)點(diǎn)擾動引力及近似真值均采用1 080階球諧函數(shù)進(jìn)行計算。

4.1.1 不同發(fā)射區(qū)域

針對射程為12 000 km，發(fā)射方位角為90°的彈道開展研究。發(fā)射點(diǎn)分別為：1)平原λ=115°，φ=35°，H=0 km；2)丘陵λ=110°，φ=27°，H=1 km；3)特大山區(qū)λ=80°，φ=42°，H=3 km. 賦值結(jié)果見表4.

表4 不同發(fā)射區(qū)域的賦值誤差統(tǒng)計結(jié)果

4.1.2 不同發(fā)射方位角

針對射程為12 000 km，發(fā)射點(diǎn)位于特大山區(qū)的彈道開展研究。發(fā)射方位角α分別為：1)正北α=0°；2)正東α=90°；3)正南α=180°；4)正西α=270°. 賦值結(jié)果見表5.

表5 不同發(fā)射方位角的賦值誤差統(tǒng)計結(jié)果

4.1.3 不同射程

針對發(fā)射點(diǎn)位于特大山區(qū)，發(fā)射方位角α=90°，射程分別為5 000 km、8 000 km和12 000 km的彈道開展研究。賦值結(jié)果見表6.

表6 不同射程彈道的賦值誤差統(tǒng)計結(jié)果

4.1節(jié)的仿真結(jié)果表明，本文建立的擾動引力快速逼近方法能夠適應(yīng)各種條件下的彈道計算，全程彈道積分無奇異且不存在實(shí)際飛行彈道超出預(yù)設(shè)逼近空間的情況。表4～表6的統(tǒng)計結(jié)果表明，在不同發(fā)射區(qū)域、不同射向和不同射程的彈道導(dǎo)彈全程彈道計算中，擾動引力三向分量的逼近誤差均值控制在10-2mgal量級，滿足單向分量逼近誤差均值不超過1 mgal的要求。

4.2 雙彈頭全程擾動引力快速賦值結(jié)果分析

本節(jié)將擾動引力快速逼近方法應(yīng)用于多頭分導(dǎo)彈道導(dǎo)彈，以考察方法的適應(yīng)性。假定導(dǎo)彈射向?yàn)檎?，發(fā)射點(diǎn)的經(jīng)度、緯度和高程分別為111°、39°和1 500 m. 假設(shè)分導(dǎo)彈頭的射程小于主彈頭，且主彈頭落點(diǎn)和分導(dǎo)彈頭落點(diǎn)縱程相差約1 000 km，橫程相差約160 km. 節(jié)點(diǎn)擾動引力采用1 080階球諧函數(shù)法賦值，并以球諧函數(shù)法計算結(jié)果作為計算參考值。考慮到分導(dǎo)過程發(fā)生在主動段結(jié)束以后，分導(dǎo)彈頭與主彈頭的彈道僅在被動段不同，因此只針對被動段擾動引力賦值情況進(jìn)行分析。被動段空域剖分參數(shù)設(shè)置如下：1)起始網(wǎng)格大小δz0=δr0=200 km，δf0=2.0°；2)相鄰網(wǎng)格邊長增幅δzi+1-δzi=δri+1-δri=20 km，δfi+1-δfi=0.5°. 經(jīng)仿真可以得到，主彈頭和分導(dǎo)彈頭的被動段彈道及其彈道管道，如圖2所示。

上述空域剖分方案對應(yīng)的網(wǎng)格總數(shù)為25個，節(jié)點(diǎn)數(shù)為88個，總存儲量為528個數(shù)據(jù)，可以滿足彈上數(shù)據(jù)存儲量要求。根據(jù)上述單元劃分，將快速賦值方法的計算結(jié)果與球諧函數(shù)方法計算結(jié)果進(jìn)行比較，可以得到主彈頭及分導(dǎo)彈頭被動段擾動引力逼近誤差，如圖3和圖4所示。

由圖3和圖4可知，在被動段飛行的大部分時段，主彈頭和分導(dǎo)彈頭擾動引力逼近誤差均可控制在1 mgal以內(nèi)，具有較高的賦值精度。僅當(dāng)飛行至被動段末端，由于漏斗形彈道管道形成的逼近空間逐漸增大，導(dǎo)致賦值精度下降，但逼近誤差仍控制在10 mgal以內(nèi)。隨著飛行時間增加，擾動引力賦值誤差對落點(diǎn)精度的影響逐漸降低，上述賦值精度滿足應(yīng)用需求。

在主彈頭和分導(dǎo)彈頭的落點(diǎn)偏差計算中，分別基于球諧函數(shù)方法和快速賦值方法求解擾動引力項(xiàng)。表7給出了兩種擾動計算方法對應(yīng)的彈頭落點(diǎn)縱向偏差、橫向偏差和總偏差。以球諧函數(shù)的計算結(jié)果為參考值，視兩種方法計算結(jié)果之差為快速賦值方法的逼近誤差對應(yīng)的落點(diǎn)偏差。結(jié)果表明，快速賦值方法逼近誤差對應(yīng)的落點(diǎn)偏差可控制在米級，具有較高精度。

表7 擾動引力計算誤差對應(yīng)的落點(diǎn)偏差

4.3 逼近誤差對落點(diǎn)精度的影響分析

針對發(fā)射點(diǎn)為λ=110°，φ=40°，H=1.5 km，射程分別為5 000 km、8 000 km和12 000 km，發(fā)射方位角為αi=-180°+i×5°(i=0,1,2,…,72)的彈道開展研究。彈道積分中的擾動引力項(xiàng)分別采用快速逼近方法和1 080階球諧函數(shù)方法計算得到，視兩種彈道計算結(jié)果對應(yīng)的落點(diǎn)偏差之差為賦值誤差引起的落點(diǎn)偏差，計算結(jié)果如圖5所示。結(jié)果表明，各種情況下賦值誤差引起的落點(diǎn)偏差均不超過8 m，滿足實(shí)際應(yīng)用要求。

綜合4.1節(jié)～4.3節(jié)的仿真結(jié)果，認(rèn)為本文建立的擾動引力快速逼近方法滿足逼近算法適應(yīng)性判斷準(zhǔn)則，能夠適應(yīng)于不同情況下的彈道計算。

5 討論

5.1 計算時間及存儲量

本文提出的擾動引力快速逼近方法主要包括引力模型構(gòu)建及飛行過程中的賦值計算，其中前者主要包括空域剖分、節(jié)點(diǎn)坐標(biāo)確定及節(jié)點(diǎn)擾動引力賦值，后者主要包括當(dāng)前計算點(diǎn)所在單元判斷及單元內(nèi)部的逼近計算。存儲數(shù)據(jù)主要包括各節(jié)點(diǎn)位置三分量及各節(jié)點(diǎn)擾動引力三分量。根據(jù)第4.1節(jié)的網(wǎng)格劃分條件，表8列出了3種射程彈道的模型構(gòu)建速度及存儲量分析結(jié)果，表9列出了不同擾動引力計算方法下射程為12 000 km彈道的計算時間。仿真計算所使用的計算機(jī)基本配置為：Celeron(R)CPU 2.53 GHz，內(nèi)存大小為512 MB. 軟件環(huán)境為Window XP操作系統(tǒng)，計算程序基于VC++6.0開發(fā)。結(jié)果表明，本文提出的方法在計算時間和存儲量方面具有很大的優(yōu)勢。

表8 重構(gòu)速度及存儲量分析

表9 不同計算方法對應(yīng)的計算時間

5.2 賦值精度及適應(yīng)性

針對某射程為12 000 km，射向?yàn)檎狈较虻膹椀?，在不考慮節(jié)點(diǎn)擾動引力賦值誤差情況下，基于網(wǎng)函數(shù)方法和型函數(shù)方法分別計算沿不同高度彈道的擾動引力。與第4節(jié)相同，仍以1 080階球諧函數(shù)方法的擾動引力計算結(jié)果作為近似真值，將網(wǎng)函數(shù)與型函數(shù)的計算結(jié)果分別與之比較，可以得到各自的賦值誤差均方差，如表10所示。

表10 網(wǎng)函數(shù)和型函數(shù)賦值誤差統(tǒng)計結(jié)果

由表10可以看出，網(wǎng)函數(shù)的賦值精度略高。從理論上講，網(wǎng)函數(shù)與型函數(shù)都是Lagrange插值的推廣，因此具有量級相當(dāng)?shù)谋平?。但與型函數(shù)等內(nèi)插計算方法相比，基于網(wǎng)函數(shù)逼近理論的賦值方法除利用當(dāng)前單元節(jié)點(diǎn)的離散數(shù)據(jù)外，還考慮了相鄰單元節(jié)點(diǎn)數(shù)據(jù)的空間趨勢信息，根據(jù)引力異常變化趨勢構(gòu)造最佳的1-網(wǎng)逼近曲線，進(jìn)而以逼近曲線為邊界曲線來調(diào)控整個曲面的形狀和趨向，構(gòu)建分塊表征模型。這一算法相當(dāng)于在離散的節(jié)點(diǎn)數(shù)據(jù)基礎(chǔ)上追加了一定的邊界條件，因此提高了精度。由于采用的數(shù)據(jù)都是已有節(jié)點(diǎn)的數(shù)據(jù)，因此無需增加額外的存儲量。

該方法在不同發(fā)射區(qū)域、不同射向和不同射程彈道中的應(yīng)用結(jié)果表明，其賦值誤差在10-2mgal量級，由此引起的落點(diǎn)偏差在8 m以內(nèi)，具有精度高、計算速度快、存儲量小和適應(yīng)范圍廣等特點(diǎn)。

6 結(jié)論

本文基于標(biāo)準(zhǔn)彈道構(gòu)建飛行管道，實(shí)現(xiàn)對沿飛行彈道附近空域的有限元剖分，首次將網(wǎng)函數(shù)逼近理論應(yīng)用于擾動引力快速賦值計算，提出了一種沿飛行彈道的擾動引力快速計算方法。在理論推導(dǎo)該方法的插值余項(xiàng)誤差基礎(chǔ)上，基于某一區(qū)域的擾動引力測量數(shù)據(jù)，數(shù)值模擬了不同因素對賦值精度的影響。以1 080階球諧函數(shù)構(gòu)建的地球引力模型為近似真值，分析了提出的快速計算方法在不同情況下彈道導(dǎo)彈全程彈道計算中的應(yīng)用情況。結(jié)果表明，該方法產(chǎn)生的賦值誤差約為10-2mgal量級，由此引起的落點(diǎn)偏差在8 m以內(nèi)，全程彈道計算時間遠(yuǎn)小于其他方法。因此，本文提出的方法能夠?qū)崿F(xiàn)沿飛行彈道任意點(diǎn)擾動引力三分量的快速賦值，其賦值精度、計算速度和存儲量均滿足彈道計算要求，具有廣泛的適應(yīng)性。