劉思語(yǔ)，黃勇，毛銀盾，賈耀紅，魯文強(qiáng),4，黃乘利,3，鄭景輝，楊鵬

1. 中國(guó)科學(xué)院上海天文臺(tái)，上海 200030 2. 中國(guó)科學(xué)院大學(xué) 天文與空間科學(xué)學(xué)院，北京 100049 3. 上?？萍即髮W(xué) 物質(zhì)學(xué)院，上海 201210 4. 國(guó)家衛(wèi)星氣象中心，北京 100081

1 引言

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

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

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

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

本文利用光學(xué)測(cè)角數(shù)據(jù)對(duì)FY-4A數(shù)據(jù)進(jìn)行定軌分析，分析了僅利用測(cè)角數(shù)據(jù)定軌以及利用測(cè)距測(cè)角數(shù)據(jù)聯(lián)合定軌的精度，并評(píng)估了僅利用測(cè)距數(shù)據(jù)的實(shí)時(shí)軌道精度。

2 觀測(cè)實(shí)驗(yàn)

2019年12月27～31日，上海天文臺(tái)GEO全球聯(lián)合觀測(cè)網(wǎng)50 cm望遠(yuǎn)鏡(以下稱50 cm望遠(yuǎn)鏡)對(duì)FY-4A進(jìn)行了光學(xué)觀測(cè)，觀測(cè)樣例見圖1和圖2。利用CCD相機(jī)觀測(cè)得到FY-4A的測(cè)角數(shù)據(jù)，直接測(cè)量量為目標(biāo)星的量度坐標(biāo)，再轉(zhuǎn)換為定軌所需要的赤經(jīng)赤緯或高度角與方位角。

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

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

圖2 利用凝視模式觀測(cè)GEO衛(wèi)星實(shí)例Fig.2 An example of stare-mode for GEO satellites

表1 望遠(yuǎn)鏡參數(shù)及測(cè)站位置

表2 光學(xué)測(cè)角弧段

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

3 測(cè)量模型

3.1 照片天體測(cè)量與底片模型

圖3 赤經(jīng)赤緯與理想坐標(biāo)幾何關(guān)系Fig.3 The geometric relationship between ideal coordinate and ascension,declination

照相天體測(cè)量技術(shù)基于底片上目標(biāo)星的量度坐標(biāo)推算其理想坐標(biāo)，并計(jì)算其赤經(jīng)赤緯，赤經(jīng)赤緯與理想坐標(biāo)的幾何關(guān)系見圖3。根據(jù)觀測(cè)結(jié)果所計(jì)算的赤經(jīng)和赤緯均在站心赤道坐標(biāo)系下表示。

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

(1)

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

(2)

聯(lián)立(1)(2)有：

(3)

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

(4)

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

(5)

本文定軌所使用觀測(cè)資料為高度角與方位角，赤經(jīng)赤緯與高度角方位角存在以下轉(zhuǎn)換關(guān)系，A和E為目標(biāo)星的方位角和高度角，S為測(cè)站所在恒星時(shí)。

3.2 高度角與方位角

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

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

方位角A對(duì)衛(wèi)星狀態(tài)的偏導(dǎo)數(shù)：

(6)

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

(7)

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

高度角E表現(xiàn)為目標(biāo)位置矢量與Zs軸的夾角：

高度角E對(duì)衛(wèi)星狀態(tài)參數(shù)的偏導(dǎo)數(shù)：

3.3 測(cè)距

測(cè)距型數(shù)據(jù)獲取的是測(cè)站與衛(wèi)星相對(duì)距離關(guān)系，可表示為：

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

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

3.4 定軌原理

在對(duì)FY-4A進(jìn)行光學(xué)觀測(cè)的過程中。假設(shè)在ti時(shí)刻獲得一組觀測(cè)量Yi:

(8)

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

4 數(shù)據(jù)分析及小結(jié)

4.1 測(cè)角數(shù)據(jù)定軌及精度分析

基于50 cm望遠(yuǎn)鏡的測(cè)角數(shù)據(jù)，對(duì)2019年12月28～31日的FY-4A軌道進(jìn)行定軌，F(xiàn)Y-4A每日17～18時(shí)進(jìn)行變軌，本文僅分析變軌前后的軌道。統(tǒng)計(jì)多天的定軌結(jié)果，高度角殘差rms為0.25＂，方位角殘差rms為0.45＂，具體定軌殘差見表3。

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

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

表4 存在測(cè)角數(shù)據(jù)弧段的軌道精度

在定軌弧段內(nèi)，包含測(cè)角數(shù)據(jù)弧段的定軌精度優(yōu)于其它弧段，結(jié)果見圖4。

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

以本次動(dòng)量輪卸載后到下次動(dòng)量輪卸載前的弧段為一個(gè)定軌周期，在一個(gè)定軌周期內(nèi)，隨著測(cè)角數(shù)據(jù)的增加，定軌結(jié)果與參考軌道的位置差異會(huì)呈現(xiàn)下降的趨勢(shì)，具體結(jié)果見圖5。

圖5 3個(gè)單天的測(cè)角數(shù)據(jù)定軌精度Fig.5 Precision of orbit determination with optical angle measurement data of three single day

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

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

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

基于測(cè)距數(shù)據(jù)的軌道精度約為30 m，本次實(shí)驗(yàn)取得的測(cè)角殘差O-C約為0.2＂，故無法利用測(cè)角殘差評(píng)估本次定軌結(jié)果。

4.2 測(cè)角和測(cè)距聯(lián)合數(shù)據(jù)定軌及精度分析

基于測(cè)距數(shù)據(jù)和50 cm望遠(yuǎn)鏡的測(cè)角數(shù)據(jù)對(duì)2019年12月29日和12月30日的FY-4A進(jìn)行聯(lián)合定軌。利用12月28日18時(shí)-12月29日6時(shí)的測(cè)距和測(cè)角數(shù)據(jù)(12 h定軌弧段)，以及12月29日5時(shí)-12月29日17時(shí)的測(cè)距和測(cè)角數(shù)據(jù)分別定軌(12 h定軌弧段)，比較12月29日5時(shí)-6時(shí)的重疊軌道，軌道精度約為15 m。利用12月29日18時(shí)-12月30日6時(shí)的測(cè)距和測(cè)角數(shù)據(jù)(12 h定軌弧段)，以及12月30日5時(shí)-12月30日17時(shí)的測(cè)距和測(cè)角數(shù)據(jù)分別定軌(12 h定軌弧段)，比較12月30日5時(shí)-6時(shí)的重疊軌道，軌道精度約為13 m，具體定軌結(jié)果見圖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

以測(cè)距+測(cè)角數(shù)據(jù)定軌結(jié)果為參考軌道，對(duì)12月29日和12月30日運(yùn)控測(cè)距軌道進(jìn)行精度評(píng)估，比較位置差異。評(píng)估結(jié)果表明，在每日衛(wèi)星機(jī)動(dòng)(機(jī)動(dòng)時(shí)間約為UTC17:30)之后，控后半小時(shí)定軌預(yù)報(bào)半小時(shí)軌道精度約為30～45 m，此后隨著測(cè)距數(shù)據(jù)的增加，定軌精度穩(wěn)步提升，控后6 h至下次機(jī)動(dòng)，軌道精度穩(wěn)定在20 m以內(nèi)，具體結(jié)果見圖10和圖11。

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

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

5 結(jié)論

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