跳躍-滑翔彈道擾動引力自適應網格快速賦值方法*

王順宏，戴陳超，李 劍，楊奇松

(火箭軍工程大學 作戰(zhàn)保障學院， 陜西 西安 710025)

隨著彈道導彈制導工具誤差的不斷減小，制導方法誤差的影響日益突出，而擾動引力是影響制導方法誤差的主要因素。對于滑翔距離超過12 000 km 的高超音速滑翔導彈，擾動引力影響滑翔段彈道所產生的落點偏差達到千米量級。因此，必須提高滑翔導彈制導計算機中擾動引力賦值精度。同時，由于彈載計算機存儲空間有限，應當盡可能減小擾動引力賦值所需的數據存儲量。另外，由于是彈上實時計算，對計算速度也提出了較高的要求，傳統(tǒng)擾動引力賦值方法(如點質量法、球諧函數法等)無法同時滿足計算精度、存儲量和計算速度的要求，因此必須針對跳躍-滑翔彈道開展高精度擾動引力快速賦值方法研究。

數值逼近方法是高精度擾動引力快速賦值的主要方法，對于主動段的擾動引力，文獻[1]通過拉格朗日插值模型，利用有限元方法逼近擾動引力；文獻[2]提出了廣義延拓逼近算法在擾動引力快速賦值中的應用；文獻[3-4]分別采用了神經網絡逼近算法和三次等距B樣條函數方法，以上方法對于主動段擾動引力均取得了較好的逼近效果；文獻[5]針對彈道導彈全程提出了“漏斗形”有限元單元構建方法，將有限元方法擴展到被動段擾動引力快速賦值中，取得了較好的效果；文獻[6]采用球諧函數換極法研究被動段擾動引力快速賦值，在保證精度的前提下，提高了賦值速度，但是仍然不能滿足彈上實時計算的要求。對于跳躍-滑翔彈道的擾動引力快速賦值，上述方法均有一定的借鑒意義，但仍存在以下不足：第一，傳統(tǒng)彈道導彈的再入段通常忽略擾動引力的影響，研究重點集中于主動段與被動段，但跳躍-滑翔彈道屬于再入段，此時導彈在臨近空間內以高超音速滑翔，其彈道特性與上述文獻所研究的彈道導彈彈道特性存在較大差異。因此上文提及的快速賦值方法在跳躍-滑翔彈道的適用性需做進一步分析。第二，文獻[2]中廣義延拓逼近算法相對于一般方法精度大幅提高，但也增加了數據存儲量，同時由于跳躍-滑翔彈道的彈道軌跡較長，若按照傳統(tǒng)方法進行網格劃分同樣會導致網格數量較多，數據存儲量加大。

1 基本思路

傳統(tǒng)基于有限元法的擾動引力快速賦值是將參考彈道附近空域進行有限元劃分，建立以參考彈道為中心的飛行管道，確定各單元節(jié)點并對節(jié)點進行賦值，而后根據導彈實時位置判斷所在單元，通過插值算法，快速計算當前位置的擾動引力值。

眾所周知，對于有限元方法而言，單元劃分的大小直接決定了求解精度以及存儲的節(jié)點數目，考慮到上文提及的跳躍-滑翔彈道的特殊性，傳統(tǒng)劃分方式將導致較大的彈上數據存儲量?；谏鲜隹紤]，提出跳躍-滑翔彈道擾動引力自適應網格賦值模型，在保證賦值精度的同時，通過增加單元格邊長，從而減少單元格劃分數量，降低彈上存儲量。基本思路如下：

1)根據發(fā)射任務，在不考慮擾動引力等干擾因素下確定一條參考跳躍-滑翔彈道；

2)根據參考彈道將附近空域進行有限元劃分，以參考彈道為中心建立較大的一級網格單元，并采用球諧函數法對一級單元節(jié)點的擾動引力賦值；

3)導彈進入跳躍-滑翔段之后，根據導彈實際位置坐標確定所在一級單元，并在一級單元內部在線生成包含該點的自適應二級單元(以下簡稱二級單元)；

4)根據一級單元節(jié)點的擾動引力值，采用基于反距離加權的廣義延拓逼近算法計算二級單元節(jié)點擾動引力值，然后根據二級單元節(jié)點值，采用拉格朗日插值法快速計算導彈當前位置的擾動引力值；

5)對于新的待求點，先判斷是否屬于當前的二級單元，若在其內部，直接通過該二級單元計算擾動引力值，反之，轉第3步。

2 自適應網格構建模型

文獻[7]中基于發(fā)射坐標系對被動段彈道進行了單元劃分，該方法的優(yōu)點在于整體坐標與局部坐標的轉換較為簡單，并且提高了搜索效率。由于本文只考慮跳躍-滑翔段彈道的擾動引力計算，因此將發(fā)射坐標系平移至跳躍-滑翔段彈道下方，其原點與起始滑翔點延地心矢徑方向在地球表面的投影重合，各坐標軸方向均不變，此時稱該坐標系為再入滑翔段坐標系，簡稱再入段坐標系，本文將在此坐標系內基于基準彈道構建單元。

2.1 一級單元構建及節(jié)點確定

在再入段坐標系內設定空域Ω并進行劃分，按照直角坐標(x,y,z)截取六面體單元Ωe，Ωe可由3組坐標區(qū)間表示為{[x1,x2],[y1,y2],[z1,z2]}，同時Ωe的8個節(jié)點坐標為(xi,yi,zi)(i=1,2,3，…)，如圖1所示。

圖1 一級單元空域劃分示意圖Fig.1 Schematic diagram of the first level grid space division

(1)

通過導彈在再入段坐標系中的位置坐標即可確定其所在的一級單元。

2.2 自適應二級單元構建及節(jié)點確定

自適應二級單元的構建是根據導彈所處的一級單元以及在再入段坐標系中的坐標在線生成的，具體步驟如下：

步驟1：根據導彈的全局坐標A(xA,yA,zA)確定其所處的一級單元格j，單元格j的坐標區(qū)間可以表示為{[xj,1,xj,2],[yj,1,yj,2],[zj,1,zj,2]}，同時定義單元格j-1的坐標區(qū)間為{[xj-1,1,xj-1,2],[yj-1,1,yj-1,2],[zj-1,1,zj-1,2]}。

步驟2：判斷A(x,y,z)與一級單元坐標區(qū)間的相對位置，若x?[xj-1,1,xj-1,2]且y∈[yj-1,1,yj-1,2],z∈[zj-1,1,zj-1,2]，則定義ox軸為二級單元格擴展方向，并且當x>xj-1,2時，則沿ox軸正方向擴展，反之沿負方向擴展。若存在2個或2個以上坐標超出單元格j-1的坐標區(qū)間，則任選一個方向作為擴展方向。

步驟3：以擴展方向為ox軸方向為例，以A點為幾何中心，在平行于oyz的平面內作邊長為Δy和Δz的長方形平面k1，各邊分別與oy軸和oz軸平行。

步驟4：將步驟3中的長方形平面沿擴展方向平移Δx得到另一長方形平面k2(若k2平面超出一級單元，則直接取一級單元邊界作為k2平面)。依次連接長方形平面的各個頂點kij(i=1,2;j=1,2,3,4)，構成一個六面體單元。

步驟5：假設新待求點B的全局坐標為(xB,yB,zB)，判斷(xB,yB,zB)與當前二級單元坐標區(qū)間的相對位置。若在該二級單元內部，則以該二級單元節(jié)點進行插值計算；若不在該二級單元內部，按照步驟2確定點B所在二級單元的擴展方向，建立新的二級單元。以此類推，直至到達當前一級單元邊界。

自適應二級單元在線生成示意圖如圖2所示，圖中坐標點1～4對應頂點k11～k14，坐標點5～8 對應頂點k21～k24。

根據上述二級單元的構建方式，可確定以A(xA,yA,zA)為幾何中心，沿ox方向構建二級單元時，對應的8個節(jié)點坐標依次為：

(2)

oy軸與oz軸的表達式可類推。

(a) 自適應二級單元沿坐標軸方向擴展示意圖(a) Schematic diagram of second level adaptive grid extending on axis direction

(b) 自適應二級單元在線生成示意圖(b) Schematic diagram of second level adaptive grid online building圖2 二級單元構建示意圖Fig.2 Schematic diagram of the second level grid structure

根據式(1)，可將二級單元節(jié)點坐標轉換為局部坐標(ξ,η,ζ)，通過一級單元節(jié)點，采用基于反距離加權的廣義延拓逼近算法，插值計算出二級單元各節(jié)點的擾動引力值，再通過拉格朗日插值法計算A點擾動引力值。

3 擾動引力賦值模型

3.1 基于反距離加權的廣義延拓逼近算法

自適應網格插值計算的關鍵在于二級單元節(jié)點擾動引力的求解精度，但由于一級單元在構建過程中所選取的邊長較長、網格較大，如果采用傳統(tǒng)的拉格朗日插值方法必然導致計算結果存在較大誤差，最終影響彈道擾動引力的計算精度。文獻[7-9]將廣義延拓逼近算法用于計算彈道擾動引力，得出結論：該算法的求解精度優(yōu)于拉格朗日算法，并且廣義延拓逼近算法求解精度隨待求點位置的變化而變化，在每個單元邊緣(靠近節(jié)點)的逼近結果優(yōu)于單元中心(遠離節(jié)點)?；谏鲜鼋Y論，本文在求解二級單元節(jié)點擾動引力的過程中，根據反距離加權插值算法，引入距離權系數矩陣，提出了一種優(yōu)化廣義延拓逼近算法。在計算過程中考慮待求點與已知節(jié)點的相對位置關系，進而提高二級單元節(jié)點擾動引力值的求解精度。

3.1.1 一般廣義延拓逼近算法原理

根據延拓逼近的思想，六面體一級單元稱為主域單元，將主域單元每個節(jié)點再分別沿x,y,z方向向外延伸，形成次域單元，如圖3所示。每個主域單元有8個主節(jié)點，與之相連的次域單元有24個次域節(jié)點。

圖3 一級單元延拓示意圖Fig.3 Schematic diagram of the first level grid extension

由于擾動引力是坐標的函數，因此在一級單元內采用階次為m的多項式來逼近待求點的擾動引力。根據文獻[10]，為了避免出現嚴重的龍格現象，增強多項式的容錯性，選取m=17，并且設a為多項式系數向量，Fi為多項式類，即

a=[a0,a1,…,a16]T

(3)

(4)

因此逼近函數可以表示為：

δi=Fia

(5)

廣域延拓逼近要求主域節(jié)點上的多項式值與擾動引力計算結果相等，而次域節(jié)點上的多項式值與擾動引力差值的平方和最小，即

(6)

其中，FI=[F1,F2,F3,…,F8]T，FJ=[F9,F10,F11,…,F32]T，s為主域與次域節(jié)點上的擾動引力值，sI=[s1,s2,…,s8]，sJ=[s9,s10,…,s32]。

引入拉格朗日乘子λ=[λ1,λ2,…,λ8]，可將模型表示為：

(7)

根據優(yōu)化原理求解出a的值，通過式(5)結合待求點坐標即可求出該點的擾動引力值。

3.1.2 算法優(yōu)化及求解

次域節(jié)點的作用在于將擬合的思想加入插值計算中，使得擾動引力插值結果與次域節(jié)點多項式值之差的平方和最小。但如果待求點位于網格邊緣位置(靠近節(jié)點)，則必然存在部分次域節(jié)點與待求點的距離較遠，此時次域各個節(jié)點的重要程度應當有所區(qū)別。因此根據反距離加權插值思想，在次域節(jié)點擬合過程中引入距離權系數矩陣。

反距離加權插值算法是一種以距離作為權重的滑動平均加權插值法，其公式如下：

(8)

式中，f*表示待求點數值，f表示已知點數值，wi表示各已知點的權重。

(9)

式中，di為待求點與已知點之間的距離，k為冪指數。

根據式(9)，建立次域節(jié)點與待求點的權系數對角矩陣：

(10)

結合式(7)，優(yōu)化廣義延拓逼近算法的模型：

(11)

根據優(yōu)化原理以及分塊矩陣求逆，解出待定系數矩陣a。

(12)

式中，

(13)

結合式(5)，擾動引力為：

(14)

結合2.2節(jié)中二級單元節(jié)點坐標的求解公式，即可求出二級單元各個節(jié)點的擾動引力值。

3.2 拉格朗日插值算法

在同等單元劃分下，廣義延拓逼近算法在精度上優(yōu)于拉格朗日插值算法，但所需要的節(jié)點數目是拉格朗日插值算法的4倍?？紤]到二級單元數量較多，如果同樣采用廣義延拓法計算二級單元內待求點擾動引力，必然會導致計算量加大。因此，在二級單元內部采用模型較為簡單的拉格朗日插值算法求解。

假設二級單元內任意一點P的全局坐標為(x,y,z)(在求解過程中節(jié)點的全局坐標轉換為局部坐標(ξ,η,ζ)，為表述方便此處使用全局坐標(x,y,z))，根據2.2節(jié)可以求出二級單元各節(jié)點坐標(xm,ym,zm)(m=1,2,…,8)。采用形函數法求解P點擾動引力，令

N(x,y,z)=L(x)L(y)L(z)

(15)

式中，L(x),L(y),L(z)分別為計算x,y,z3個方向的拉格朗日插值基函數。以擴展方向為ox軸正方向為例，假設二級單元節(jié)點編號如圖2所示，則8個節(jié)點對應的形函數分別為：

(16)

設二級單元節(jié)點處擾動引力為δm(m=1,2,…,8)，則二級單元內一點P的擾動引力計算式為：

(17)

上述跳躍-滑翔彈道擾動引力賦值模型通過優(yōu)化廣義延拓逼近算法求解二級單元節(jié)點值，采用拉格朗日插值法計算自適應二級單元內部實際彈道點擾動引力值。在保證計算精度與計算速度的同時，可以減少彈上數據存儲量，適合于彈上擾動引力的實時計算。

3.3 模型誤差分析

由于在計算二級單元格節(jié)點數值時采用了最優(yōu)平方逼近的擬合思想，因而無法推算其誤差的解析表達式，因此本節(jié)采用仿真的方式對自適應網格賦值模型的誤差進行分析說明。

根據插值算法截斷誤差公式可知，當待求點位置距離插值節(jié)點越近，則插值結果的誤差越小，精度越高。因此所建立的賦值模型的誤差分析思路為：

1)分析優(yōu)化廣義延拓逼近算法賦值精度與待求點位置的關系；

2)分析拉格朗日插值算法賦值精度與待求點位置的關系；

3)結合兩種算法的誤差變化曲線，分析本文賦值模型的誤差變化曲線。

在分析算法賦值精度與待求點位置關系過程中，不同位置的待求點選取方法為：

1)在大地直角坐標系oexsyszs中的任意位置建立一個正六面體一級單元格，邊長為L。

2)在一級單元格內部建立局部坐標系oxyz，原點o為單元格中心點，ox軸、oy軸及oz軸的指向與地心大地直角坐標系一致。

3)在局部坐標系內由原點出發(fā)，按照步長Δl向外擴展構建單元格，單元格邊長記為li，新建單元格的頂點即為待求點，如圖4所示。

圖4 待求點選取示意圖Fig.4 Schematic diagram of the point selection

由于三方向擾動引力的誤差具有相似性，因此以zs軸方向擾動引力為例進行仿真計算。在仿真計算中取L=100 km，Δl=5 km，采用EGM2008模型計算一級單元節(jié)點擾動引力以及內部單元格節(jié)點擾動引力值，將同一單元格8個節(jié)點的標準誤差值作為該距離下待求點的誤差。

通過優(yōu)化廣義延拓逼近算法求解待求點擾動引力值，各內部單元格標準誤差如表1、圖5所示。

表1 優(yōu)化廣義延拓逼近算法不同位置待求點標準誤差

圖5 優(yōu)化廣義延拓逼近標準誤差與待求點位置關系Fig.5 Optimization generalized extension approximation standard error and point position

通過拉格朗日插值算法求解待求點擾動引力值，各內部單元格標準誤差如表2、圖6所示。

通過對兩種算法誤差曲線的分析，可以得出以下結論：

1)隨著待求點與一級單元格節(jié)點距離的增加，其擾動引力值的求解精度逐漸降低，誤差值逐漸增大。

2)隨著距離的增大，誤差的變化率逐漸增大。

表2 拉格朗日插值算法不同位置待求點標準誤差

圖6 拉格朗日插值標準誤差與待求點位置關系Fig.6 Lagrange interpolation standard error and point position

本文的賦值模型采取了兩級單元格分步逼近的方式，根據圖5及圖6兩種算法的誤差曲線可以近似采用如下公式表示：

(18)

其中：Δδi(i=1,2)表示兩種算法的標準誤差；ai,bi(i=1,2)表示兩種算法誤差的待定系數；xl表示待求點與其所在一級單元格最近節(jié)點的距離。

假設某一一級單元格I1內存在待求點A，其在一級單元格I1內對應的距離xl記為xA，因此當采用優(yōu)化廣義延拓逼近算法求解時，其標準誤差為：

(19)

假設待求點A的拓展二級單元為I2，且單元格I2的中心點與一級單元格重合，此時二級單元格8節(jié)點在一級單元格I1內對應的距離xl相等，記為x1,待求點A在二級單元格I2內對應的距離xl記為x2。

當采用自適應網格賦值模型求解時，其標準誤差為：

(20)

將式(19)與式(20)相減得：

(21)

式中由于誤差計算采用的是標準誤差，因此b2>0。當x1≥x2時，

(22)

由表1及表2可知，優(yōu)化廣義延拓逼近的誤差小于拉格朗日逼近，結合圖5及圖6可知，a1-a2<0，且a1與a2的差距較小，因此3a1-a2>0，R>0，即自適應網格賦值模型的誤差小于優(yōu)化廣義延拓逼近。

4 仿真分析

4.1 基本仿真條件

一級單元劃分以參考彈道所在再入段坐標系空域為基礎進行劃分，一級單元以及自適應二級單元均劃分為正六面體。一級單元節(jié)點擾動引力采用球諧函數法進行賦值，采用EGM2008模型2156階位系數[11]。參考跳躍-滑翔彈道以遠程高超聲速滑翔式再入飛行器為對象，氣動參數采用美國波音公司1998年設計的帶控制翼的錐形體再入機動飛行器CAV-H的參數擬合得到[12]。

4.2 賦值精度分析

設定滑翔彈道的起始高度為70 km，起始速度為6500 m/s，為了使得仿真結果更具一般性，取步長為45°，將再入方位角由0°遞增至315°，選此8個再入方位角的典型跳躍-滑翔彈道進行仿真計算[13]，將8條彈道對應的誤差結果取平均值作為最終結果進行精度分析?？紤]到彈道末端擾動引力對落點精度影響很小，因此以高度20 km作為滑翔彈道終端條件，一級單元劃分選擇相鄰節(jié)點間Δx=Δy=Δz=100 km，同時二級單元選定Δx=Δy=Δz=20 km。圖7給出了正北方向彈道擾動引力逼近效果(由于篇幅有限，僅給出x軸方向逼近效果圖)，參考值為采用球諧函數法計算得到的擾動引力值，所采用的數據與節(jié)點賦值數據相同。圖8給出了不同逼近方法標準差對比。表3給出了在同等單元格劃分下，不同賦值方法對三方向擾動引力的逼近精度。

圖7 擾動引力逼近效果圖Fig.7 Accuracy of gravity anomaly approximating

圖8 不同賦值方法標準誤差對比Fig.8 Comparison of standard error of different assignment methods

通過分析上述圖表的仿真結果可得：

1)由圖7的逼近效果可以看出，所建模型的賦值精度較高，誤差不隨時間的增加而積累，整體逼近效果較好。

2)由圖8可知，優(yōu)化廣義延拓逼近的標準誤差小于傳統(tǒng)廣義延拓逼近，優(yōu)化廣義延拓自適應網格的標準誤差小于廣義延拓自適應網格，因此基于反距離加權的優(yōu)化廣義延拓算法進一步提高了傳統(tǒng)廣義延拓逼近算法的求解精度，有更好的逼近效果。

3)自適應網格插值與傳統(tǒng)網格插值相比較，自適應網格的標準誤差沒有因為二次插值而增加，同時在加入優(yōu)化廣義延拓之后，其標準誤差與傳統(tǒng)廣義延拓自適應網格相比有所下降，逼近精度提高。

表3 不同賦值方法逼近精度比較

考慮到擾動引力計算的主要目的是進行擾動引力落點影響修正，因此給出了優(yōu)化廣義延拓自適應網格的逼近誤差對應落點偏差的大小，其橫向偏差為1.225 4 m，縱向偏差為4.338 2 m[14]。

4.3 彈上存儲量分析

彈上存儲量的大小與單元格劃分大小直接相關，在保證逼近精度的前提下，單元格劃分越大，則所需要的單元格數量越少，因此彈上需要存儲的數據量越少；反之，單元格劃分越小，則所需單元格數量越多，從而彈上存儲量越大。

將網格劃分為正六面體，邊長取L，針對一般廣義延拓逼近法和優(yōu)化廣義延拓自適應網格模型，通過選取不同的網格劃分大小，分別分析兩種賦值方法的逼近精度。表4為四種網格大小下賦值算法的絕對誤差最大值以及標準誤差。圖9表示三方向標準誤差隨網格大小的變化情況，其中初始邊長L0=100 km，步長ΔL=50 km，最大網格邊長Lf=500 km。

表4 不同單元格大小下算法的逼近精度比較

(a) x

(b) y

(c) z圖9 標準誤差隨網格大小變化情況Fig.9 Variation of standard error with the size of the grid

由上述圖表可知：

1)兩種賦值方法的逼近誤差隨著網格邊長的增加而增大，但自適應網格模型的標準誤差的變化率總體上低于一般廣義延拓逼近方法。

2)當標準誤差相等時，自適應網格賦值方法對應的單元格邊長較長，圖9(c)中，當標準誤差為10×10-5m/s2時，一般廣義延拓的單元格邊長約為325 km，而自適應網格賦值的單元格邊長約為450 km，單元格邊長增幅超過38%。因此在同等精度要求下，與一般賦值方法相比，優(yōu)化廣義延拓自適應網格賦值方法的最大網格邊長更長，可以減少參考彈道的網格劃分數量，從而減小彈上數據存儲量。

4.4 適應性分析

自適應網格插值方法依賴于基準彈道，基準彈道的形狀可能會影響單元的構建，進而對擾動引力賦值精度造成影響。因此有必要分析在不同滑翔距離、滑翔方向下自適應網格賦值方法的適應性。

通過改變起始速度，分別針對滑翔距離為9000 km、13 000 km以及16 000 km的三條跳躍滑翔彈道進行分析。圖10表示起始滑翔方向為-180°～180°情況下三條彈道擾動引力逼近誤差對應落點偏差的大小。

圖10 不同滑翔距離彈道擾動引力逼近 誤差對應落點偏差Fig.10 Drop point deviation corresponding to the gravity anomaly approximation error under different gliding distance

由圖10可知，滑翔距離越遠，擾動引力逼近誤差對應落點偏差越大，但均保持在5 m之內。由此可知，在各個起始滑翔方向下，擾動引力逼近精度均很高，可以滿足導彈精度要求，自適應網格模型具有較好的適應性。

4.5 賦值速度分析

擾動引力的賦值包括兩個方面：地面數據準備階段和彈上實時賦值階段。在普通配置的微機上(CPU主頻2.53 GHz，內存512 MB)，對于4.2節(jié)的跳躍-滑翔彈道，地面數據準備時間為16.8 s，彈上單點賦值時間為0.19 ms。因此自適應網格模型雖然增加了計算量，但計算時間仍然滿足快速諸元計算以及彈上實時計算需要。

5 結論

本文提出了一種自適應網格賦值模型，并根據反距離加權理論，優(yōu)化了廣義延拓逼近算法。通過與傳統(tǒng)賦值方法的對比分析以及適應性分析，得出以下結論：

1)在同等大小的網格劃分下，優(yōu)化廣義延拓自適應網格模型的逼近精度高于一般賦值方法。

2)隨著單元格邊長的增加，優(yōu)化廣義延拓自適應網格賦值模型的誤差增加平緩，在同等精度要求下，該賦值模型的最大單元格邊長大于一般賦值方法，從而減小了單元格劃分數量，降低了彈上數據存儲量。

3)該賦值模型能夠適應不同滑翔方向以及不同滑翔距離的彈道擾動引力快速賦值，逼近誤差對應落點偏差小于5 m。同時模型的計算速度滿足快速諸元計算以及彈上實時計算需要。