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

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

引言

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

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

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

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

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

1.2 變量設(shè)計(jì)

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

2 條帶分解

2.1 計(jì)算條帶分解的參考直線

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

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

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

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

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

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

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

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

2.2 確定條帶邊界線

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

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

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

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

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

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

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

第2個(gè)條帶邊界線方程如下：

第j個(gè)條帶邊界線方程如下：

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

2.3 確定條帶的觀測(cè)起始點(diǎn)和結(jié)束點(diǎn)

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

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

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

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

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

作一條過(guò)頂點(diǎn)v1并垂直于參考直線l0的直線L，直線方程如下：

計(jì)算每個(gè)條帶中的所有頂點(diǎn)qp到直線L的距離，如圖4所示：

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

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

圖4 條帶兩端頂點(diǎn)計(jì)算示意圖Fig.4 Vertex calculation figure of strip

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

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

3 成像規(guī)劃

計(jì)算條帶stripj的觀測(cè)起始點(diǎn)的開(kāi)始可見(jiàn)時(shí)刻，并作為條帶stripj的實(shí)際觀測(cè)開(kāi)始時(shí)刻tjs。根據(jù)衛(wèi)星的地速vg計(jì)算條帶stripj的觀測(cè)時(shí)長(zhǎng)，具體計(jì)算如下：

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

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

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

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

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

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

4 仿真驗(yàn)證

4.1 仿真場(chǎng)景設(shè)計(jì)

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

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

（1）仿真場(chǎng)景1

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

表2 目標(biāo)區(qū)域1頂點(diǎn)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）仿真場(chǎng)景2

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

表3 目標(biāo)區(qū)域2頂點(diǎn)Table 3 Vertex of disaster aera 2

4.2 仿真

（1）仿真場(chǎng)景1計(jì)算結(jié)果

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

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

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

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

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

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

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

5 結(jié)論

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