許友善，林廣宏，趙 晗，謝亞恩

（1.中國船舶集團有限公司第七二二研究所，湖北 武漢 430200；2.中國船舶及海洋工程設(shè)計研究院，上海200011；3.哈爾濱工程大學(xué)航天與建筑工程學(xué)院，黑龍江 哈爾濱 150001）

引言

我國是自然災(zāi)害多發(fā)的國家之一，其災(zāi)害種類多、分布廣、發(fā)生頻率高、受災(zāi)損失嚴重。在軌衛(wèi)星具有觀測頻次高、遙感范圍廣等特點，在災(zāi)害監(jiān)測評估中發(fā)揮著重要的作用［6］。衛(wèi)星遙感一直是森林草原火災(zāi)、雪災(zāi)、洪澇、地震等自然災(zāi)害的災(zāi)后監(jiān)測和評估的重要手段。

由于受災(zāi)區(qū)域在災(zāi)后的監(jiān)測條件復(fù)雜，傳統(tǒng)的地面監(jiān)測、氣象衛(wèi)星監(jiān)測等手段很難實現(xiàn)對受災(zāi)區(qū)域的全面、高精度的監(jiān)測，使得災(zāi)后救援難以有針對性地進行。因此，可以充分利用遙感衛(wèi)星對受災(zāi)區(qū)域進行監(jiān)測［1-5］。由于受災(zāi)區(qū)域面積較大，衛(wèi)星幅寬有限，因此，用遙感衛(wèi)星對受災(zāi)區(qū)域進行監(jiān)測需要首先對受災(zāi)區(qū)域進行區(qū)域分解，使其能被衛(wèi)星多次推掃覆蓋。傳統(tǒng)的衛(wèi)星在軌觀測任務(wù)規(guī)劃通常為點目標規(guī)劃，隨著敏捷衛(wèi)星的發(fā)展，衛(wèi)星在軌機動能力的提升為大區(qū)域的條帶拼幅成像提供了技術(shù)支持［7］。因此文中針對受災(zāi)區(qū)域災(zāi)后監(jiān)測，給出了基于條帶分解的敏捷衛(wèi)星災(zāi)后快速監(jiān)測方法，能夠有效分析單顆敏捷衛(wèi)星對受災(zāi)區(qū)域的覆蓋觀測時長，為各種自然災(zāi)害的災(zāi)后救援工作提供了數(shù)據(jù)支撐。

1 數(shù)學(xué)模型

1.1 集合和參數(shù)設(shè)計

定義T={v1,v2,...,vN}表示待觀測的受災(zāi)區(qū)域，vi表示目標區(qū)域的第i個頂點，N表示目標區(qū)域的頂點數(shù)量；S={Si|Si=(a,e,i,Ω,ω,tp)}，表示衛(wèi)星的集合，每個衛(wèi)星對應(yīng)一組軌道六要素，a,e,i,Ω,ω,tp分別表示衛(wèi)星的軌道半長軸、軌道偏心率、軌道傾角、升交點赤經(jīng)、近地點幅角、近地點時刻；LLA={LLAi|LLAi=(loni,lati),1 ≤i≤Norbit}表示預(yù)報的衛(wèi)星軌道星下點經(jīng)緯度數(shù)據(jù)集合，Norbit表示預(yù)報的衛(wèi)星軌道點數(shù)；PV={PV1,PV2,...,PVNorbit}表示預(yù)報的衛(wèi)星軌道位置和速度數(shù)據(jù)。Width表示衛(wèi)星的視場幅寬；Cover表示相鄰條帶邊緣所需的重疊寬度。

1.2 變量設(shè)計

m表示進行條帶劃分時，將目標區(qū)域沿垂直于星下點的方向均等分成m個條帶；Strip={strip1,strip2,...,stripm}表示劃分后的條帶集合；stripj={vs,ve}表示第j個條帶劃分后得到的觀測起始點和結(jié)束點；W={W1,W2,...,Wm} 表示計算得到的所有條帶的成像時間窗口；Wj={tjs,tje} 表示衛(wèi)星對條帶stripj的實際成像時間窗口，tjs、tje分別表示實際成像時間窗口的開始時刻、結(jié)束時刻；TM表示2次條帶成像任務(wù)間機動時長；γj,k為二進制變量，表示條帶stripj的執(zhí)行情況，γj,k= 1表示條帶stripj在第k軌被執(zhí)行，否則γj,k= 0［8］。

2 條帶分解

2.1 計算條帶分解的參考直線

計算目標區(qū)域每個頂點與同緯度下每一軌衛(wèi)星星下點之間的經(jīng)度差，第k軌的經(jīng)度差之和計算如下。

則可以得到所有軌道對應(yīng)的經(jīng)度差之和的集合{sum1,sum2,...,sumk,...,summ} ，sumk的最小值對應(yīng)的軌道K即為參考軌道。

確定參考軌道后，可以根據(jù)目標區(qū)域所有頂點與參考軌道的經(jīng)度差{Δlon1,Δlon2,...,ΔlonN} 判斷區(qū)域與參考軌道之間的位置關(guān)系。若所有Δlon均大于零，則區(qū)域在參考軌道右側(cè)且不與參考軌道相交，若所有Δlon均小于零，則區(qū)域在參考軌道左側(cè)且不與參考軌道相交，若Δlon既有大于零的值也有小于零的值，則參考軌道穿過受災(zāi)區(qū)域，若有Δlon等于零，則表明有頂點位于參考軌道之上。

在區(qū)域所有頂點中尋找緯度最大值latmax和最小值latmin，則可以得到中間緯度。

在參考軌道第K軌中尋找latmiddle對應(yīng)的星下點LLAmiddle作為條帶分解的參考點，如圖1中的點LLAmiddle。

圖1 參考直線示意圖Fig.1 Reference line of aera

找到LLAmiddle兩側(cè)的星下點LLA1、LLA2，將其轉(zhuǎn)換到高斯坐標系得到P1=(x1,y1)、P2=(x2,y2（)后文中提到的x,y坐標均為高斯坐標系下坐標）則可以得到參考直線的斜率。

因此條帶分解的參考直線可以寫為：

2.2 確定條帶邊界線

確定條帶邊界線首先需要計算邊界的平移量Δx，再根據(jù)邊界平移量確定最終的條帶數(shù)量m。圖2 給出了條帶邊界平移量示意圖，根據(jù)給定的衛(wèi)星幅寬Width計算得到條帶邊界的平移量如下。

圖2 條帶邊界平移量示意圖Fig.2 Displacement of strip edge

計算所有頂點到參考直線的距離，找到距離最小值dmin和最大值dmax，從而計算得到條帶的數(shù)量如下

確定條帶邊界線首先需要檢查頂點順序并重新排列頂點。若參考直線不穿過受災(zāi)區(qū)域，則判斷v1是否為距離參考直線最近的點，若不是則找到最近的頂點vi作為第一個頂點并將原來的順序進行如下調(diào)換（原頂點是按照順時針順序給出的）：

若參考直線穿過受災(zāi)區(qū)域，則參考圖3，首先找到參考直線兩側(cè)距離最遠的2個頂點v1,v4，在兩點之中找到一個距離相對較近的點v4作為第一個點，同樣按照上述方式對頂點順序進行調(diào)整。

圖3 參考直線穿過受災(zāi)區(qū)域示意圖Fig.3 Reference line crossing disaster aera

在得到邊界平移量和第一個頂點之后，即可得到每個條帶邊界的直線，此處v1指的是調(diào)整順序后得到的第一個頂點。Δx的正負按照左加右減的原則確定，例如受災(zāi)區(qū)域在參考直線左側(cè)，或參考直線穿過受災(zāi)區(qū)域且第一個頂點在直線右側(cè)，條帶邊界線從第一個頂點開始依次向左平移Δx，則Δx為正。同理，若受災(zāi)區(qū)域在參考直線右側(cè)，或參考直線穿過受災(zāi)區(qū)域且第一個頂點在直線左側(cè)，條帶邊界線從第一個頂點開始依次向右平移Δx，則Δx為負。第一個條帶邊界線穿過第一個頂點，平行于參考直線，直線方程如下：

第2個條帶邊界線方程如下：

第j個條帶邊界線方程如下：

由于劃分共得到m個條帶，則條帶邊界共有m+ 1個。

2.3 確定條帶的觀測起始點和結(jié)束點

條帶的觀測起始點和結(jié)束點確定的垂直于參考直線的2 條邊界線我們稱其為起始邊線和結(jié)束邊線，條帶的起始邊線和結(jié)束邊線可能經(jīng)過條帶中存在的受災(zāi)區(qū)域頂點，也可能經(jīng)過條帶邊界線與受災(zāi)區(qū)域邊界線的某一個交點，因此要確定條帶的起始邊線和結(jié)束邊線，首先需要找到條帶中兩端邊界線所有可能的經(jīng)過的點。

首先找到條帶中存在的受災(zāi)區(qū)域的頂點。依次判斷每個條帶與每個頂點之間的位置關(guān)系，得到每個條帶中的頂點，若頂點位于某一條帶邊界線上，則兩側(cè)條帶所考慮的受災(zāi)區(qū)域頂點中均包括該頂點，最終可以得到每個條帶中的頂點集合Vertex={vertex1,vertex2,...,vertexm} 。其中，vertexj={v1,...,vi} ，vertexj也可能為空集，即條帶中不存在頂點。

在確定條帶中存在的受災(zāi)區(qū)域的頂點后，再找到所有區(qū)域邊界線與條帶邊界線之間的交點。根據(jù)集合Vertex={vertex1,vertex2,...,vertexm} 可以確定每個頂點所在的條帶，假設(shè)2個相鄰的頂點vi和vi+1分別位于條帶s1和s2中，則頂點vi和vi+1之間的條帶邊界線為

式中(s1,s2]表示條帶邊界線s1+ 1,s1+ 2,...,s2，2 個相鄰頂點vi和vi+1之間共s2-s1個條帶邊界線。因此與受災(zāi)區(qū)域邊界edgei有交點的條帶邊界線集合={strips1+1,...,strips2} 。

計算受災(zāi)區(qū)域邊界edgei與條帶邊界線={strips1+1,...,strips2} 的交點，記錄到每個條帶的頂點集合Vertex={vertex1,vertex2,...,vertexm} 中。則頂點集合中包括每個條帶中存在的受災(zāi)區(qū)域頂點和條帶邊界線與受災(zāi)區(qū)域邊界的交點。

作一條過頂點v1并垂直于參考直線l0的直線L，直線方程如下：

計算每個條帶中的所有頂點qp到直線L的距離，如圖4所示：

式中，d的值可能小于零。

從大到?。紤]d的正負）對計算得到的某一條帶中的所有d進行排序，找到d的最大值和最小值對應(yīng)的2 個點，即為條帶起始和結(jié)束兩端直線需要經(jīng)過的點。以圖4 中的條帶strip2為例，計算條帶中的4個頂點與直線L的距離，由于頂點位于L兩側(cè)，因此d有正負，此處正負只表示方位。對4個頂點的d進行排序，則可以確定條帶最上方和最下方的2個頂點為v5、q3。分別過v5、q3作垂直于l0的直線，并計算直線與條帶邊界線的交點，得到條帶的4個頂點{point1,point2,point3,point4} ，取條帶的起始邊線和結(jié)束邊線的中點為條帶兩端的觀測起始點和結(jié)束點，則條帶的觀測起始點和結(jié)束點可以表示為

圖4 條帶兩端頂點計算示意圖Fig.4 Vertex calculation figure of strip

最后判斷星下點的走向，以圖1 中的LLAmiddle為t0時刻，則LLA1={lon1,lat1} 、LLA2={lon2,lat2} 分別對應(yīng)t0- 1和t0+ 1時刻。

若lat1＜lat2，則衛(wèi)星上行，則選取point01、point02中緯度較小的點為觀測起始點，另一點為觀測結(jié)束點，反之衛(wèi)星下行，選取point01、point02中緯度較大的點為起始點，另一點為觀測結(jié)束點。由于上述計算是在高斯坐標系中進行的，因此需要將得到的觀測起始點和結(jié)束點轉(zhuǎn)換為經(jīng)緯度，得到條帶stripj的觀測起始點和結(jié)束點stripj={vs,ve}。最終計算m個條帶的觀測起始點和結(jié)束點，得到Strip={strip1,strip2,...,stripm}用作成像任務(wù)規(guī)劃的輸入。

3 成像規(guī)劃

計算條帶stripj的觀測起始點的開始可見時刻，并作為條帶stripj的實際觀測開始時刻tjs。根據(jù)衛(wèi)星的地速vg計算條帶stripj的觀測時長，具體計算如下：

式中P1,P2分別表示條帶stripj的起始觀測點和結(jié)束觀測點在地固坐標系中的坐標。因此，當(dāng)條帶的實際觀測開始時刻確定后，其實際觀測的結(jié)束時刻為：

衛(wèi)星的可見性如圖5所示，本文假設(shè)衛(wèi)星是敏捷衛(wèi)星，軌道偏心率為0，衛(wèi)星最近可對星下點進行觀測，最遠可觀測到與衛(wèi)星地球表面的切點。衛(wèi)星可見性約束可以表述為0 ＜θ＜θmax，其中RE、RS分別表示地球半徑和衛(wèi)星半長軸。

圖5 衛(wèi)星可見性示意圖Fig.5 Satellite visibility figure

根據(jù)上述可見性約束，依次計算條帶Strip={strip1,strip2,...,stripm}的可見性。

遍歷衛(wèi)星的LLA數(shù)據(jù)，計算衛(wèi)星與目標間的地心角θ，若0 ＜θ＜θmax，則當(dāng)前時刻衛(wèi)星對目標可見，并以當(dāng)前時刻作為條帶stripj的實際觀測開始時刻tjs，并可以通過式（17）計算得到條帶stripj成像結(jié)束時刻。重復(fù)上述方法依次計算每個條帶的實際成像時間。但是在計算每個條帶的實際成像時間時，還需要考慮任務(wù)間的沖突，即條帶stripl的成像時段對于已規(guī)劃的任務(wù)序列中的任意其他條帶stripj應(yīng)滿足tls-TM≥tje或tle+TM≤tjs。

最終，遍歷所有條帶，計算得到整個受災(zāi)區(qū)域的所有條帶的實際成像時間窗口W={W1,W2,...,Wm} 。

4 仿真驗證

4.1 仿真場景設(shè)計

表1 給出了仿真用到的衛(wèi)星軌道參數(shù)。衛(wèi)星的軌道預(yù)報時長為24 h，步長為1s，圖6 中的曲線即為衛(wèi)星的所有星下點軌跡。衛(wèi)星的視場幅寬Width= 60 km，相鄰條帶邊緣所需的重疊寬度Cover= 2 km。

表1 衛(wèi)星軌道參數(shù)Table 1 Orbital parameters of satellite

（1）仿真場景1

考慮到四川地區(qū)自然災(zāi)害頻發(fā)，因此在四川附近地區(qū)選取一個五邊形區(qū)域作為受災(zāi)區(qū)域1，表2給出了目標區(qū)域1的頂點，圖6給出了受災(zāi)區(qū)域1與軌道的示意圖。

表2 目標區(qū)域1頂點Table 2 Vertex of disaster aera 1

圖6 軌道與受災(zāi)區(qū)域1示意圖Fig.6 Orbit and disaster aera 1

圖7 軌道與受災(zāi)區(qū)域2示意圖Fig.7 Orbit and disaster aera 2

（2）仿真場景2

為驗證算法在處理受災(zāi)區(qū)域與衛(wèi)星下點不同位置關(guān)系下的準確性，設(shè)計一種受災(zāi)區(qū)域位于星下點兩側(cè)的仿真場景，表3給出了目標區(qū)域的頂點，圖6給出了受災(zāi)區(qū)域與軌道的示意圖。

表3 目標區(qū)域2頂點Table 3 Vertex of disaster aera 2

4.2 仿真

（1）仿真場景1計算結(jié)果

圖8 給出了條帶計算結(jié)果，圖9 給出了條帶實際成像時段規(guī)劃結(jié)果示意圖。在不考慮軌道機動的情況下，以本文算例給出的軌道初始條件，得到的所有條帶的最早成像時刻為19842 s（5.511 h），成像結(jié)束時刻為21 003 s（5.834 h），從開始成像到成像結(jié)束共耗時1 161 s（0.323 h），其中成像時長共1 004 s（0.279 h），任務(wù)間的機動及等待成像共耗時157 s。

圖8 受災(zāi)區(qū)域1條帶分解結(jié)果Fig.8 Strips splitting result of disaster aera 1

圖9 受災(zāi)區(qū)域1成像時段Fig.9 Imaging time result of disaster aera 1

（2）仿真場景2結(jié)果

圖10 給出了條帶計算結(jié)果，圖11 給出了條帶實際成像時段規(guī)劃結(jié)果示意圖。在不考慮軌道機動的情況下，以本文算例給出的軌道初始條件，得到的所有條帶的最早成像時刻為19 991 s（5.553 h），成像結(jié)束時刻為20 951 s（5.820 h），從開始成像到成像結(jié)束共耗時960 s（0.267 h），其中成像時長共820 s（0.228 h），任務(wù)間的機動及等待成像共耗時140 s。

圖10 受災(zāi)區(qū)域2條帶分解結(jié)果Fig.10 Strips splitting result of disaster aera 2

圖11 受災(zāi)區(qū)域2成像時段Fig.11 Imaging time result of disaster aera 2

5 結(jié)論

針對災(zāi)后小范圍區(qū)域觀測，氣象衛(wèi)星則難以滿足精度要求的問題，文中給出了基于條帶分解的受災(zāi)區(qū)域衛(wèi)星遙感觀測方法，提供了災(zāi)后小范圍區(qū)域的觀測思路。通過將受災(zāi)區(qū)域進行條帶分解，使敏捷遙感衛(wèi)星能夠?qū)κ転?zāi)區(qū)域進行單軌連續(xù)觀測或多軌連續(xù)觀測，從而達到對受災(zāi)區(qū)域的災(zāi)后應(yīng)急觀測的目的。本文通過仿真算例，驗證了本文給出的基于條帶分解的區(qū)域目標的成像時段的規(guī)劃方法，算法可行有效，為災(zāi)后小范圍區(qū)域的災(zāi)后監(jiān)視提供了參考。