呂林立，肖歆昕，馮冠華，李文皓,*

1. 中國科學院 力學研究所，北京 100190 2. 中國科學院大學 工程科學學院，北京 100049

隨著商業(yè)航天的蓬勃發(fā)展，巨型星座網絡近年來在建設、組網、服務等多方面取得快速進展。例如，美國SpaceX公司的星鏈(Starlink)計劃部署40 000多顆衛(wèi)星，目前已經部署了超過1 300顆，超過了其他低軌衛(wèi)星在軌數量的總和，其他如美國的Kuiper計劃、英國的OneWeb、中國的“鴻雁”“鴻云”以及各類部署組網的Telesat和Galaxy Space計劃等。

在巨型星座之前，對地和對天覆蓋性已經得到了廣泛研究，并用于對地觀測分析(如NASA的EOS計劃)、空間機器人遙操作的中繼分析等。對于巨型星座網絡，其主要預期服務包括：對地觀測、天基互聯網通訊服務等。這些服務的實現要求對地有效覆蓋、星間有效互聯，因此巨型星座的對地、對天覆蓋性準確分析是其產生預期能力的前提。此外，巨型星座的建立需要多次部署，逐漸組網，包括建設和服務的過程中，也存在因部分衛(wèi)星失效而產生補網需求，因此，動態(tài)和高效地分析巨型星座的覆蓋性非常重要。其計算效率直接影響著星座的設計效率。用于覆蓋性分析的覆蓋帶法，僅能分析圓軌道星座，基于三角剖分的覆蓋性分析方法，該類方法僅能分析較簡單覆蓋區(qū)域，也有基于映射的二維圖法。對于復雜星座構型，由于星座的覆蓋性能難以用解析的方法求解，網格法為廣泛使用的方法。網格法最早由Morrison提出，該方法的基本思路是將目標區(qū)域的包圍盒按經緯度、距離、面積等方式劃分為網格，以網格中心點代表網格，按給定步長計算衛(wèi)星空間位置極點，由地心角法計算并記錄衛(wèi)星對網格點的覆蓋情況。

由于網格法的分析是對星座中的各個衛(wèi)星，均逐個網格遍歷求解，因此隨著星座中衛(wèi)星數量的增加，其覆蓋性計算的消耗必然會上升，對于10，甚至10量級的星座，使用網格法的覆蓋性分析對計算能力和計算效率的需求會直線上升。

本文提出一種覆蓋性分析方法：并行墨卡托投影圖疊加法(PMPS)，該方法將星座在指定空間天球的覆蓋域，通過圖像畸變投影至二維墨卡托圖上，使用圖像面元并行賦值和圖像疊加原理，分析星座覆蓋性，從而高效提升巨型星座多重覆蓋性的分析效率。

1 衛(wèi)星覆蓋模型

1.1 星地覆蓋模型

假設地球為正球體，衛(wèi)星對地覆蓋不考慮距離約束，則星座中的單顆衛(wèi)星對地覆蓋域在地球表面為以星下點為圓心的球冠，覆蓋模型如圖1所示。

圖1中，設為衛(wèi)星對地通訊/觀測視場半錐角;為地球半徑;為衛(wèi)星當前軌道高度，則衛(wèi)星覆蓋域對應半地心角′，以及覆蓋域球面半徑′滿足

(1)

圖1 單星對地覆蓋模型示意圖Fig.1 Single satellite to ground coverage model

1.2 同高度衛(wèi)星的星間覆蓋模型

星間覆蓋主要考慮地球曲率遮擋和大氣影響，對于同軌道高度的衛(wèi)星，單顆衛(wèi)星在同軌道高度空間天球上的覆蓋域為以該衛(wèi)星為圓心的球冠，如圖2所示。根據球面幾何關系，該球冠對應的半地心角″，以及球面半徑″滿足

(2)

式中：為大氣層厚度。

圖2 同高度星間覆蓋示意圖Fig.2 Inter satellite coverage on same height

2 PMPS覆蓋性分析方法

PMPS方法具體步驟為

以選定色值初始化給定分辨率墨卡托底圖。

單星覆蓋域邊界空間幾何關系曲線計算。針對星座中任一顆衛(wèi)星，根據式(1)或式(2) 計算其覆蓋域邊界空間幾何關系曲線。

單星覆蓋圖獲取。將空間幾何關系曲線進行墨卡托投影，并對區(qū)域內的像元進行批賦值，獲取其對地、對天覆蓋區(qū)域。

遍歷獲取覆蓋圖。按照步驟2和步驟3遍歷獲取所有衛(wèi)星的對地、對天覆蓋圖。

星座覆蓋圖獲取。通過給定透明度和色值編碼方式，依次疊加步驟4所獲取的單星覆蓋圖至墨卡托底圖。

星座覆蓋性分析。基于圖像疊加原理和并行處理算法，對星座覆蓋圖處理分析獲得覆蓋性矩陣。結合目標區(qū)域的權重矩陣和面積矩陣，獲取大規(guī)模層疊后的全星座覆蓋性分析結果。

2.1 墨卡托覆蓋圖

定義墨卡托覆蓋圖(簡稱：覆蓋圖)為星座內各衛(wèi)星在地面或天球上覆蓋域的墨卡托投影疊加圖。

2.1.1 單星覆蓋圖

圖3為單顆衛(wèi)星(=1,2,…,)對地覆蓋圖。

圖3 單星覆蓋圖Fig.3 Coverage of single satellite

由式(1)或式(2)可求得相應覆蓋模型下覆蓋域球面半徑。采用大地主題正解算法，解算以星下點為圓心的覆蓋域空間幾何關系曲線。選定色彩模型(如三原色光模型(RGB)、印刷四分色模型(CMYK)等)及色值進行圖像渲染，獲得衛(wèi)星的覆蓋圖(=1,2,…,)，為衛(wèi)星數量。編碼方式滿足

(3)

式中：為覆蓋域內像元。

2.1.2 星座覆蓋圖

以選定的透明度將各衛(wèi)星的覆蓋圖逐個疊加，其中將色值″疊加至′獲得的疊加原理如式(4)，由此得星座覆蓋圖。

=(1-)′+α″

(4)

由式(4)可知星座覆蓋重數與覆蓋圖色值一一對應。

2.1.3 各重覆蓋矩陣

式(4)可給出重覆蓋對應的色值(=1,2,…)，根據和星座覆蓋圖定義各重覆蓋矩陣(=1,2,…)，賦值方法為

(5)

2.2 星座覆蓋性分析

各重覆蓋面積及各重覆蓋率是衡量星座覆蓋性的重要指標。本文就這2類覆蓋性指標應用PMPS方法進行求解。

2.2.1 面積矩陣

定義維面積矩陣，其各元素為星座覆蓋圖中，各像元實際地理面積。由墨卡托投影特點可得矩陣：

(6)

(7)

式中：為覆蓋域球面半徑；=為圖像分辨率；為星座覆蓋圖第行像元對應地理緯度。

2.2.2 目標域權重矩陣

定義維目標域權重矩陣，中元素的位置包含于目標區(qū)域內時，為1，否則為0。滿足

(8)

式中：為第行列元素。

2.2.3 覆蓋面積與覆蓋率

定義為目標區(qū)域重覆蓋面積，滿足

(9)

式中：*為哈達瑪積。

定義為目標區(qū)域重覆蓋率，滿足

(10)

2.3 算法實現

本文算法實現包括如下幾個步驟：

初始化。

賦值。對星座參數，圖像疊加透明度、地球半徑、目標時刻等參數進行賦值。

衛(wèi)星軌道外推。

星下點求解。

單星覆蓋圖投影疊加。

計算覆蓋率及覆蓋面積輸出結果。由式(9)、式(10)計算出覆蓋面積及覆蓋率，并輸出結果。

具體算法流程如圖4所示。

圖4 PMPS流程圖Fig.4 Flow chart of PMPS

3 計算消耗分析

3.1 GPA計算消耗

GPA的計算過程主要分為3步：

1) 網格點劃分。

2) 對個網格點、顆衛(wèi)星進行逐點逐星覆蓋判定。

3) 統(tǒng)計分析結果。因此，其計算消耗′滿足

′=′+′+′

(11)

式中：′為網格劃分計算消耗；′為單星單點覆蓋性判定計算消耗；′為統(tǒng)計分析計算消耗。

由式(11)乘積項′可知，GPA網格點數與衛(wèi)星數量的乘積影響了算法計算消耗。在給定網格點數下，GPA算法中′、′兩項為定值，而′項隨著衛(wèi)星數量變化線性增長。

3.2 PMPS計算消耗

PMPS分析過程主要分為2步：

1) 對顆衛(wèi)星逐星求解覆蓋域并疊加至初始化的墨卡托底圖。

2) 統(tǒng)計分析覆蓋率。

設為單星覆蓋域空間幾何關系曲線計算時間消耗，該值為常數，為空間幾何關系曲線墨卡托投影時間消耗，該值為常數。設單星覆蓋域內像元縱坐標最大值、最小值間差值為，則基于掃描線法的覆蓋域內色彩空間編碼賦值時間消耗a()=，其中為線填充效率當量；設單星覆蓋域邊界像元數量為，則基于邊界的區(qū)域布爾疊加時間消耗b()=，其中為布爾運算效率當量。單星覆蓋域求解計算消耗″滿足

″=++a()+b()=+++

(12)

式中：、均為常數；a()、b()項與圖像分辨率有微弱關聯。

(13)

式中：(=0,1,…)為重覆蓋像元數;p(,)為圖像并行處理加速比，其理論值可由Amdahl定律給出，在數字圖像處理領域，多核環(huán)境下的并行加速比已被廣泛研究。本文后續(xù)算例中，用未引入并行計算的PMPS算法時間消耗與GPA對比分析，也即p(,)=1，因此，并行算法對效率影響本文暫不做詳細討論。

根據上述過程，PMPS計算時間消耗″滿足

″=″()+″+″()

(14)

式中：″()為初始化墨卡托底圖時間消耗，該值受圖像分辨率影響。

因此：

(15)

由式(12)可知，對于給定的覆蓋域，單星覆蓋域的計算消耗″為常數，式(14)中″項隨著衛(wèi)星數量線性增長。由式(13)可知，給定計算器核數下，覆蓋分析時間消耗″()與圖像分辨率相關，當圖像大小確定時，″()、″()為常數。式(14)可知，PMPS算法中衛(wèi)星數量和圖像分辨率對計算消耗的影響相互獨立，該算法突破了計算精度與星座衛(wèi)星數量的關聯約束。

GPA與PMPS計算消耗比值滿足式(16)，該值反映了PMPS計算效率比GPA提高的倍數。

(16)

圖5 提高倍數與衛(wèi)星數量級關系曲線Fig.5 Relation curve between multiple of increase of K and order of magnitude of number of satellites

未引入并行計算的PMPS方法不一定比GPA有明顯優(yōu)勢。

4 仿真校驗

本實驗中計算機CPU為Intel(R)Core(TM)i7-6700 K CPU@4.00 GHz,運行內存為4 GB,顯卡為NVIDIA GeForce GTX 1060 6 GB，操作系統(tǒng)為Windows 10 64 Bit，開發(fā)語言為Python3.7。用未引入并行計算的PMPS算法時間消耗與GPA進行對比分析。

由于傳統(tǒng)網格法計算精度隨網格數量增加逐漸提高，校驗前先確定給定計算精度下GPA網格數量及PMPS圖像分辨率。該過程基于對構型為1 584/24/1的Starlink子星座計算結果完成。衛(wèi)星天線掃描半錐角設為40°，軌道傾角為53°，軌道高度為550 km。

確定計算精度后，算例1對現有典型星座或其子星座進行仿真實驗，研究未引入并行計算的PMPS算法與GPA算法間計算誤差及計算消耗；算例2在Starlink子星座基礎上進行構型變化，研究衛(wèi)星總數量恒定、星座構型變化下2種算法間的計算誤差及計算消耗；算例3采用Starlink的種子衛(wèi)星參數進行衛(wèi)星數量變化，研究衛(wèi)星數量變化下2種算法間的計算誤差及計算消耗。本仿真校驗中均采用Walker星座，構型為，其中為衛(wèi)星數量，為軌道平面數，為相位因子。

4.1 計算精度確定

本文采用移動方差法確定GPA計算結果趨于穩(wěn)定時的網格數。取移動方差窗口步長為5，直至方差小于0.1認為計算結果趨于穩(wěn)定。結果見表1。并對選定網格數進行網格無關性驗證，結果見表2。在確定GPA計算精度后，逐漸提高PMPS中覆蓋圖分辨率，將計算結果和GPA結果進行對比，兩者結果相對誤差在5%以內時，認為滿足精度要求，結果見表3。

表1 移動方差Table 1 Moving variance

續(xù)表1

表2 網格無關性驗證Table 2 Grid independence verification

表3 PMPS與GPA相對誤差Table 3 Relative error of PMPS and GPA

4.2 仿真算例

對現有星座的計算誤差及計算消耗

在現有典型星座或其子星座覆蓋性分析算例中，僅考慮物理遮擋因素下的全球覆蓋性。星座參數見表4。未引入并行計算的PMPS與GPA分別計算星座對地全球1～5重覆蓋率，算法間誤差如圖6所示，圖中(=1,2,3,4,5)表示重覆蓋率誤差。2種算法計算消耗見圖7。

由圖6可知，2種算法計算誤差在2%以內。由圖7可以看出在不同星座下，未引入并行計算的PMPS計算消耗均在400～500 s范圍內波動。本算例中GPA對于數量在10量級內的星座，計算消耗相對較低，衛(wèi)星數量在10量級時，GPA計算消耗高于未引入并行計算的PMPS，這是由于未引入并行計算的PMPS的消耗項″()與圖像分辨率呈正相關，在較大圖像分辨率下，″()較大。

星座構型變化的計算誤差及計算消耗

表4 現有典型星座或子星座參數

圖6 不同星座未引入并行計算的PMPS與GPA間的計算誤差Fig.6 Calculation errors between PMPS without parallel computing and GPA for different constellations

圖7 不同星座下未引入并行計算的PMPS與GPA計算消耗Fig.7 Time consumption of PMPS without parallel computing and GPA for different constellations

本算例中取衛(wèi)星總數量為1 584顆，位于高度550 km的圓軌道，傳感器錐角為40°。對表5所示16個構型的星座進行對地全球1～5重覆蓋率計算。2種算法間的計算誤差見圖8，2種算法計算消耗見圖9。

由圖8可知2種算法的計算誤差保持在1%內。由圖9可知GPA算法計算消耗在其均值2 492.57 s的2.73%范圍內波動，未引入并行計算的PMPS算法計算消耗在其均值514.13 s的2.72%范圍內波動，因此可基本排除星座構型對計算消耗的影響。

表5 星座構型Table 5 Constellation configurations

圖8 不同構型下未引入并行計算的PMPS與GPA間的計算誤差Fig.8 Calculation errors between PMPS without parallel computing and GPA for different constellation configurations

圖9 不同構型下未引入并行計算的PMPS與GPA計算消耗Fig.9 Time consumption of PMPS without parallel computing and GPA for different configurations

衛(wèi)星數量變化下的計算誤差及計算消耗

本算例中各Walker星座均采用24軌道面，星座衛(wèi)星數量由288顆以步長72依次增至1 584顆。具體星座構型見表6。分別用2種算法計算各星座對地全球1～5重覆蓋率，計算結果誤差如圖10所示，計算消耗如圖11所示。

由圖10可知，2種算法計算各星座1～5重覆蓋率誤差在2%以內。由圖11可知2種算法符合線性增長趨勢，對數據進行線性擬合可知，在該計算精度下，GPA計算消耗隨衛(wèi)星數量的增長率為1.643 4，未引入并行計算的PMPS計算消耗隨衛(wèi)星數量的增長率僅為0.027 6，因此未引入并行計算的PMPS計算消耗對衛(wèi)星數量的敏感度明顯低于GPA。

表6 衛(wèi)星數量及其星座構型

圖10 不同衛(wèi)星數量下未引入并行計算的PMPS與GPA間的計算誤差Fig.10 Calculation errors between PMPS without parallel computing and GPA for different number of satellites

圖11 不同衛(wèi)星數量下未引入并行計算的PMPS與GPA計算消耗比較Fig.11 Time consumption of PMPS without parallel computing and GPA for different number of satellites

圖12 給定精度下未引入并行計算的PMPS相對GPA提高倍數與衛(wèi)星數量級關系曲線Fig.12 Relation curves between multiple of increase and order of magnitude of number of satellites of PMPS without parallel computing relative to that of GPA at given accuracy

5 結 論

提出一種高效的星座覆蓋性分析方法PMPS，通過計算效率的理論分析，及現有星座和構型參數變化下的星座計算分析得出以下結論：

1) 給定計算精度下，隨著衛(wèi)星數量的增加，GPA與未引入并行計算的PMPS計算消耗比值逐漸增加，并最終趨于極限值，由式(16)該極限值和GPA計算精度選取呈正相關。

2) 在GPA網格數為10量級、2種算法計算誤差小于2%時，未引入并行計算的PMPS計算效率比GPA最大提高1～2個量級，當GPA網格數量更大時，如10量級，該計算效率將提升2～3個量級。

3) 后續(xù)選擇合適的處理器核數、引入圖像壓縮技術，將降低PMPS計算精度對計算消耗的影響系數，有效減小式(14)中″()項。使PMPS計算效率相較于GPA在較小量級星座下便開始提高。

4) 本文工作可以為其他覆蓋計算方法的設計提供參考，為解決星座設計問題和星座調度問題提供依據。