地面對量子衛星信號捕獲及粗跟蹤過程的仿真研究*

叢 爽 汪海倫 陳 鼎

1.中國科學技術大學自動化系，合肥2300272. 北京衛星信息工程研究所，天地一體化信息技術國家重點實驗室，北京100086

量子定位是以量子光為信息載體，在量子衛星和地面站之間進行數據傳輸的一種定位方法。量子光的糾纏特性，決定了量子定位與傳統定位相比具有定位精度高、保密性好等優勢，但同時也存在精確對準困難和保持穩定鏈路難度大等問題。ATP(Acquisition, Tracking and Pointing)系統采用粗精跟蹤的復合嵌套，通過對量子衛星發射信標光的捕獲、瞄準和跟蹤，可以使地面站與量子衛星間建立穩定的信標光鏈路和通信光鏈路，解決精確對準困難和保持穩定鏈路難度大的問題[1]。其中，粗跟蹤系統負責在初始階段掃描和捕獲大范圍的信標光，并引導信標光進入精跟蹤視場；精跟蹤系統負責在粗跟蹤精度的基礎上減少跟蹤誤差和補償由衛星運動引起的超前瞄準誤差。

ATP系統的正常工作需要以量子衛星軌道信息作為其控制器的輸入信號，控制器驅動電機使光學天線的視軸始終指向量子衛星，維持量子光鏈路的穩定。在ATP系統的工作中，量子衛星的軌道信息一般是通過星歷表計算或者GPS得到[2]，在此基礎上對已知的衛星位置進行坐標轉換，才可以作為ATP系統控制器的輸入信號。此外，還需根據軌道信息計算得到的通過時間決定ATP系統的開始工作時刻。由美國分析圖形有限公司(Analytical Graphics Inc, AGI)研制開發的STK(Satellite Tool Kit)軟件是國際航天領域中先進的系統分析軟件，具有強大的計算分析和仿真顯示功能[3]，可以根據已知的軌道參數進行仿真，直接給出ATP系統控制器所需的輸入信號，但是其中的變換過程沒有詳細給出。

本文在完成衛星運行軌道設計的基礎上，首先計算出衛星通過地面站上方可視范圍的通過時間，然后推導了衛星軌道信息在WGS-84坐標系下的位置參數，變換到地球東北天坐標系下，再變換到地面站載荷設備坐標系下，計算出地面接收端的方位角和俯仰角隨時間變化曲線(即ATP系統控制器的輸入信號)的計算過程，之后以地面站位于安徽省合肥市為例，給出3次左邊轉換的計算結果，并在STK中以相同的參數對衛星軌道進行軟件仿真，給出三維和二維衛星通過地面站上方的仿真動畫，最后將仿真動畫與控制器輸入信號、控制器捕獲和粗跟蹤階段的仿真結果在Matlab的GUI中直觀地顯示出來。

1 軌道參數確定及時間計算

唯一地確定1條人造地球衛星的軌道一般需要6個參數：軌道半長軸，偏心率，軌道傾角，升交點赤經，近地點角距和真近點角。在ATP系統的工作中，量子衛星的軌道信息一般是由衛星測控部門通過軌道預報算法得到[3]，本文我們是在設定衛星軌道為圓軌道，并利用衛星軌道高度，以及衛星軌道平面滿足經過地面上空區域某點的條件，在MATLAB環境下完成衛星軌道計算的。所以，所采用的參數為確定衛星軌道平面的2個相交向量(6個參數)，加上衛星軌道高度，一共7個參數，其中，軌道高度r為衛星繞地球運行的軌道距離地球表面的高度，一般用近地點與遠地點均值表示。參考“墨子號”量子衛星的實際軌道的長半軸為584km，短半軸為488km[4]，為了簡化模型，量子定位衛星的軌道高度設計為r=500km的圓軌道。

下面我們還需要衛星經過地面端上方區域的時間。圖1為衛星軌道經過地面端正上方的過程及其時間，其中，圖1(a)中實線代表地球表面，虛線代表衛星運動軌道，水平線上方部分的虛線圓弧為衛星經過地面的可見區域；圖1(b)是衛星經過地面端上空區域情況，其中，AB段為星地ATP的捕獲做好機動準備，主要通過控制地面端方位軸俯仰軸的轉動與衛星自身姿態實現；在BC段完成衛星對地面端的捕獲的單向瞄準過程，實施方法是由地面端發射信標光(Transmitting beacon light)，讓衛星端接收信標光，通過這一過程來實現該功能；然后在CD段完成星地間的雙向瞄準(Mutual pointing)，主要通過衛星端向地面端發射信標光，地面端捕獲信標光來實現；DE段為測試時間(Test time)，進行基于量子光通信的星地間量子定位實驗的測試。衛星經過地面端上方區域時間為衛星過境的總長AE段時間。

衛星沿其軌道環繞一整圈所需要的時間，稱之為軌道周期T：它可以由軌道周期公式獲得[5]：

(1)

其中，R為地球半徑，取6378km；r=500km，GM是地球引力常數，取3.986×1014m3/s2。

對于經過地面正上方的勻速圓周運動的衛星軌道，其過境時間總長Ta與周期T成正比：

Ta/T=2γ/2π

(2)

其中，γ為地面端地平線垂直的地球半徑與地面端地平線和軌道交點的軌道半徑之間的夾角，它可以根據圖1(a)中所顯示出的反三角函數關系計算出來：

γ=arccos(r/(r+R))=
arccos(6378/(6378+500))=22.0°

(3)

根據圖1(a)中所表示出的衛星運動軌道(整個虛線圓弧)與衛星經過地面的可見區域(水平線上方部分的虛線圓弧)之比，求出衛星過境時間總長Ta為Ta=γT/π=692.4s。以DE段為130°為例，可以計算出星地間量子衛星經過地面端的可見時間為：TDE=(130°/180°)Ta=500.1s。

圖1 衛星軌道參數

2 方位角和俯仰角推導

2.1 坐標定義及計算過程

GPS定位系統的地心坐標系是WGS-84坐標系，該坐標系是右手坐標系，該坐標系以地球質心作為坐標原點，x軸為指向零子午面和赤道交點的方向，z軸為指向協議地極方向。地球東北天坐標系也是右手坐標系，在該坐標系中，原點是地面端中心，x軸和y軸在當地水平面上，x軸在當地水平面上指向正東方向(東)，y軸在當地水平面上指向北方向(北)，z軸指向上天頂(天)，這也是東北天坐標系命名的由來。地面端載荷設備坐標系根據地面端的載荷設備安裝角度來定義，若載荷設備坐標系的x軸與地球東北天坐標系的x軸重合，地面端載荷坐標系與東北天坐標系只有一個繞x軸旋轉角度的差別。

在星地量子定位中，一般是采用ATP中二維轉臺的方位角與俯仰角來控制光學天線的轉動，它們是基于載荷設備坐標系而定義的，以出射光沿-y軸方向定義為基準零位，則方位角繞z軸逆時針轉動的方向為正，俯仰角繞x軸指向天的方向為正。所以有必要將由衛星在WGS-84坐標系中的位置矢量，通過坐標系轉化為載荷設備坐標系中的位置矢量，然后通過解耦獲得二維轉臺的方位角與俯仰角。

圖2給出了計算地面端方位角、俯仰角的流程。下面將根據圖2依次推導：1)將WGS-84坐標下的位置矢量轉換到地球東北天坐標系下；2)將由1)得到的位置矢量轉換到地面端載荷設備坐標系下；3)在地面端載荷設備坐標的情況下將三維空間函數解耦為方位軸與俯仰軸的關于時間的函數；4)最終根據具體條件計算結果。

圖2 方位角與俯仰角求取流程圖

設在WGS-84坐標系給出的衛星軌道位置的經度L、緯度B和高度H，則在WGS-84坐標系中直角坐標為[6]：

(4)

1)WGS-84坐標系到地球東北天坐標系轉換

該過程主要分為4步驟進行：a)WGS-84坐標系繞z軸旋轉L；b)繞y軸旋轉-B；c)繞y軸旋轉90°；d)繞z軸旋轉90°。該坐標變換過程用矩陣的乘積可表示為式(5)[6]，其中L和B分別表示經度和緯度：

(5)

2)東北天坐標系到地面載荷坐標系轉換矩陣

令東北天坐標系轉換到地面端載荷設備坐標系的轉換矩陣為Med，在該矩陣中，地面載荷設備的安裝角度與其密切相關，假若重合的一軸為x軸，Med如式(6)所示。

(6)

其中，?為將東北天坐標系繞x軸旋轉后與地面端載荷設備坐標系重合所需的角度。

3)方位角和俯仰角的求取

在WGS-84坐標系下，令OSw為衛星的位置矢量，OF為地面端的位置矢量。令GS為經過坐標變換后的衛星在地面端載荷設備坐標系下的位置矢量，其表達式如下所示，并可以由式(5)和(6)帶入得出：

GS=Med·Mdw(OSw-OF)

(7)

然后將衛星在地面端載荷設備坐標系下的位置矢量GS解耦為方位角As和俯仰角Es[6]：

(8)

其中，xGS,yGS,zGS為位置矢量GS的x,y,z坐標值。

2.2 計算結果及分析

在衛星距地面高度R=500km，地球半徑為r=6378km已知的情況下，地面端經緯度表示為(L，B)，地面端高度為H。假定地面端位于安徽省合肥市，其經緯度為(117.27°，31.86°)，地面端高度為0.5km。衛星軌道即為在WGS-84坐標系中，以原點為圓心，過地面端上方500km一點的空間圓。由已知條件，通過公式(4)可以算出地面端在WGS-84坐標系中的坐標，選擇與OF垂直的向量n，并求與n垂直且相互垂直的2個單位向量u和v。根據空間圓的參數方程，當θ∈(-π,π)，衛星軌道可以由參數方程表示為[6]：

(9)

計算得到的衛星軌道(x,y,z)w就是衛星的位置矢量OSw，再經過式(5)的Mdw和式(6)Med的坐標轉換，其中，坐標轉換矩陣Med中角度?由在地面載荷設備實際安裝后進行標定，此處取?=12°。最終可以根據式(8)得到方位角As和俯仰角Es，在Matlab環境下實現的計算結果如圖3(a)所示，其中，實線為方位角As，點劃線為俯仰角Es，整個輸入信號長度為692s，虛線框為實際可以進行星地間量子定位實驗的時間段，大約為500s。圖3(b)為方位角與俯仰角之間的關系。從圖3中可以看出方位角的變化范圍為[66°，247°]，俯仰角的變化分為 [0°，78°]。最大俯仰角沒有達到預計的90°，主要原因在于本文采用的俯仰角計算方法存在局限性，根據公式(8)，俯仰角由第2個式子求反正切得到，理論上只有反正切的值無窮大時，才能得到90°的俯仰角，實際中并不能得到。此外，在編寫運算過程中，有一部分長度參數是以km為單位，經過多次運算會存在一定誤差，影響了最后俯仰角的結果。

2.3 可覆蓋空間的仿真

為了得到衛星與地面端相對位置的2D和3D視圖，需要在STK中創建一個與上述例子相同的STK場景(Scenario)并運行仿真，該場景包含一個地面站(Facility)、一個衛星(Satellite)及其對應的傳感器(Sensor)[7]，需要進行以下步驟設計：

1)創建地面站：合肥地面站具體參數為經度(Longitude)117.27°，緯度(Latitude)31.86°，高度(Altitude)0.5km。

2)創建衛星：由于量子定位衛星的軌道高度設計為500km，在STK軟件里衛星軌道參數的設定中，選擇從地球表面到橢圓的最大值(Apogee Altitude)和從地球表面到橢圓的最小值(Perigee Altitude)，都設為500km。衛星軌道參數中的軌道傾角設定為109°，根據衛星與地面站之間的訪問關系，2018年7月15日上午九時許兩者時間可以產生訪問。根據訪問可持續的時間，本文的場景運行時間設置為2018年7月15日9：19到9：31。

3)創建傳感器：分別在地面站和衛星上創建傳感器，傳感器的類型(Sensor Type)選擇圓錐形(Simple Conic)，地面站上傳感器的圓錐角(Cone Angle)設置為60°，衛星上傳感器的圓錐角設置為45°。

圖3 地面端的方位角與俯仰

圖4為采用STK軟件設計的衛星通過地面端上方路徑的仿真結果在同一時刻的截圖，其中，(a)為三維動畫仿真截圖，(b)為二維動畫仿真截圖。圖4(a)中一條圍繞地球一圈細線是衛星繞地球運動的軌道，軌道上的實心方點是衛星目前所在位置，軌道下面的方點是地面站所在位置，以衛星為頂點的圓錐面是衛星端傳感器圓錐角設置為45°時可覆蓋空間的邊界，在圖4(b)的二維視圖中，是以衛星為圓心的扇形顯示，其中大的(紅色)橢圓是地面端傳感器圓錐角設置為60°時可覆蓋空間的邊界。衛星端傳感器的可覆蓋區域和地面端傳感器的可覆蓋區域重疊時，STK軟件中的二維和三維仿真界面中才會顯示兩端傳感器的可覆蓋空間的邊界。換句話說，只有當兩者具有交集的部分時，衛星與地面端之間才能夠進行鏈路的建立。

3 仿真演示平臺設計

ATP系統的工作過程可分為捕獲、粗跟蹤和精跟蹤3個過程，其中捕獲過程是根據GPS得到地面光學天線初始指向的方位角和俯仰角，然后驅動二維轉臺，通過掃描過程對準信標光，并實施捕獲動作，當雙方互相捕獲到對方的信標光時，捕獲過程結束；粗跟蹤過程是根據接收到的信標光位置不斷調整二維轉臺，保持星地兩者在相對運動中維持粗跟蹤精度穩定[8]?；贕UI[9]對量子定位系統仿真演示平臺中捕獲和粗跟蹤部分進行設計，捕獲和粗跟蹤部分所涉及的所有參數選擇及輸入、輸出信號顯示被設計到輸入信號區、捕獲和粗跟蹤區2個區域中，得到的動畫仿真演示平臺如圖5所示?？梢栽谳斎胄盘枀^域根據所選擇初始指向誤差的坐標，展示衛星與地面端相對關系的二維、三維動畫和經過坐標轉換后的方位軸和俯仰軸的角度變化曲線；在捕獲和粗跟蹤部分選擇開始捕獲的位置，獲得相應的捕獲與粗跟蹤的仿真結果。

圖4 仿真動畫

3.1 初始指向誤差與開始捕獲位置的設置

圖5中的初始指向誤差是光學天線在已知衛星運動軌道的情況下，對衛星進行初始指向后光學天線的指向與衛星之間的x-y方向的誤差，我們在由初始預指向精度0.5(換算成直徑為8.7mrad的不確定區域內，共設置了6種初始指向誤差[-0.8889,-0.9851]，[-2.2558,1.5323]，[1.1981,3.9819]，[-3.2734,-1.8682]，[1.6977,-3.7585]和[3.0300,2.7300]；其選擇的控件有：1個坐標軸控件，6個單選按鈕和1個按鈕。6個單選按鈕用來選擇不同的初始指向誤差，按鈕用來在捕獲和粗跟蹤區域內左上角的坐標軸控件中顯示期望的方位軸信號和俯仰軸信號的圖像。

圖5 基于GUI的動畫仿真演示平臺

與開始捕獲位置選擇相關的控件位于捕獲和粗跟蹤區域中，包含：1個坐標軸控件(axes1)，3個單選按鈕和1個按鈕。3個單選按鈕用來選擇開始捕獲的位置，我們設定有3個位置：位置1是從衛星剛剛進入可見區域，即從輸入信號的時刻為0；位置2是衛星經過地面端上方時，即從輸入信號總長的1/2位置開始捕獲；位置3是衛星進入可見區域后，且距離到地面端上方一半的位置，即從輸入信號總長的1/4位置開始捕獲。捕獲按鈕中需添加捕獲階段仿真的回調函數。

根據選定的掃描路徑式編寫掃描路徑的函數，設不確定區域為8.7mrad，粗跟蹤探測器視場為3mrad，粗跟蹤探測器的幀頻為100fps，掃描步長為倍的粗跟蹤探測器半徑，約1.59mrad，函數的終止條件是目標與粗跟蹤探測器的坐標歐式距離小于粗跟蹤探測器的視場半徑[8]。捕獲部分的仿真結果在捕獲按鈕按下后，會動態地顯示在坐標軸控件(axes2)中，捕獲仿真結束后，3個靜態文本中會顯示捕獲仿真的結果，即所需要的捕獲步長和此時方位軸和俯仰軸的位移。捕獲部分的仿真結果同時受初始指向誤差和開始捕獲位置2個選擇的影響，其中，初始指向誤差影響捕獲仿真的步長，開始捕獲位置影響捕獲仿真的位移偏移量。

3.2 粗跟蹤部分的演示

粗跟蹤部分采用2個電流-速度-位置三閉環粗跟蹤控制器，分別控制二維轉臺的方位軸和俯仰軸，其期望信號即2.2節計算得出的方位軸和俯仰軸隨時間變化的函數。每個三閉環粗跟蹤控制器都有5個參數：速度環比例系數Kvp、速度環積分系數Kvi、位置環比例系數Kpp、位置環積分系數Kpi、位置環微分系數Kpd。根據文獻[10]的調參結果，我們取方位軸控制器參數為Kvp=1.5、Kvi=0.15、Kpp=45、Kpi=7和Kpd=5，俯仰軸的控制器參數為Kvp=0.01、Kvi=0.1、Kpp=10、Kpi=1.6和Kpd=0.01。

與粗跟蹤部分相關的控件有：2個坐標軸控件(axes3和axes4)和1個按鈕，如圖5所示，其中，2個坐標軸控件分別用來顯示粗跟蹤過程方位軸和俯仰軸的跟蹤誤差及仿真結果，粗跟蹤按鈕中需添加粗跟蹤階段仿真的回調函數，需要注意的是粗跟蹤按鈕必須在單擊捕獲按鈕后，并且等捕獲仿真結果在坐標軸控件中完全顯示后才能單擊，否則由于回調函數中初始化操作的存在，會打亂坐標軸控件中的仿真結果。同時，捕獲階段仿真結束后，為了直觀顯示開始捕獲位置的選擇，會用1條綠色豎線在輸入信號的坐標軸控件(axes1)中標記捕獲完成時的位置，由于捕獲時間相對于整個輸入信號時間來說非常短，平均只占萬分之一，因此捕獲完成時的位置可以近似等于開始捕獲的位置。粗跟蹤階段的仿真是建立在捕獲階段完成的基礎上，因此粗跟蹤仿真時間是從捕獲完成時刻開始，左下角的坐標軸控件(axes3)和右下角的坐標軸控件(axes4)會同時同步顯示粗跟蹤階段方位軸和俯仰軸在時間序列上的誤差和在空間序列上的誤差，同時在輸入信號的坐標控件中顯示第2條綠色豎線，動態標記粗跟蹤仿真進行的時刻。

根據上述設計，我們得到圖5右邊的4幅實驗結果圖，圖5中間兩幅圖從上到下分別為在選擇了初始指向誤差的坐標情況下所對應的輸入信號及其粗跟蹤后的跟蹤誤差，其中紅色實線為俯仰軸信號，藍色虛線為方位軸信號。圖5右邊兩幅圖從上到下分別為信號的捕獲過程，以及粗跟蹤誤差平面指示圖，其中，上圖中的中心紅色小虛線圓圈為捕獲的開始位置，藍色大虛線圓圈為直徑為0.5°(8.7mrad)的不確定區域，小的實線圓圈為捕獲掃描過程，每一個小圓圈代表直徑為3mrad的粗跟蹤探測器視場的一次捕獲，與紅色點(目標點)相交的粗藍色圓圈為捕獲到信號并開始粗跟蹤。圖5右下圖中的紅色大圓圈為直徑為0.5mrad的粗跟蹤性能范圍，中間的黑色軌跡為粗跟蹤過程中的信號變化軌跡?？梢钥闯?，在所設計的粗跟蹤控制系統參數下，粗跟蹤控制完全能夠達到0.5mrad的粗跟蹤精度性能。

4 結論

設計了基于GUI的動畫仿真演示平臺中的捕獲和粗跟蹤部分，在Matlab環境中實現了地面端接收衛星發射信號的方位角和俯仰角的計算，并將其結果作為ATP系統的期望輸入信號，完成了在6種初始指向誤差和3種開始捕獲時刻條件下，ATP系統捕獲和粗跟蹤部分的仿真。同時，在STK軟件中對衛星在軌運動進行三維動畫仿真，將仿真動畫集成在基于GUI的動畫仿真演示平臺中。