尤 浩 王 寧 趙久奮 彭 泉

（1.火箭軍工程大學 西安 710025）（2.中國國際工程咨詢有限公司 北京 100080）

1 引言

隨著我國經(jīng)濟的日益發(fā)展，人口越來越稠密，地面的各種建筑物和設施不斷增加，加上作戰(zhàn)理論的豐富發(fā)展以及武器裝備數(shù)量的增多，多導彈協(xié)同作戰(zhàn)的方式已經(jīng)逐漸成為目前戰(zhàn)爭的主要方式，多導彈協(xié)同作戰(zhàn)帶來的直接影響就是導彈殘骸數(shù)量增多，相應地面規(guī)避區(qū)域人員、物資撤離難度加大、范圍廣，由此，對導彈殘骸落區(qū)安全問題的研究日益必要［1］。

保障導彈殘骸落區(qū)安全可以從將殘骸落點控制在人煙稀少的地方和提高落區(qū)劃定的精度，減小落點散布的范圍這兩方面入手［2］。

目前，關(guān)于殘骸落區(qū)劃定的文獻較多，但對具體定量計算殘骸落區(qū)和地面規(guī)避區(qū)域之間安全性的文獻則相對較少。文獻［3］通過坐標變換關(guān)系得出的運動學方程，它只討論導彈運動的幾何性質(zhì)，不涉及產(chǎn)生運動和改變運動的原因，算出殘骸落區(qū)精度較差。文獻［4］則是以牛頓力學定律為基礎的積分法，沒有考慮空氣阻力的影響，殘骸落區(qū)覆蓋范圍較大。文獻［5］對某發(fā)射衛(wèi)星中心的歷次火箭殘骸落點數(shù)據(jù)進行了統(tǒng)計分析，基于分級預警思想，對不同預警區(qū)的危險系數(shù)進行了計算，但是該方法定義的預警區(qū)域范圍較廣，相應地面規(guī)避區(qū)域人員撤離工作量較大。文獻［6］分析了殘骸在下落過程中的動力學特性，提出了帶空氣阻力修正的殘骸落點算法，較好地縮小了導彈殘骸落區(qū)分布范圍。文獻［7］通過引進柵格翼氣動性能分析法，大幅縮小了殘骸落點的散布范圍。

通過上述文獻可知，現(xiàn)階段，導彈殘骸落區(qū)精度已經(jīng)有了較大水平的提高，較大地減少了規(guī)避區(qū)域人員資產(chǎn)撤離的損失，但是目前還需要對所有在殘骸落區(qū)范圍內(nèi)的規(guī)避區(qū)域進行人員物資撤離，當參與導彈數(shù)量很多時，損失任然很大。

針對上述問題，本文詳細提出了一種導彈殘骸落區(qū)和地面規(guī)避區(qū)域的位置關(guān)系判斷的方法，在此基礎上，通過構(gòu)建二維正態(tài)分布模型定量快速計算導彈殘骸落區(qū)多邊形和規(guī)避區(qū)域多邊形的交集區(qū)域概率，并將交集區(qū)域概率與設置閾值比較，從而選擇性地選擇相關(guān)規(guī)避區(qū)域人員和物資撤離，極大地減小了人力物力損失。

2 坐標的地球投影變換

2.1 坐標系的空間轉(zhuǎn)換

假設地球是半徑為R的一個空間球體，構(gòu)建空間直角坐標系如圖1所示：圓點位于空間圓球的中心，Z軸指向空間圓球的北極，X軸指向起始子午面與赤道的交點，Y軸位于赤道面上且按右手系與X軸呈90°夾角。某點在空間中的坐標可用該點在此坐標系的各個軸上的投影來表示［8］。

圖1 空間直角坐標系示意圖

設某點的經(jīng)度為L，緯度為B，則該點的地理坐標和空間直角坐標的轉(zhuǎn)換關(guān)系為

在投影變換前，需要將導彈殘骸落區(qū)多邊形和規(guī)避區(qū)域多邊形的地理坐標系內(nèi)的經(jīng)緯度數(shù)據(jù)統(tǒng)一采用上述方法轉(zhuǎn)換到三維空間直角坐標系。

2.2 殘骸落區(qū)切平面方程的建立

在圖1建立的空間直角坐標系內(nèi)，已知球面上一點A=(x0，y0，z0)的坐標，則該點的切平面方程可用平面點法式方程表示，點法式的法線為

投影面的選擇影響計算精度，只要在導彈殘骸落區(qū)內(nèi)任意選擇一點做切平面，計算的誤差相差不大。本文采用輸入殘骸落區(qū)多邊形的第一個點做投影切平面，比較方便簡單。

2.3 空間直角坐標系到切平面的投影

做直角坐標系上任意一點B=(xB，yB，zB)到切平面的垂線，則切平面上垂直投影點為B'，則：

投影點B'坐標計算如下：

2.4 切平面內(nèi)的二維坐標系建立

投影變換的目的是把三維數(shù)據(jù)變換到二維平面內(nèi)，便于計算，因此本文在投影的切平面內(nèi)擬建立一個二維坐標系。二維坐標系的原點設為切點A=(x0，y0，z0)，設地球北極點在平面內(nèi)的投影點為N=(xN，yN，zN)，取方向為X軸方向。不妨設在A點存在一個從球內(nèi)指向球外的Z軸，則可由右手定則確定Y軸的方向。

Y軸方向：

對于切平面上任意一點B'=(xB'，yB'，zB')，其在二維坐標系下的坐標可表示為

通過上述幾步計算即可以將地理坐標系以經(jīng)緯度表示的點，投影到一個二維坐標系下的平面。由此，可快速計算出導彈殘骸落區(qū)多邊形和地面各規(guī)避區(qū)域多邊形在投影平面上的二維坐標點。

3 快速判斷兩個多邊形的交集區(qū)域

求解兩個多邊形的交集區(qū)域?qū)嶋H意義在于求解殘骸落區(qū)與規(guī)避區(qū)域在二維坐標系內(nèi)的相交區(qū)域［7～8］。

如圖2所示，平面多邊形1與平面多邊形2相交，交集區(qū)域多邊形表示為k1b2k4k3k2e1?？梢园l(fā)現(xiàn)交集區(qū)域多邊形都是由不同折線段首尾拼接而成的。

通過觀察可知，交集區(qū)域多邊形的每條折線段都是由相交的點和相包含的點組成。由此，本文將快速判斷兩個多邊形的交集區(qū)域分解為以下三個過程：圖形預處理，求取交集折線段，拼接交集折線段。

圖2 兩多邊形交集

3.1 圖形預處理

本文交集區(qū)域折線段的求取采用兩個多邊形互相遍歷邊的方式。對于一般的多邊形，每條邊上的兩個點關(guān)于另一個多邊形存在以下四種關(guān)系。

1）起點和終點都在多邊形外。如圖2中的b1c1、c2d2、c1d1等；

2）起點在多邊形外，終點在多邊形內(nèi)。如圖2中的a2b2、d1f1等；

3）起點在多邊形內(nèi)，終點在多邊形外。如圖2中的f1a1、b2c2等；

4）起點和終點都在在多邊形內(nèi)。如圖3中的g2h2、b2c2等。

圖3 兩多邊形交集

上述四種位置關(guān)系，包含了所有邊與另一個多邊形的位置關(guān)系。然而這四種關(guān)系又可以細分為多種情況，容易造成復雜多邊形計算不準確，主要原因是一條邊與另一個多邊形的交點個數(shù)大于1，理論上可以根據(jù)不同“關(guān)系”，設計由交點“奇偶位置”選出組成交集多邊形折線段規(guī)則，這樣需要對四種關(guān)系分別設計規(guī)則，造成算法顯得冗長復雜，下面本文將給出一種簡單高效的處理方法。

本文采用“圖形預處理”的方法解決該問題，具體預處理方法如下：如果一條邊與另一個多邊形交點個數(shù)大于1，則在每兩個相鄰交點中點處添加一個“虛擬的頂點”。則圖2經(jīng)過預處理后兩個多邊形為a1b1c1v1d1f1和a2b2c2d2v2。圖3經(jīng)過預處理得到新的多邊形分別為a1b1c1v1d1e1v3v4和a2b2c2d2e2v2f2g2h2i2。加入的虛擬頂點用vi表示。通過上述預處理得到新的多邊形，保證每條邊與另一個多邊形交點個數(shù)最大為1。預處理后的交集區(qū)域如圖4所示。

圖4 預處理后交集區(qū)域

如圖4所示，經(jīng)過預處理的兩個多邊形，如果存在交集區(qū)域，那么其中一個多邊形至少必有一個頂點在另一個多邊形內(nèi)部。判斷兩個多邊形頂點的包含關(guān)系都需要互相判斷一遍。

3.2 求取交集折線段

經(jīng)過預處理后，再來討論前文中的四種關(guān)系。1）起點和終點都在多邊形外。該線段全在多邊形外，不做處理。

2）起點在多邊形外，終點在多邊形內(nèi)。這種關(guān)系表示遍歷方向從多邊形外進入多邊形內(nèi)，生成一條新的“折線段”保存后續(xù)值。

3）起點在多邊形內(nèi)，終點在多邊形外。這種關(guān)系表示遍歷方向從多邊形內(nèi)走出到多邊形外，結(jié)束一條“折線段”的記錄。

4）起點和終點都在多邊形內(nèi)。該條邊在多邊形內(nèi)，將起點和終點加入已經(jīng)存在的折線段中。

求取圖形折線段過程如下。

輸入：經(jīng)預處理后的平面二維點表示的兩個多邊形區(qū)域Q1，Q2。

1）按照某個順序依次遍歷多邊形Q1中的每個頂點v；

2）以v點為終點，v點的前一頂點為起點，判斷起點和終點是否在多邊形Q2中。

3）參考上述四種關(guān)系。如果起點在Q2外，終點在Q2外，該線段全在Q2外，不做處理；如果起點在Q2外，終點在Q2內(nèi)，遍歷方向從多邊形外進入多邊形內(nèi)，生成一條新的“折線段”，依次保存“進入交點”和v點；如果起點在Q2內(nèi)，終點在Q2外，遍歷方向從多邊形內(nèi)穿到多邊形外，依次保存v點和“穿出交點”，保存折線段，結(jié)束該條折線段；起點在多邊形內(nèi)，終點在多邊形內(nèi)，該條邊在多邊形內(nèi)，將起點和終點加入已經(jīng)存在的折線段中。

4）結(jié)束多邊形Q1的頂點遍歷；

5）采用類似方法，依次遍歷Q2中的所有頂點，計算各折線段；

6）輸出組成兩個多邊形的交集區(qū)域的所有折線段。

3.3 拼接交集折線段

在“求取交集折線段”部分已經(jīng)計算出了組成交集區(qū)域的所有折線段，下一步就是將所有折線段內(nèi)的點按照先后順序拼接成一個多邊形。

方法如下。

輸入：由上面計算得到的多條待拼接折線段組成的集合。

1）從折線段集合中任意選取一條折線段作為初始“拼接交集折線段”，取該條折線段的一個端點為起點，則另一個端點為終點；

2）用“拼接交集折線段”的終點與集合中剩余折線段的起點和終點比較，如果與某條折線段的起點重合，則正向拼接；如果與某條折線段的終點重合，則反向拼接，形成新的“拼接交集折線段”；

3）判斷“拼接交集折線段”的起點和終點是否相等，即判斷拼接的多邊形是否已經(jīng)閉合，如果閉合，保存該多邊形，轉(zhuǎn)下一步；否則，轉(zhuǎn)2）；

4）判斷集合是否為空，即所有的折線段是否都已經(jīng)拼接，如果為空，轉(zhuǎn)下一步；如果不為空，轉(zhuǎn)1）繼續(xù)新的多邊形拼接；

5）輸出所有拼接完成的多邊形。

至此完成了兩個多邊形交集區(qū)域的算法設計。當殘骸落區(qū)多邊形確定之后，只需要將規(guī)避區(qū)域多邊形與其進行上述計算即可。當規(guī)避區(qū)域是多個時，只需將每一個多邊形與其計算即可［9］。

4 殘骸落區(qū)和規(guī)避重合區(qū)域概率計算

利用計算機計算二維正態(tài)分布的概率涉及到復雜的二重積分計算，而且隨著計算區(qū)域的復雜多變，這種積分計算更加復雜耗時［10］。積分計算的一個優(yōu)點是精度非常高，然而大多數(shù)情況我們并不關(guān)心小數(shù)點第四位后的情況，因此本文采用一種近似的計算方法［11～12］，滿足精度要求的同時使得算法更加簡單，計算速度更快。

4.1 二維正態(tài)分布密度公式

4.2 二維正態(tài)分布概率計算

對二重積分區(qū)域D內(nèi)的概率為

本文采用的近似計算方法采用上述的思想，設想將積分區(qū)域D細分成無數(shù)個小網(wǎng)格，每個小網(wǎng)格的面積為gridArea，在每個小網(wǎng)格內(nèi)選擇一個特征點，計算該特征點對應的概率密度f(xi，yi)，則小網(wǎng)格對應概率Pgrid≈gridArea·f(xi，yi) ，小網(wǎng)格面積gridArea相當于長方體的底面積，f(xi，yi)相當于長方體的高，只要網(wǎng)格分得足夠細，有。上述方法的理論基礎是將做二元級數(shù)展開后，只考慮低階項，忽略高階項。

下面將詳細討論本文采用的近似方法，設μ1=0，σ12=1，μ2=0，σ22=1，ρ=0為標準二維正態(tài)分布，密度計算公式為fs(x，y)。參考圖5，計算的步驟如下。

輸入：殘骸落區(qū)多邊形區(qū)域C、殘骸落區(qū)與邊區(qū)域的交集多邊形區(qū)域D

1）根據(jù)輸入的殘骸落區(qū)多邊形C求解其在X軸和Y軸上的最大最小值，建立一個最小外包矩形M，左下點為(xmin，ymin)，右上點為 (xmax，ymax)。

3）建立到標準分布的轉(zhuǎn)換方程，將上述平面上的值轉(zhuǎn)換到標準二維正態(tài)分布中得：

對于標準分布，由于σ=1，所以上式最后的σ都略去。上式使得xmin對應，xmax對應；ymin對應，ymax對應。對于原平面二維坐標系內(nèi)的任意一點都可以轉(zhuǎn)換到標準二維正態(tài)分布坐標系。

4）將最小外包矩形M在X軸和Y軸方向上分別N等分，于是共有N×N個方格，每個方格的，面積為dxdy；在標準二維正態(tài)坐標系中，。

5）取每個方格的中點(xci，yci)為該方格的特征點，將其轉(zhuǎn)換到標準坐標系下：(xci_standard，yci_standard) ，如果該中點(xci，yci) 在在積分區(qū)域D內(nèi)（并非殘骸落區(qū)C，而是殘骸落區(qū)與規(guī)避區(qū)域所有的交集區(qū)域為積分區(qū)域D），計算該方格區(qū)域的概率為fs(xci_standard，yci_standard)dxstandarddystandard。

6）疊加所有計算的方格概率，即為積分區(qū)域D的概率。同理，如果將積分區(qū)域D劃分為若干個小區(qū)域Di，則可以按上述方法計算出每個區(qū)域Di的概率。

至此，完成關(guān)于導彈殘骸落區(qū)和相應規(guī)避區(qū)域重合區(qū)域的二維正態(tài)分布的概率計算。

圖5 平面二維系和標準二維正態(tài)分布坐標系

5 算例分析

程序執(zhí)行環(huán)境為MatlabR2014a，設計安全閾值為book=68,ebook=73。以經(jīng)度122，緯度26為導彈殘骸落區(qū)分布中心，協(xié)方差矩陣為，隨機生成一組導彈殘骸落區(qū)二維高斯分布經(jīng)緯度坐標點，隨機生成規(guī)避區(qū)域的有序頂點經(jīng)緯度坐標，使得規(guī)避區(qū)域在殘骸落區(qū)內(nèi)外都有分布，仿真結(jié)果如圖6所示；通過上述坐標轉(zhuǎn)換和地球投影方法較快算出了殘骸落區(qū)多邊形C和各規(guī)避區(qū)域多邊形Gi在投影平面上的二維坐標，而后快速判斷交集區(qū)域范圍Ai，取N=1000，仿真結(jié)果如圖7所示；并在構(gòu)建的二維正態(tài)分布模型中定量快速計算殘骸落區(qū)多邊形和規(guī)避區(qū)域多邊形的交集區(qū)域概率Pi，仿真分布圖如圖8所示，最終計算結(jié)果如表1所示。

圖6 殘骸落區(qū)和各規(guī)避區(qū)域經(jīng)緯度分布

圖7 殘骸落區(qū)和各規(guī)避區(qū)域在平面二維系分布

圖8 殘骸落區(qū)在平面二維系下標準正態(tài)分布

表1 交集區(qū)域計算概率

通過上述算例計算結(jié)果可知，規(guī)避區(qū)域G3對應的交集區(qū)域A3的概率要大于閾值Y，即導彈殘骸落入該規(guī)避交集區(qū)域概率較大，因此，為安全起見，需要對該區(qū)域的人員物資進行轉(zhuǎn)移；同理，因?qū)棜埡÷淙胧Ｏ乱?guī)避區(qū)域概率極小，低于設定閾值，即該規(guī)避區(qū)域交集區(qū)域是安全的，相應區(qū)域人員物資則無需轉(zhuǎn)移。相比較傳統(tǒng)的轉(zhuǎn)移全部規(guī)避區(qū)域人員物資而言，本文提出的算法不僅能夠快速定量評估出殘骸落區(qū)安全性，還極大地減小了實際過程中的轉(zhuǎn)移工作量。

6 結(jié)語

本文對導彈殘骸落區(qū)范圍和地面規(guī)避區(qū)域交集區(qū)域安全性進行了研究。首先，推導了地球投影變換的相關(guān)算法過程，建立了從地理坐標系到切平面二維坐標系的轉(zhuǎn)化過程。在此基礎上，提出了一種判斷確定兩多邊形交集區(qū)域范圍的算法，快速計算確定了相應殘骸落區(qū)和規(guī)避區(qū)域在二維坐標系上的交集區(qū)域。最后，通過構(gòu)建二維正態(tài)分布模型定量快速計算導彈殘骸落區(qū)多邊形和規(guī)避區(qū)域多邊形的交集區(qū)域概率，并將交集區(qū)域概率與設置閾值比較。算例結(jié)果顯示，該算法易于實現(xiàn)，計算精度較高。當作戰(zhàn)導彈數(shù)量多的情況下，該方法具有重要的工程應用價值。