基于光學測角數據的風云四號同步軌道衛星精密定軌

劉思語，黃勇，毛銀盾，賈耀紅，魯文強,4，黃乘利,3，鄭景輝，楊鵬

1. 中國科學院上海天文臺，上海 200030 2. 中國科學院大學 天文與空間科學學院，北京 100049 3. 上?？萍即髮W 物質學院，上海 201210 4. 國家衛星氣象中心，北京 100081

1 引言

風云四號氣象衛星于2016年12月11日由西昌衛星發射中心成功發射系列第一顆衛星(以下稱FY-4A)并于2017年9月25日正式交付用戶使用，是中國第二代靜止軌道氣象衛星，用于在地球同步軌道帶上進行高效、穩定、連續的氣象觀測[1]。

FY-4A目前采用無線電測距技術實現對衛星的精密定軌，歲差章動模型參考IAU1976，N體攝動模型參考DE405歷表，采用數值求解二階微分方程組的KSG定步長線性多步積分器求解動力學方程，重疊1 h軌道的內符合精度優于20 m。FY-4A衛星定點于東經99°，地面五個跟蹤站(北京，佳木斯，廣州，烏魯木齊，騰沖)同時進行觀測，數據采樣率為1 s，觀測模式為每小時或每半小時連續觀測10 min，每日軌控一次，獲取得到S波段與L波段的測距信號，對測距數據經過系統差標校后，利用雙頻改正方法進行電離層延遲改正，每個測站配備自動氣象站，為大氣折射修正提供氣壓、溫度、濕度等參數。

對于大部分地球同步軌道衛星(GEO)定軌，國內現階段除了光學測角手段之外，還包括無線電測距、激光測距(SLR)、雷達等技術[2]。宋葉志等基于多臺站雙頻雙程測距數據計算風云四號氣象衛星軌道，通過軌道重疊分析，在非變軌時期獲得20 m的定位精度[3]；劉凱等采用一發多收的轉發式測定軌方法測量北斗導航衛星的軌道，精度達到10 m[4]；卿蕓等利用星載GNSS和SLR數據對北斗C08、C10地球同步軌道衛星進行聯合定軌，定軌精度優于0.1 m[5]；Zongbo HUYAN等基于星間鏈路測量對北斗導航衛星進行定軌，48 h內軌道預報精度優于2.4 m[6]；張志斌等使用口徑1 m的“微型天線”在上海、都勻和烏魯木齊三個測站組網，監測“亞太六號”地球同步軌道衛星，軌道精度可達百米量級[7]；Yong Huang等利用地基VLBI數據對地球同步軌道衛星定軌，VLBI時延精度可達3.6 ns，定軌精度優于20 m[8]。

地基光學觀測手段相比上述測量技術具有運行成本低、結果更直觀等優勢，并且作為一種無源觀測手段，在探測非合作目標時擁有更廣泛的應用范圍[9]。目前地面的光學觀測設備主要成像終端為CCD或CMOS相機，通過照相天體測量的方式獲取目標星的位置。照相天體測量技術需要預知參考星的理想坐標、參考星和目標星在CCD或CMOS像面上的量度坐標，利用每一組參考星的理想坐標和量度坐標構建底片模型，再通過底片模型和目標星的量度坐標反解出目標星的理想坐標和天球坐標。GEO相對地面測站幾乎靜止不動，受到地球自轉影響，地球同步軌道衛星與背景恒星存在運動差異(大約15＂/s)，因此相機采用凝視模式觀測并獲得GEO圓形星象的同時，背景恒星的星象圖像會拖長，無法精確確定參考星的位置中心。于涌等利用CCD漂移掃描技術解決星象拖長的問題[10]。望遠鏡對目標精密跟蹤的情況下，星像固定于CCD視場。星像所產生的電荷在同一像素區域內逐漸積累。當積累的能量明顯超過周圍的背景噪聲時，信噪比高的星象會顯現出來。CCD漂移掃描技術實現電荷在像素間的行轉移，當電荷轉移速度與目標天體視運動速度在CCD像面上的投影同步且目標運動方向與電荷轉移方向一致時，能夠在固定望遠鏡指向下實現對運動目標的電荷跟蹤，得到圓形星象。嚴丹等利用相似的方式觀測某顆GEO和木衛六，分別在赤經赤緯上取得0.32＂、0.28＂和0.100＂、0.105＂的觀測精度[11]。2011年西班牙皇家天文臺基于光學測角數據對H1D GEO衛星精密定軌，取得100 m的外符合精度[12]。

本文利用光學測角數據對FY-4A數據進行定軌分析，分析了僅利用測角數據定軌以及利用測距測角數據聯合定軌的精度，并評估了僅利用測距數據的實時軌道精度。

2 觀測實驗

2019年12月27～31日，上海天文臺GEO全球聯合觀測網50 cm望遠鏡(以下稱50 cm望遠鏡)對FY-4A進行了光學觀測，觀測樣例見圖1和圖2。利用CCD相機觀測得到FY-4A的測角數據，直接測量量為目標星的量度坐標，再轉換為定軌所需要的赤經赤緯或高度角與方位角。

圖1 利用CCD漂移掃描技術觀測恒星背景實例Fig.1 An example of CCD drift scanning technology for sidereal stars

50 cm望遠鏡安裝了支持時延積分功能的Alta-F9000相機和GPS/北斗雙模時統，具備CCD漂移掃描功能具體參數見表1。2019年12月27～31日該望遠鏡對FY-4A進行觀測，CCD相機每一輪觀測拍攝7幀圖像，首末2幀以漂移掃描模式拍攝背景恒星，中間5幀以凝視模式跟蹤FY-4A。相機快門開啟時時間系統記錄開始時刻，快門關閉時記錄結束時刻。每2幀圖像之間的時間間隔約為6.7 s，其中5 s為曝光時間，1.7 s為CCD電荷讀取時間。5 s的曝光時間能保證望遠鏡可以觀測15星等的GEO目標。光學測角要求測站處于陰影且被測目標處于可見光照射，保證良好天光背景，因此觀測時間多集中于天氣狀況良好的晴朗夜間[13]，具體觀測弧段見表2。

圖2 利用凝視模式觀測GEO衛星實例Fig.2 An example of stare-mode for GEO satellites

表1 望遠鏡參數及測站位置

表2 光學測角弧段

初始圖像經過中值濾波去除高頻信號，保留低頻信號，抑制背景的不均勻。再通過連通域星象法檢測圖像中星象，若某星象的信噪比小于設置門限，則去除該星象。同時基于二維修正矩定心法消除背景噪聲的影響，計算星象中心。最后利用相鄰幀圖像差分和航跡關聯法識別GEO衛星[14]。

3 測量模型

3.1 照片天體測量與底片模型

圖3 赤經赤緯與理想坐標幾何關系Fig.3 The geometric relationship between ideal coordinate and ascension,declination

照相天體測量技術基于底片上目標星的量度坐標推算其理想坐標，并計算其赤經赤緯，赤經赤緯與理想坐標的幾何關系見圖3。根據觀測結果所計算的赤經和赤緯均在站心赤道坐標系下表示。

在處于焦平面上的照相底片取一點C，以C為原點，以焦平面為坐標平面構建笛卡爾坐標系C-ξη，η為赤緯圈在底片上的投影，以此坐標系為理想坐標系。CO軸為望遠鏡的光軸，C'(α0,δ0)為CO軸延長線與天球的交點,其中α0和δ0為該點的赤經和赤緯，|CO|為焦距F?，F有一顆目標星在焦平面上的投影點為S,其對應的在天球上的點位為S'，P為北極所在位置。若設焦距F為單位變量，在三角形OSC中則有：

(1)

式中：φ為OC和OS夾角；θ為CS與η軸夾角。在球面三角形C'S'P中，存在以下轉換關系：

(2)

聯立(1)(2)有：

(3)

量度坐標系設立于CCD圖像，其坐標平面為CCD圖像平面，即CCD芯片所在平面。第一軸與第二軸為矩形像元的兩條邊，第三軸垂直于坐標平面。理想坐標系的坐標平面即為望遠鏡焦平面，X軸和Y軸分別平行于赤經圈與赤緯圈，Z軸垂直于坐標平面。設備安裝中無法保證CCD芯片與焦平面平行，導致整套光學系統不共軸。理想坐標系的原點為望遠鏡的光心，量度坐標系的原點為相關量度設備在照相底片上的位置中心，二者無法完全一致。不同的望遠鏡焦距和底片尺寸會導致理想坐標系和量度坐標系的坐標軸比例尺無法嚴格相等。CCD像元由于制作工藝的因素可能導致像元不是嚴格的矩形，量度坐標系的兩軸無法嚴格正交，可能與理想坐標系的兩軸不平行，除此之外量度坐標系受到光學系統像差等影響，因此通過多項式建立量度坐標(x,y)和理想坐標(ξ,η)的關系[15]：

(4)

底片模型參數aij和bij通過背景恒星的理想坐標和量度坐標建立觀測方程，基于最小二乘法求解獲得。由(4)可建立赤經赤緯與量度坐標系的轉換關系(5)[16]：

(5)

本文定軌所使用觀測資料為高度角與方位角，赤經赤緯與高度角方位角存在以下轉換關系，A和E為目標星的方位角和高度角，S為測站所在恒星時。

3.2 高度角與方位角

方位角和高度角是基于測站坐標系下的概念，測站坐標系的原點位于觀測站站心，即測量設備跟蹤天線的旋轉中心。站心所在位置的地平面為基本平面，主方向由站心指向正北方向。笛卡爾系下的測站坐標系，YS軸在基本平面內指向東方，XS軸指向主方向，ZS軸垂直于基本平面指向上方。

現設測站坐標系下有一目標位置矢量為rtp(xtp,ytp,ztp),方位角A表現為目標位置矢量在基本平面上的投影與Xs軸的夾角[17]：

方位角A對衛星狀態的偏導數：

(6)

式中：r為衛星在J2000.0天球坐標系中的位矢；rb為衛星在地固坐標系中的位矢。式(6)具體展開為：

(7)

式中：HG為J2000.0天球坐標系到地固坐標系的轉換矩陣，M為地固坐標系到測站坐標系的轉換矩陣。HG包括歲差、章動、自轉和極移的修正，λ和φ分別為站心的經度和大地緯度,式(7)中M為：

高度角E表現為目標位置矢量與Zs軸的夾角：

高度角E對衛星狀態參數的偏導數：

3.3 測距

測距型數據獲取的是測站與衛星相對距離關系，可表示為：

式中：ρ為測站相位中心和星載天線相位中心的距離；ρu為上行信號傳播距離；ρd為下行信號傳播距離。ρu和ρd具體可表示為：

式中：t為測站接收下行信號的時刻；Δt1為上行光行時；Δt2為下行光行時。r(t-Δt2)為衛星發送下行信號時，衛星的位置矢量。R(t-Δt1-Δt2)為測站發送上行信號時，測站的位置矢量。R(t)為測站接受下行信號時的位置矢量。Δρ是誤差修正，包括相對論修正、電離層修正、大氣折射修正等。εu和εd分別為上行信號傳播距離和下行信號傳播距離的測量誤差。上行光行時Δt1和下行光行時Δt2無法直接獲得，可先迭代計算下行光行時，推導得到衛星發送下行信號的時刻t-Δt2，再經過一次計算上行光行時，推導得到發送上行信號的時刻t-Δt1-Δt2。具體的迭代過程如下：

3.4 定軌原理

在對FY-4A進行光學觀測的過程中。假設在ti時刻獲得一組觀測量Yi:

(8)

設x(t)=X(t)-X*(t)，由(8)式得：

4 數據分析及小結

4.1 測角數據定軌及精度分析

基于50 cm望遠鏡的測角數據，對2019年12月28～31日的FY-4A軌道進行定軌，FY-4A每日17～18時進行變軌，本文僅分析變軌前后的軌道。統計多天的定軌結果，高度角殘差rms為0.25＂，方位角殘差rms為0.45＂，具體定軌殘差見表3。

表3 高度角與方位角擬合殘差

將有測角數據的定軌弧段和參考軌道重疊比較位置差異，參考軌道基于FY-4A測距數據精密定軌獲得，測距數據由北京、佳木斯、烏魯木齊、廣州以及騰沖的測站提供，采樣率為一秒一個點。在有測角數據的弧段，軌道精度總體優于50 m。具體定軌精度見表4。

表4 存在測角數據弧段的軌道精度

在定軌弧段內，包含測角數據弧段的定軌精度優于其它弧段，結果見圖4。

T1、T2-Period with angel-data圖4 12月29～30日測角數據定軌精度Fig.4 Precision of orbit determination with optical angle measurement data from 29th December to 30th December

以本次動量輪卸載后到下次動量輪卸載前的弧段為一個定軌周期，在一個定軌周期內，隨著測角數據的增加，定軌結果與參考軌道的位置差異會呈現下降的趨勢，具體結果見圖5。

圖5 3個單天的測角數據定軌精度Fig.5 Precision of orbit determination with optical angle measurement data of three single day

以FY-4A運控軌道為參考軌道，選擇12月29日和12月30日有測角數據的弧段，計算測角殘差O-C，結果見圖6、圖7。高度角與方位角殘差精度在絕大多數的弧段優于0.5＂，高度角平均O-C為0.23＂，方位角平均O-C為0.34＂，具體結果見圖6和圖7。

圖6 12月29日測角殘差O-CFig.6 Residual O-C of orbit determination for FY-4A on 29th December

圖7 12月30日測角殘差O-CFig.7 Residual O-C of orbit determination for FY-4A on 30th December

基于測距數據的軌道精度約為30 m，本次實驗取得的測角殘差O-C約為0.2＂，故無法利用測角殘差評估本次定軌結果。

4.2 測角和測距聯合數據定軌及精度分析

基于測距數據和50 cm望遠鏡的測角數據對2019年12月29日和12月30日的FY-4A進行聯合定軌。利用12月28日18時-12月29日6時的測距和測角數據(12 h定軌弧段)，以及12月29日5時-12月29日17時的測距和測角數據分別定軌(12 h定軌弧段)，比較12月29日5時-6時的重疊軌道，軌道精度約為15 m。利用12月29日18時-12月30日6時的測距和測角數據(12 h定軌弧段)，以及12月30日5時-12月30日17時的測距和測角數據分別定軌(12 h定軌弧段)，比較12月30日5時-6時的重疊軌道，軌道精度約為13 m，具體定軌結果見圖8和圖9。

圖8 12月29日1 h軌道重疊精度Fig.8 Precision of orbit determination of 1 h overlap on 29th December

圖9 12月30日1 h軌道重疊精度Fig.9 Precision of orbit determination of 1 h overlap on 30th December

以測距+測角數據定軌結果為參考軌道，對12月29日和12月30日運控測距軌道進行精度評估，比較位置差異。評估結果表明，在每日衛星機動(機動時間約為UTC17:30)之后，控后半小時定軌預報半小時軌道精度約為30～45 m，此后隨著測距數據的增加，定軌精度穩步提升，控后6 h至下次機動，軌道精度穩定在20 m以內，具體結果見圖10和圖11。

圖10 12月29日測距和測角聯合定軌精度Fig.10 Precision of orbit determination with optical angel measurement data and ranging data for FY-4A on 29th December

圖11 12月30日測距和測角聯合定軌精度Fig.11 Precision of orbit determination with optical angel measurement data and ranging data for FY-4A on 30th December

5 結論

基于上海天文臺GEO全球聯合觀測網50 cm望遠鏡的光學測角數據對FY-4A進行定軌分析，剔除動量輪卸載變軌的弧段，單天測角數據定軌的高度角殘差rms優于0.25＂，方位角殘差rms優于0.45＂。與基于測距數據提供的運控軌道相比，在有測角數據的弧度，軌道精度優于50 m。在單天定軌弧段內，隨著測角數據量的增加，軌道的精度呈現上升趨勢。利用測角和測距數據對FY-4A進行聯合定軌，以多組12 h的數據為定軌弧段，兩組數據存在1 h的相同弧段，重疊弧段的精度優于15 m。以一個單天的數據為定軌弧段，高度角平均O-C為0.23＂，方位角平均O-C為0.34＂。在FY-4A動量輪卸載后，控制后半小時軌道精度在30-45 m，隨著測量數據的增加，定軌精度穩定提升，控后6 h至下次變軌前，軌道精度優于20 m。