空間目標的雷達定軌坐標轉換問題

黃曉斌， 張 燕， 肖 銳

（空軍預警學院，武漢 430019）

0 引言

空間目標軌道確定是空間目標監視雷達的一項重要任務，主要的數據處理步驟（如圖 1 所示）是：①坐標轉換；②初軌確定；③軌道改進；④編目。在數據處理過程中，涉及到大量的天文學和數學知識。雖然國內外已有一些開源定軌軟件可供借鑒，如 OrbFit[1]、ORSA[2]、NEOPROP[3]、Novas[4]、ODTBX[5]等，但這些軟件都是針對天體軌道問題而研發的，且主要是處理光學觀測數據，并不適合處理雷達觀測數據。為此，本文開發了 1 套雷達定軌函數庫（RadarOrbDet），能夠有助于雷達技戰指標的快速實現。RadarOrbDet 庫按照上述4 個處理步驟開發，本文主要介紹其坐標轉換模塊。

本文通過討論雷達定軌中涉及的基本坐標系、坐標轉換的數學原理，給出程序設計思路，最后利用衛星工具包（satellite tool kit,STK）軟件驗證開發模塊的有效性。

圖1 雷達定軌數據處理流程

1 坐標系簡介

雷達觀測基于地球坐標系，空間目標軌道是基于天球坐標系，這就涉及到地球坐標系與天球坐標系之間的轉換。

協議天球坐標系由國際天文學聯合會（The International Astronomical Union, IAU）和國際地球自轉和參考系服務組織（The International Earth Rotation and Reference Systems Service, IERS）發布，目前采用的是國際天球參考系（international celestial reference frame, ICRS）。依據坐標原點的不同，ICRS 可分為太陽系質心天球參考系（barycentric celestial reference system, BCRS）和地球質心天球參考系（geocentric celestial reference system, GCRS）。BCRS 用于計算行星的運動軌道及編制星表；GCRS用于計算衛星軌道及編制衛星星歷。ICRS 由國際天球參考框架（international celestial reference frame,ICRF）來實現。在1997 年IAU 第23 屆大會上，通過并決定自1998-01-01 起，在天文研究、空間探測、大地測量以及地球動力學等領域中采用ICRS[6]。

協議地球坐標系由國際地球參考系（international terrestrial reference system, ITRS）實現。GCRS 是1 個相當好的準慣性系，衛星的軌道計算一般都是在 GCRS 中進行。這就必須涉及到GCRS 與ITRS 間的坐標轉換問題[7]。

2 坐標轉換原理

在雷達定軌中，需要將在站心地平坐標系下的雷達觀測數據轉換到GCRS 坐標系中。

2.1 站心地平坐標系Xh與ITRS 坐標系XGO的轉換

站心地平坐標系Xh與 ITRS 坐標系XGO的定義如表1 所示。

表1 站心地平坐標系與ITRS 坐標系的定義

站心地平坐標系與ITRS 坐標系之間的轉換公式為

2.2 ITRS 與GCRS 坐標系的轉換

ITRS 與GCRS 的轉換早期是基于春分點的。目前，IERS 2010 年建議使用基于無旋轉原點（nonrotating origin, NRO）的轉換方法[9]。基于 IAU 2006/2000A-CIO 模型的轉換流程[9]如圖2 所示。

圖2 “IAU 2006/2000A-CIO based”坐標轉換流程

轉換過程中涉及到2 個中間坐標系：地球中間坐標系（terrestrial intermediate reference system, TIRS）和天球中間坐標系（celestial intermediate reference system, TIRS）它們的定義見表2。

表2 TIRS 坐標系與CIRS 坐標系的定義

表 2 中：TIP（terrestrial intermediate pole）為地球瞬時自轉軸；TIO（terrestrial intermediate origin）為地球中間零點；CIP（celestial intermediate pole）為天球瞬時自轉軸；CIO（celestial intermediate origin）為天球中間零點。

在t時刻，ITRS 和 GCRS 的轉換是 2 個 3 維直角坐標系間的轉換，可以寫成

式中：M(t)、RCIO(t)和W(t)分別是由于 CIP 在 GCRS中的運動（歲差章動）、地球的自轉以及 CIP 在ITRS 中的運動（極移）引起的旋轉矩陣。它們的具體表達式見文獻[9]。

3 程序設計

圖3 給出了坐標轉換的程序設計流程圖，其轉換過程分為 3 個部分：①將目標在測站地平坐標系的極坐標轉換為ITRS 坐標系下的直角坐標；②根據轉換時刻以及地球定向參數（Earth orientation parameters，EOP）數據形成 ITRS 到 GCRS 的轉換矩陣；③將轉換矩陣與ITRS 下的直角坐標相乘得到GCRS 下的直角坐標。轉換流程中涉及一些時間系統的轉換可參見文獻[7]。圖3 中的“TT 時間”表示地球時（Terrestrial Time）。

圖3 坐標轉換程序設計流程

在轉換過程中有如下幾個問題應該注意：

1）按照式（1）的要求，應該輸入站點的天文經緯度和地理經緯度。前者用來進行坐標旋轉，后者用來計算站點的大地直角坐標。但通常只給出站點的地理經緯度，在進行坐標旋轉時，用地理經緯度替代天文經緯度，2 者雖有細微的差別，但在絕大多數場合可忽略不計，如果需要極高坐標轉換精度時應加以區分。

2）轉換過程中需要地球定向參數（Earth orientation parameters, EOP）數據，它可從網站下載獲得[10]，要注意表中數據的有效時段和數據條目的時間間隔（通常為1 d）。因此，在給定轉換時刻時，通常需要依據該表插值來獲取更為精確的數據，而當轉換時刻超出EOP 數據表有效時間范圍時，一般可取最接近該轉換時刻的數據來替代，但當超出時間范圍間隔過大時，其誤差應加以考慮。

3）當需要極高坐標轉換精度時，在用EOP 數據表插值計算極移和世界協調時（coordinated universal time，UTC）與 1 類世界時（universal time，UT1）參數時，需要額外考慮海潮和固體潮的修正[7]，但同時會導致轉換速度的下降。

4 軟件接口調用及有效性驗證

4.1 接口調用示例

RadarOrbDet 庫坐標轉換功能的使用方法如下：

步驟①：打開VS2015，新建控制臺應用程序，并命名為 radarorddet_Test，輸出模式設為 64 位Debug，具體方法可參考文獻[11]。

步驟②：從百度網盤下載RadarOrbDet 開發包[12]，并拷貝到工程目錄中，如圖4 所示。

圖4 測試工程目錄結構

步驟③：在radarorbdet_Test.cpp 文件中輸入如圖5 所示的代碼。

圖5 函數接口調用示例

4.2 接口函數有效性驗證

在 STK 11.2 版中仿真 1 顆衛星，其軌道參數設置如圖 6 所示，測站位置參數設置如圖 7所示。

圖6 衛星軌道參數設置

圖7 測站位置參數設置

仿真時間段設在2016-09-15 UTCG（格林尼治協調世界時） 04:00:00 至 2016-09-16 UTCG 04:00:00（STK 沒有最新的EOP 數據，故設此時間段保證STK 自帶的EOP 數據有效）。利用STK 的報表功能輸出測站對衛星的觀測時刻、觀測距離、方位和俯仰數據，同時輸出STK 內部計算所得的衛星GCRS 坐標，坐標數據的采樣時間段為2016-09-15 UTCG04:00:00 至 2016-09-15UTCG 04:02:00，間隔1 s，總共121 點數據，圖8(a)顯示了距離測量數據，圖8(b)顯示了方位測量數據，圖8(c)顯示了俯仰測量數據。

利用 RadarOrbDet 軟件包計算 GCRS 坐標（XROD,YROD,ZROD），采用式（3）描述的 GCRS 坐標相對誤差與STK 的輸出進行比較，圖9 顯示了相對誤差曲線（實線）。從圖9 中可以看出相對誤差在 1×10-7水平，經分析表明，誤差的主要原因是STK 報表以文本方式輸出數據，導致數據產生精度截斷誤差，由此可知軟件包算法是有效的。此外，在圖9 中還繪制了俯仰角曲線（虛線），可以看出當俯仰角越小時，轉換的誤差越大，其原因是俯仰角越小時，截斷誤差對坐標轉換誤差的影響越大。

圖8 測站觀測的距離、方位和俯仰數據

圖9 GCRS 坐標轉換相對誤差

5 結束語

本文介紹了 RadarOrbDet 庫中坐標轉換模塊的基本數學原理，給出了程序設計思路，描述了調用方法，最后用STK 驗證了算法的有效性。下一步將繼續介紹初軌確定模塊的功能。