徐曉東 趙建亭 許春雷

(江蘇自動化研究所，江蘇 連云港 222006)

0 引言

射程在數千千米的遠程彈道導彈在沿地球外部上空飛行時，時刻受到地球引力場的作用。由于高速飛行的彈道導彈彈載計算機處理參數的速度和時間非常有限，因此選擇計算擾動引力的方法，除要滿足一定的計算精度外，還要盡可能的簡便，以減少計算量、提高計算速度［1］。由于斯托克斯積分方法和點質量法能夠滿足導彈制導與導航的精度，在彈道導彈擾動引力計算的工程中廣泛應用。但是這兩種方法具有計算量大、不能滿足快速計算要求的不足;而與梯度法、球諧函數展開法相比較，它們又具有計算精度高和計算范圍不受限制的優(yōu)點，適用于全程彈道導彈彈道擾動引力的計算［2］。實際應用中一般都采用先在地面計算標準彈道附近點的擾動引力，然后通過標準彈道擾動引力插值或逼近來計算實際彈道的擾動引力［3］。

在滿足工程需要精度要求的前提下，為了達到彈道導彈擾動引力快速計算的要求，本文在采用斯托克斯積分法計算地球外部空間擾動引力的基礎上，提出了簡化后的快速計算模型，并根據簡化模型的特點，運用并行計算技術來提高計算速度、減少計算時間。仿真試驗驗證了快速計算模型的有效性。

1 斯托克斯積分原理

采用斯托克斯積分方法計算地球外部空間擾動引力，就是將用斯托克斯積分公式表示的地球外部空間擾動引力位對任一方向求偏導數，以求其擾動引力在該方向上的分量。

假設地球為圓球體時，采用斯托克斯積分公式表示的擾動引力位公式為:

式中:R為地球的平均半徑;r為球心到P點的距離;ρ為球面上面積圓dσ到P點的距離;ψ為面積圓dσ與P點間對應的地心極角;φs、λs分別為P點的地心緯度和經度、分別為面積圓dσ的地心緯度和經度;函數S(r，ψ)為廣義斯托克斯函數;Δgσ為球面重力異常值［4－5］。

地球擾動引力位示意圖如圖1所示。

圖1 地球擾動引力位示意圖Fig.1 Sketch map of the earth gravitational potential of disturbance

2 快速計算模型

2.1 地心彈體坐標系

假設地球是圓球體，坐標原點取在地心，X軸為地心與彈體中心連線的延長線，Y軸在彈道導彈P點的速度V方向與地心構造的平面內，且與X軸垂直，Z軸與X、Y軸構成右手直角坐標系。由此構成的地心彈體坐標系示意圖如圖2所示。

圖2 地心彈體坐標系Fig.2 The coordinate system of the geocentric elastomer

2.2 斯托克斯積分簡化模型

利用簡化的斯托克斯積分法求解外部空間擾動引力模型時，首先將地球以適當的經緯度差(根據實際探測的地面重力異常值區(qū)域)進行劃分。在劃分的基礎上，在地心彈體坐標系中，沿著X軸逆向觀察，若X軸正方向經過的球體點恰好在某一經緯度差區(qū)域內，則以此區(qū)域作為起始區(qū)域;若X軸正方向經過的球體點恰好在幾個經緯度差區(qū)域邊界上，則以這幾個相鄰的區(qū)域作為起始區(qū)域，定義經過的輪次為第一輪，經緯度差區(qū)域的個數為N1。接著以X軸為軸心，沿著Y軸和Z軸向外八個方向進行延伸，除了起始區(qū)域外，首先經過的經緯度差區(qū)域定義為第二經緯度區(qū)域，定義輪次為第二輪，經緯度差區(qū)域個數為N2。然后繼續(xù)向外八個方向延伸，除了起始區(qū)域和第二經緯度區(qū)域外，先經過的經緯度差區(qū)域定義為第三經緯度區(qū)域，經緯度差區(qū)域個數為N3。依次類推，直至當繼續(xù)外推時出現與上一輪次經緯度差區(qū)域相重疊的區(qū)域時，將上一輪次經過的經緯度差區(qū)域定義為結束區(qū)域，經過的輪次定義為M，經緯度差區(qū)域個數為NM。

將上述斯托克斯積分法求解外部空間擾動引力積分方程離散化為:

式中:下標r，e，n分別代表天北東坐標系中的三個坐標方向分量。

球面上的經緯度差區(qū)域的面積為:

式中:φi為該區(qū)域中心的緯度;φ為該區(qū)域的緯度、為取自第i經緯度區(qū)域中隨機選擇一個經緯度差區(qū)域中心的地心經度和地心緯度(由于假設地球為圓球，則可以認為這些經緯度差區(qū)域的r，R，ρ，ψ相等);Δλ、Δφ分別為每一個經緯度差區(qū)域四邊的經緯度差。

3 并行計算技術

在采用斯托克斯積分法求解彈道導彈外部空間擾動引力的過程中，隨著經緯度差的減少，求解的精度會增加，但同時計算工作量以指數倍增加，存在計算工作量大的問題。為了解決這一問題，本文采用并行計算技術來提高計算的速度，以滿足實時性的要求。

主從模式是一種比較常用的并行計算架構模式，它有一個控制進程，稱為主進程，其余的進程稱為從進程。主從架構模式的體系結構如圖3所示。

圖3 主從模式并行機體系結構Fig.3 Architecture of the master-slave parallel machine

負載均衡問題是影響并行效率的主要因素。在本文的任務分配中，采用地心彈體坐標系左右半球平均分配的原則，使每個處理器中的計算量基本相同，大體能夠滿足負載均衡的要求。在由三個CPU構成的主從模式并行機構中，程序設計結構框圖如圖4所示。

圖4 程序設計結構框圖Fig.4 Structural block diagram of program design

MPI是一個用于開發(fā)基于消息傳遞并行程序的標準，它提供了一個實際可用、可移植、高效和靈活的消息傳遞接口庫。MPI支持非阻塞通信的方式，即 MPI_Isend(buf，count，datatype，dest，tag，comm，request)能夠實現計算和通信的重疊［6］，節(jié)省了傳輸數據的時間。同時，MPI還支持基于非阻塞通信模式的廣播方式和收集方式。試驗證明，與普通的發(fā)送(MPI_Send)和接收(MPI_Recv)操作相比，采用廣播和收集的方式收發(fā)數據更有效。

4 仿真計算

本文采用5°×5°的經緯度差劃分方法，將整個球體劃分成(360/5)×(360/5)塊經緯度差區(qū)域，各個經緯度差區(qū)域球面上的平均重力異常值Δg(i，j)σ通過數據模擬實現，地球外部空間任一點 P(r，φs，λs)簡化的取值為P(7 378 140，0，0)，R 取值為地球的平均半徑6 371 km，則地球外部空間與球平面的交點為(6 371 000，0，0)。以此交點為中心，沿地心彈體坐標系X軸向周圍八個方向進行延伸，這八個經緯度差區(qū)域即是第二經緯度區(qū)域。以此類推，可以確定其他的經緯度區(qū)域。傳統(tǒng)計算模型與快速計算方法結果比較如表1所示。

表1 傳統(tǒng)/快速計算法的比較Tab.1 Comparison of the fast calculation method and traditional calculation method

表1計算結果表明，簡化的快速計算模型與傳統(tǒng)的計算方法相比，計算時間降低了33.52%，由此證明了快速計算簡化模型的有效性。

為了進一步降低簡化的快速計算模型的計算時間，本文引入了如圖4所示的并行計算技術。本文在輸入相同數據的情況下，比較并行與串行程序的執(zhí)行時間，并計算相應的加速比和效率。串行/并行計算技術比較結果如表2所示。

表2 串行/并行計算技術的比較Tab.2 Comparison of serial/parallel computation technology

上述計算結果表明:并行計算方法能夠獲得較高的加速比和效率，具有比較好的快速計算效果。

5 結束語

遠程彈道導彈飛行過程中時刻受到擾動引力的影響，而計算外部空間擾動引力的斯托克斯積分方法本身存在計算量大、模型復雜的不足，不能滿足實時計算的要求。針對上述情況，本文提出了斯托克斯解算擾動引力的快速計算模型。通過仿真可知，在滿足精度要求的情況下，該模型可以大幅度提高計算速度。同時，針對斯托克斯快速計算模型在計算擾動引力的過程中各個經緯度差區(qū)域獨立性強的特點，在Linux環(huán)境下應用MPI技術實現了斯托克斯積分法的并行計算。仿真驗證并行效果明顯，能得到較好的加速比和并行效率，從而有效減少擾動引力的計算時間，提高計算的效率。彈道導彈在飛行過程中速度很快，對擾動引力的計算速度提出了很高的要求，研究基于插值和補償的外部空間擾動引力的快速計算問題，以滿足實際應用的需要，將有待于進一步的探討。

［1］王慶賓，周世昌，王世忠，等.彈道主動段全射向擾動引力快速逼近方法［J］.測繪科學技術學報，2010，4(2):79 －81.

［2］張毅，肖龍旭，王順宏，等.彈道導彈彈道學［M］.2版.長沙:國防科技大學出版社，1999:73－76.

［3］王繼平，王明海，陳摩西.彈道導彈主動段擾動引力的一種逼近算法［J］.航天控制，2008(3):30 －38.

［4］ Difrancesco D，Meyer T，Christensen A，et al.Gravity gradiometry:today and tomorrow［C］//SAGA Biennial Technical Meeting and Exhibition，Swaziland，2009:80 －83.

［5］張赤軍，駱鳴津，王新勝，等.地球內外擾動物質引起高程異常的分析［J］.大地測量與地球動力學，2010，12(6):42 －45.

［6］都志輝.高性能計算并行編程技術-MPI并行程序設計［M］.北京:清華大學出版社，2001:99－124.