周歡， 丁智堅， 鄭偉

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

0 引言

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

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

針對以上問題，本文提出了一種基于網函數逼近理論的沿飛行彈道擾動引力快速重構方法，推導了該逼近方法的余項誤差，分析了影響賦值精度的因素，驗證了該方法在各種情況下的適應性。

1 沿飛行彈道的空域剖分

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

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

2 擾動引力快速逼近方法

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

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

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

(1)

記主動段某六面體單元的8個頂點為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-網，定義在1-網上的擾動引力值為fi(x,y,z)(i=0,1,…,11)。根據網函數逼近理論，已知1-網函數，可由三維1-網插值方法求取單元內任意一點的值。令l(3)為三維1-網插值算符，則

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

(2)

據此，可求得單元內任意一點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é)點擾動引力可由傳統擾動引力計算方法獲得，而fi(x,y,z)(i=0,1,…,7)則通過對每條棱所對應的兩節(jié)點插值獲得。對于沿彈道方向的棱fi(x,y,z)(i=8,9,10,11)，為保證飛行管道的平滑，納入相鄰網格數據的空間趨勢信息；同時為避免因多點插值引起的龍格現象，采用加權3點插值的方式求取。以f9(x,y,z)的求解為例，

(5)

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

3 逼近精度分析

3.1 插值余項及誤差估計

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

(6)

式中：

(7)

(8)

η1、ζ1為Oyz平面上的給定點坐標，ζ2、ξ2為Oxz平面上的給定點坐標，ξ3、η3為Oxy平面上的給定點坐標，ξ0、η0、ζ0為Oxyz空間坐標系下的給定點坐標。

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

3.2 逼近精度影響因素分析

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

3.2.1 改變網格大小

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

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

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

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

擾動引力三分量逼近誤差的統計結果見表1.

表1 不同網格大小的逼近誤差統計

3.2.2 改變網格形狀

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

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

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

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

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

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

表2 不同網格形狀的賦值誤差統計

3.2.3 改變數據空間變化趨勢

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

表3 數據變化趨勢對賦值誤差統計結果的影響

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

4 快速逼近方法的適應性分析

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

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

2)實際飛行彈道不超出預先構建的“漏斗形”彈道管道逼近空間；

3)根據彈道導彈射前引力模型構建的一般要求，擾動引力單向分量的逼近誤差均值應不大于1 mgal；

4)根據彈道導彈對落點精度的一般要求，由擾動引力二次逼近誤差導致的落點偏差應不超過100 m.

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

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

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

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

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

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

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

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

4.1.3 不同射程

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

表6 不同射程彈道的賦值誤差統計結果

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

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

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

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

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

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

表7 擾動引力計算誤差對應的落點偏差

4.3 逼近誤差對落點精度的影響分析

針對發(fā)射點為λ=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)的彈道開展研究。彈道積分中的擾動引力項分別采用快速逼近方法和1 080階球諧函數方法計算得到，視兩種彈道計算結果對應的落點偏差之差為賦值誤差引起的落點偏差，計算結果如圖5所示。結果表明，各種情況下賦值誤差引起的落點偏差均不超過8 m，滿足實際應用要求。

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

5 討論

5.1 計算時間及存儲量

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

表8 重構速度及存儲量分析

表9 不同計算方法對應的計算時間

5.2 賦值精度及適應性

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

表10 網函數和型函數賦值誤差統計結果

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

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

6 結論

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