張 瑋,王守斌,程承旗,陳 波,李海東
(1. 北京大學(xué) 地球與空間科學(xué)學(xué)院,北京 100871;2. 61646部隊 北京 100192;3. 北京大學(xué) 工學(xué)院,北京 100871)
時間窗口是指衛(wèi)星經(jīng)過地面目標的上空時,可對該目標進行觀測的時間范圍。時間窗口的計算具有重要意義,在觀測衛(wèi)星調(diào)度中有著重要作用[1-2]。
相對于點目標而言,區(qū)域目標的時間窗口計算更加困難:區(qū)域目標形狀可能不規(guī)則,時間窗口計算通常只能利用空間采樣點或區(qū)域邊界線進行近似求解,而這種求解的瓶頸在于時間效率。
當前,衛(wèi)星對區(qū)域目標的時間窗口計算主要有如下兩類方法:①基于空間采樣點的區(qū)域目標時間窗口計算方法。較為典型的有:跟蹤傳播法[3]、星下點大圓近似迭代法[4]、視函數(shù)法[5]、龐加萊映射解析法[6]等。當采樣點劃分尺度和時間步長足夠細時,跟蹤傳播法的時間窗口計算結(jié)果非常精確,但缺點是時耗很高。星下點大圓近似迭代法的算法精確度較低,同時存在無法收斂的可能。視函數(shù)法的精度很高,但區(qū)域目標空間范圍較大時,采樣點數(shù)量劇增,時間窗口的實時計算時效較差。龐加萊映射法基于數(shù)值計算,能精確獲得時間窗口,但計算過程涉及大量積分運算,算法復(fù)雜度高,對區(qū)域目標進行計算時效率較低。②基于邊界線段的區(qū)域目標時間窗口求解方法。較為典型的方法是星下點軌跡-邊界求交法[1,2,7],其優(yōu)勢是計算衛(wèi)星對大范圍的區(qū)域目標的時間窗口時相對基于空間采樣點的計算方法的算法效率更高,但缺點在于計算時耗與區(qū)域目標邊界的數(shù)量成正比,如果多邊形區(qū)域本身范圍不大但邊界數(shù)量很短碎時,反而時耗會比基于空間采樣點的方法更多;同時,該方法不適用于地球橢球體的情況。
上述兩種方法均無法存儲衛(wèi)星地面覆蓋信息,因此衛(wèi)星訪問目標計算的過程是相互獨立的,如圖1所示,其中存在大量實時、重復(fù)的覆蓋計算,導(dǎo)致耗時嚴重。
圖1 現(xiàn)有時間窗口計算方法:實時覆蓋和獨立過程計算處理Fig.1 The existing methods in time window computing access time:the real-time and independent coverage computing process
本文擬提出一種適用于衛(wèi)星訪問區(qū)域目標的時間窗口計算方法。這種方法將衛(wèi)星覆蓋信息提前存儲下來,衛(wèi)星訪問區(qū)域目標的時間窗口實時進行復(fù)雜計算轉(zhuǎn)化為覆蓋信息的查詢、整合,因此既能保持較高的計算精度,還具備更高的時間效率。同時,多星并行計算過程轉(zhuǎn)變?yōu)獒槍^(qū)域目標的一次覆蓋查詢過程,計算耗時相對穩(wěn)定,因此更具優(yōu)勢。
傳統(tǒng)衛(wèi)星時間窗口計算通常采用等經(jīng)緯度、等面積或等距離等方式去生成網(wǎng)格,再用網(wǎng)格點集合來近似擬合區(qū)域目標。本文同樣也需要借助網(wǎng)格,但立足點不同:①傳統(tǒng)方法使用的是分布均勻的離散點,網(wǎng)格是生成這些點的方法,而本文方法需要的是網(wǎng)格本身,因為網(wǎng)格代表了固定的小尺度地表區(qū)域;②傳統(tǒng)方法用網(wǎng)格點來表征目標,本文方法用網(wǎng)格來存儲衛(wèi)星覆蓋的時空信息。本文選擇北京大學(xué)程承旗教授團隊提出的全球等經(jīng)緯度剖分網(wǎng)格體系(Geographical coordinate global Subdivision grid with One - dimension - integer on Two to n-th power,GeoSOT)[8]作為衛(wèi)星時空覆蓋網(wǎng)格的空間基準,剖分編碼原理如圖2所示。GeoSOT的優(yōu)勢在于:①全球尺度無縫無疊;②網(wǎng)格粒度適宜,且多層級嵌套;③網(wǎng)格編碼唯一且便于計算機存儲、計算。
圖2 衛(wèi)星地面覆蓋區(qū)域剖分成同一層級網(wǎng)格Fig.2 The satellite coverage area subdivided into the same layer of GeoSOT grids
本文采用時間剖分編碼[9]作為衛(wèi)星覆蓋時間信息存儲載體。時間剖分編碼借鑒了四叉樹倒排聚合的做法,將時間剖分形成同時支持單尺度和多尺度表達的64位二進制1維編碼形式,既方便數(shù)據(jù)庫存儲和索引,又能基于編碼進行時間聚合和跨度運算。同時,這種編碼還具有十進制表達數(shù)值越小、時間越早,數(shù)值越大、時間越晚的特點,因此適合于時刻排序過程。
該方法需要計算瞬時時刻衛(wèi)星地面覆蓋區(qū)域。目前,衛(wèi)星地面覆蓋計算方法研究較多,主要有球面三角形法和矢量法兩大類。本文主要采用矢量法[10]進行衛(wèi)星地面覆蓋計算。矢量法的優(yōu)勢在于計算精度較高,同時兼容地球球體和地球橢球體的情形,且還兼容不同視場形狀的衛(wèi)星傳感器。
主要思路是:①提前計算衛(wèi)星地面覆蓋,剖分化后用網(wǎng)格記錄衛(wèi)星的時空覆蓋信息,存儲至數(shù)據(jù)庫中;②區(qū)域目標訪問分析計算轉(zhuǎn)化為按邊界網(wǎng)格遍歷查詢數(shù)據(jù)庫的過程,尋找衛(wèi)星覆蓋區(qū)域目標邊界的所有時刻,生成時間窗口;③多星分析區(qū)域目標時,只需一次遍歷查詢即可獲得所有衛(wèi)星訪問區(qū)域目標的所有時間窗口。技術(shù)路線流程圖如圖3所示。
圖3 技術(shù)路線流程圖Fig.3 Technical route fl ow chart
衛(wèi)星地面覆蓋預(yù)存儲是指將仿真時域按等間隔劃分成一系列時間區(qū)間,依次計算時間區(qū)間的邊界時刻下衛(wèi)星地面覆蓋區(qū)域,并將覆蓋區(qū)域網(wǎng)格化,將該時刻的所有覆蓋網(wǎng)格的網(wǎng)格編碼附上該時刻生成的時間剖分編碼,存入數(shù)據(jù)庫,最終形成衛(wèi)星覆蓋數(shù)據(jù)庫的過程。
2.2.1 仿真時域的劃分
假設(shè)仿真時域的起始時刻為Tbegin,終止時刻為Tend,設(shè)定時間步長為Δt,那么該時域可以分割成數(shù)量為N的小時域。從第0段開始到第N-1段,時間段的起始和終止時刻分別為:[Tbegin,Tbegin+Δt],[Tbegin+Δt,Tbegin+2Δt],…,[Tend-Δt,Tend]。其中,前一個小時域的終止時刻是后一個小時域的起始時刻,那么將起始時刻為Tbegin、終止時刻為Tend以及所有分割時刻集合,可得
此時,T被稱為仿真時域的時刻集。
2.2.2 衛(wèi)星瞬時地面覆蓋區(qū)域計算與剖分化
遍歷T中的時刻:
Step1:求解該時刻下衛(wèi)星的地面覆蓋區(qū)域,結(jié)果用區(qū)域邊界點的經(jīng)緯度序列CoverLine _pt_seq表示。
Step2:選擇符合衛(wèi)星應(yīng)用精度需求的最高網(wǎng)格層級(即基礎(chǔ)層級,一般可選擇第14層級,該層級網(wǎng)格大小在赤道附近約為4 km左右),將計算得到的衛(wèi)星傳感器地面覆蓋區(qū)域剖分至該網(wǎng)格層級,如圖2所示,得到一個由網(wǎng)格組成的集合Cs,t->e;用該集合來表征衛(wèi)星傳感器此時刻的地面覆蓋,即:
Step3:為了盡可能減少數(shù)據(jù)量,將Cs,t->e中的網(wǎng)格聚合,形成多尺度網(wǎng)格集合Cs,t->e,multi,如圖4所示。
圖4 地面覆蓋網(wǎng)格聚合Fig.4 Satellite coverage grids merged into multi-layer coverage grids
2.2.3 衛(wèi)星覆蓋時空網(wǎng)格編碼入庫
遍歷T中的時刻Ti,都能得到與之對應(yīng)的Cs,t->e,multi。將Cs,t->e,multi編碼化可以得到多層級的編碼集合Gs,t→e:
同時,將Ti編碼化為Tcode,那么結(jié)合Gs,t→e中的每個空間網(wǎng)格編碼元素,都能形成一個
衛(wèi)星覆蓋數(shù)據(jù)庫的表結(jié)構(gòu)字段應(yīng)當包括:
1)衛(wèi)星編號(satID):INT32類型;
2)傳感器編號(senID):INT32類型;
3)衛(wèi)星傳感器名稱(Name):String類型;
4)空間網(wǎng)格編碼(gcode):UNSIGNED LONG LONG或 __INT64類型;
5)時間剖分編碼(tcode):UNSIGNED LONG LONG或 __INT64類型;
6)其他屬性碼。
其中,數(shù)據(jù)庫中每個數(shù)據(jù)條目代表的含義是:某顆衛(wèi)星的某個傳感器在某時刻(或分割時域,取決于時間精度)能夠覆蓋某個網(wǎng)格。此時,衛(wèi)星傳感器地面覆蓋網(wǎng)格的映射關(guān)系全部打散、存儲至數(shù)據(jù)庫的一張表中,這張表的條目數(shù)量就是所有衛(wèi)星在所有時刻的所有網(wǎng)格編碼的總和。在此基礎(chǔ)上,基于空間網(wǎng)格編碼建立數(shù)據(jù)庫正序普通索引。
區(qū)域目標的時間窗口分析的基本做法是,先將區(qū)域目標的邊界點序列轉(zhuǎn)換成基礎(chǔ)層級的網(wǎng)格,再將網(wǎng)格代入數(shù)據(jù)庫中查詢,返回與網(wǎng)格相關(guān)的所有時間編碼,最后將時間編碼整合生成時間窗口。
2.3.1 區(qū)域目標邊界剖分化
以點序列輸入的區(qū)域目標邊界后,將點序列內(nèi)相鄰點兩兩連接,形成若干條邊界線段,依次遍歷所有線段,將線段轉(zhuǎn)換成基礎(chǔ)層級下的線段網(wǎng)格集合,如圖5所示;最后將所有集合合并,得到區(qū)域目標的邊界網(wǎng)格編碼集合。
圖5 區(qū)域目標邊界剖分化Fig.5 The area target boundary differentiation
2.3.2 數(shù)據(jù)庫查詢
將邊界網(wǎng)格編碼集合中的網(wǎng)格編碼依次代入衛(wèi)星覆蓋數(shù)據(jù)庫中查詢,獲得覆蓋所有邊界網(wǎng)格時的時間編碼。由于數(shù)據(jù)庫存儲的是多層級網(wǎng)格時空覆蓋信息,因此如果該網(wǎng)格能被衛(wèi)星覆蓋,可能存在兩種情況:
1)該網(wǎng)格編碼存在于數(shù)據(jù)庫中;
2)該網(wǎng)格編碼不存在于數(shù)據(jù)庫,但該網(wǎng)格的多層級父網(wǎng)格存在于數(shù)據(jù)庫中。
因此,應(yīng)當使用代入網(wǎng)格編碼,如下兩個條件滿足其一即返回的查詢條件進行查詢:
條件1. 等值判斷,空間編碼與查詢輸入相等,此時該輸入網(wǎng)格存在于數(shù)據(jù)庫中;
條件2. 異或求交判斷,查詢輸入與空間編碼的經(jīng)向編碼與緯向編碼進行一次異或運算,并對兩次異或運算結(jié)果進行求交運算,結(jié)果值為1。
按2.4.1節(jié)得到的邊界網(wǎng)格集合是順序排列的??紤]到衛(wèi)星在仿真時間內(nèi)進出區(qū)域目標的次數(shù)(假設(shè)為N次)是有限的,衛(wèi)星對區(qū)域邊界的覆蓋集中在2N個邊界部分,在這些部分之外的區(qū)域邊界不被覆蓋。因此,在衛(wèi)星覆蓋數(shù)據(jù)庫查詢中,可以采用跳躍查詢的方式縮減查詢次數(shù)。具體做法是:
Step 1:設(shè)置合適的跳躍步長step,step=2n,設(shè)當前查詢網(wǎng)格的集合索引值為record,設(shè)回退狀態(tài)布爾值為isStep為false;
Step 2:從第0個網(wǎng)格開始遍歷,即record初始化為0:
Step ①:將record指向的網(wǎng)格編碼代入數(shù)據(jù)庫中查詢;
Step ②:如果返回為空,則判斷record+step是否大于邊界網(wǎng)格集合的元素個數(shù),如果不是,則record+=step,進入Step①,如果是,則record指向最后一個元素,進入Step①;如果返回不為空,則判斷isStep是否為true,如果為true,則進入回退搜索臨近最遠的非空返回的索引值,如圖6所示,如果為false,則記錄該網(wǎng)格的查詢返回時間編碼,如果此時record指向最后一個元素,則進入Step3,否則record+1,進入Step①;
Step 3:遍歷結(jié)束,將時間編碼集合并去重。
圖6 回退搜索過程示意圖Fig.6 The processing of the rollback search method
2.3.3 生成時間窗口
獲得所有覆蓋邊界網(wǎng)格的時間編碼后,按值大小排序。此時,相鄰元素之間的時間差值存在3種情況:①差值等于時間步長(如1 s);②差值略大于時間步長(如10 s);③差值遠大于時間步長(如30 min)。情況①說明相鄰元素均屬于同一個時間窗口,較早者可能為時間窗口的起始時刻;情況②說明相鄰元素屬于同一個時間窗口,較早者為衛(wèi)星進入?yún)^(qū)域時覆蓋邊界的最晚時刻,較晚者為衛(wèi)星離開區(qū)域時覆蓋邊界的最早時刻;情況③說明相鄰元素不屬于同一個時間窗口,較早者為前一個時間窗口的終止時刻,較晚者為后一個時間窗口的起始時刻。
基于上述情況,生成時間窗口的主要做法是利用區(qū)域邊界的經(jīng)緯度極值,生成區(qū)域目標的最小外包矩形,計算外包矩形的較長對角大圓弧長度,再除以衛(wèi)星星下點軌跡運動平均速度,以此預(yù)估衛(wèi)星一次覆蓋區(qū)域目標可能出現(xiàn)的時間窗口最大值△T;以該預(yù)估值為臨界值,將時間編碼集合分簇,然后獲得各簇的最小和最大編碼,作為時間窗口的起始時刻和終止時刻,如圖7所示。
圖7 時刻分簇生成時間窗口示意圖Fig.7 Time clustering generation time windows
基于單星的地面覆蓋預(yù)存儲結(jié)果,分析計算兩個假想?yún)^(qū)域目標的時間窗口,并以跟蹤傳播法作為比較對象,對比分析本文方法的結(jié)果精度和時間效率;同時,對多星地面覆蓋進行存儲,計算一個區(qū)域目標的時間窗口,對比分析本文方法和傳統(tǒng)方法多星計算的時間效率。
實驗的硬件環(huán)境和軟件環(huán)境分別見表1和表2。
表1 硬件環(huán)境Tab.1 Hardware Environment of the simulation experiment
表2 軟件環(huán)境Tab.2 Software Environment of the simulation experiment
使用20顆衛(wèi)星進行實驗,其中單星分析所用的衛(wèi)星軌道根數(shù)見表3。仿真區(qū)域目標為兩個,目標1為中國大陸區(qū)域(不含中國臺灣、海南等島嶼、海域),邊界點數(shù)量為4 626,目標2為中國南海部分區(qū)域,邊界點經(jīng)緯度坐標見表4,不考慮衛(wèi)星的姿態(tài)角變化,假定視場中心點投影到地表處與星下點始終重合,同時設(shè)定水平半視場角和垂直半視場角均為15°。仿真時域的起始時刻為2018年5月21日9點整(UT0),終止時刻為2018年5月22日9點整(UT0),時域跨度為1 d(即86 400 s)。
表3 衛(wèi)星軌道根數(shù)Tab.3 Orbital elements of the satellite in the single-satellite experiment
表4 仿真區(qū)域目標邊界點經(jīng)緯度Tab.4 Geographical coordinates of the boundary points of the area targets
3.3.1 單星計算實驗結(jié)果與分析
在衛(wèi)星覆蓋預(yù)存儲階段,選定時間步長為1 s,覆蓋網(wǎng)格基礎(chǔ)層級為GeoSOT第14層級(2′網(wǎng)格,網(wǎng)格大小對應(yīng)赤道附近4 km左右),計算衛(wèi)星覆蓋時空網(wǎng)格并完成入庫。輸入目標1和目標2的邊界點序列,分析兩者的時間窗口,計算結(jié)果見表5。
表5 目標1和目標2的時間窗口計算結(jié)果Tab.5 Time window calculation results of target 1 and target 2
將目標2的起始時刻和終止時刻轉(zhuǎn)換為時間編碼,代入數(shù)據(jù)庫中查詢相關(guān)網(wǎng)格編碼,再將網(wǎng)格繪制成圖,與STK二維標準輸出窗口結(jié)果進行比對可知,本文采用的數(shù)據(jù)庫查詢方法是可靠的,如圖8所示。
圖8 本文方法結(jié)果與傳統(tǒng)方法結(jié)果的對比Fig.8 Comparison of the results of this method with the traditional method
以STK分析結(jié)果為標準值,采用絕對誤差Δ和相對誤差δ[1]作為分析結(jié)果精度的指標。
絕對誤差的計算公式為:
式(4)、(5)中,N為時間窗口的個數(shù),tSi,true和tEi,true分別是第i個時間窗口真實值的起始時刻和終止時刻,tSi和tEi分別是通過本文方法得到的第i個時間窗口的起始時刻和終止時刻。
絕對誤差反映了計算所得的時間窗口與實際時間窗口在數(shù)值上的絕對差別。相對誤差則反映了誤差相比于時間窗口的大小。
取跟蹤傳播法的時間步長為1 s,計算時間窗口,將所得結(jié)果與本文方法進行對比,計算時耗與結(jié)果精度見表6和表7。
表6 本文方法與跟蹤傳播法的誤差對比Tab.6 Time error comparison between this method and the tracking propagation method
表7 本文方法與跟蹤傳播法的計算效率對比Tab.7 Comparison of computational efficiency between this method and the tracking propagation method
由表6可知,與跟蹤傳播法相比,本文方法的絕對誤差為秒級誤差,與其他快速計算方法持平[11]。誤差主要取決于如下3個方面:
1)時間步長。時間步長越小,仿真時域的劃分就越細,相鄰時刻計算得到的衛(wèi)星地面覆蓋區(qū)域重合度越高,單個覆蓋網(wǎng)格包含的時間覆蓋信息就越豐富(原理類似于視頻幀數(shù));反之,則單個覆蓋網(wǎng)格包含的時間覆蓋信息就越少。當時間步長取值為1 s時,衛(wèi)星在1 s時間內(nèi)的覆蓋動態(tài)變化信息就無法獲取,那么網(wǎng)格只能記錄1 s時間能否被衛(wèi)星覆蓋,至于這1 s內(nèi)是否都能被覆蓋,無從得知。因此,時間步長直接決定了時間窗口分析精度的分辨率,即基于1 s步長的預(yù)存儲覆蓋數(shù)據(jù)進行區(qū)域目標查詢分析的結(jié)果絕對誤差不可能小于1 s。
2)網(wǎng)格基礎(chǔ)層級。本文提出的方法本質(zhì)上是將衛(wèi)星覆蓋能力提前計算、存儲至剖分網(wǎng)格中,剖分網(wǎng)格是衛(wèi)星時空覆蓋信息的載體。網(wǎng)格劃分越細,對衛(wèi)星瞬時覆蓋區(qū)域的擬合就越好,基于區(qū)域邊界網(wǎng)格的時間窗口查詢分析結(jié)果也越收斂、逼近真實時間窗口。
3)衛(wèi)星瞬時覆蓋的計算精度。衛(wèi)星瞬時覆蓋的結(jié)果精度也會影響本文方法的結(jié)果精度。瞬時覆蓋精度越高,數(shù)據(jù)庫存儲的時空覆蓋信息就越接近真實情況,最終經(jīng)過查詢檢索得到的時間窗口也越精確。
表7可知,本文提出的方法計算效率相對于傳統(tǒng)的跟蹤傳播法更高,主要原因是該方法將復(fù)雜覆蓋計算變成了預(yù)存儲過程,面向具體區(qū)域目標時,只需要判斷網(wǎng)格與網(wǎng)格之間的包含關(guān)系,而GeoSOT采用了64位二進制Z序網(wǎng)格編碼的形式,網(wǎng)格包含關(guān)系都可用編碼位運算解決,因此區(qū)域目標的時間窗口計算過程被分離、簡化,得以提升時間效率。
3.3.2 多星計算實驗結(jié)果與分析
將20顆衛(wèi)星地面覆蓋存儲,對目標2進行多星計算。以跟蹤傳播法為對比,計算效率結(jié)果如圖3-2所示,橫坐標為衛(wèi)星數(shù)量,縱坐標為傳統(tǒng)方法計算耗時與本文方法耗時的比率。
圖9 多星分析效率比較Fig.9 Computational efficiency comparison of the multi-satellite analysis
由圖9可知,隨著衛(wèi)星數(shù)量的增加,本文方法的時間效率比傳統(tǒng)方法越來越高,其原因在于多星地面覆蓋都用網(wǎng)格進行存儲,而基于邊界網(wǎng)格進行覆蓋時刻查詢時,可以做到一次遍歷、同時獲取所有衛(wèi)星的時刻,進而得到所有衛(wèi)星的時間窗口,因此用本文方法進行多星計算時,時間損耗是較為穩(wěn)定的;而傳統(tǒng)方法進行多星計算時,需要分別對每顆衛(wèi)星計算訪問目標的時間窗口,因此時間損耗隨著衛(wèi)星數(shù)量的增加而增加。因此,本文方法更適合于多星計算的情形。
本文提出了一種基于網(wǎng)格的衛(wèi)星訪問區(qū)域目標的時間窗口計算方法,即將衛(wèi)星地面覆蓋提前計算、剖分化形成時空覆蓋網(wǎng)格信息,存儲至數(shù)據(jù)庫中,在計算區(qū)域目標的時間窗口時,將區(qū)域邊界剖分得到的網(wǎng)格遍歷代入數(shù)據(jù)庫查詢,尋找衛(wèi)星覆蓋區(qū)域目標邊界的所有時刻,生成時間窗口。該方法具有衛(wèi)星覆蓋信息“一次存儲、反復(fù)查詢”的技術(shù)優(yōu)勢,仿真試驗證明,該方法比傳統(tǒng)方法更適用于多星時間窗口計算,因此可為多星聯(lián)合快速調(diào)度分析提供技術(shù)支撐。該方法在衛(wèi)星覆蓋計算及網(wǎng)格存儲方面,未考慮敏捷衛(wèi)星姿態(tài)角變化、傳感器側(cè)擺等因素,同時需要消耗一定的計算機存儲空間。因此,下一步工作需要改進時空覆蓋網(wǎng)格存儲的方法,以便兼容衛(wèi)星姿態(tài)角、傳感器側(cè)擺等復(fù)雜情況下的衛(wèi)星覆蓋信息存儲,同時還要進一步壓縮數(shù)據(jù)庫大小,以便節(jié)省更多的空間開支。