張子林

（1.惠州大亞灣經(jīng)濟(jì)技術(shù)開發(fā)區(qū)國土資源測繪所，廣東 惠州 516081）

國際陸地參考系（ITRS）是國際應(yīng)用最廣泛、影響最深遠(yuǎn)的坐標(biāo)系統(tǒng)，其具體實(shí)現(xiàn)的坐標(biāo)框架稱為國際陸地參考框架（ITRF）。很多國家本土獨(dú)立坐標(biāo)系都以ITRF為標(biāo)準(zhǔn)建立，比如我們國家的CGCS2000坐標(biāo)[1-3]，對應(yīng)于ITRF97框架的2000.0歷元。由于地球內(nèi)部（地殼）不斷運(yùn)動(dòng)，造成地表站點(diǎn)的位移發(fā)生變化，坐標(biāo)系統(tǒng)也就需要不斷地維護(hù)更新核心站點(diǎn)的坐標(biāo)，每隔一段時(shí)間（一般以年為單位）就需發(fā)布新的坐標(biāo)框架。至今，ITRF已經(jīng)發(fā)布了超10個(gè)版本，最新的框架版本是ITRF2014。

不同ITRF框架版本間有嚴(yán)密的數(shù)學(xué)模型轉(zhuǎn)換[4]， 然而，由于對地殼具體運(yùn)動(dòng)程度不明確，框架內(nèi)不同歷元轉(zhuǎn)換困難。有不少學(xué)者進(jìn)行了深入研究，魏子卿院士提出4種方法解決歷元轉(zhuǎn)換難題，從精度、實(shí)用角度做了詳細(xì)闡述；在Trimble開發(fā)的CenterPoint RTX后處理服務(wù)，也明確提供了全球地殼板塊的改正模型；GAMIT/GLOBK軟件在tables中提供了ITRF核心的速度場模型，能夠?qū)崿F(xiàn)不同歷元 轉(zhuǎn)換。

本文結(jié)合多年研究與工作經(jīng)驗(yàn)，總結(jié)、梳理前人的研究成果，提出一整套詳盡、可操作的實(shí)用方法，實(shí)現(xiàn)我國大陸地區(qū)ITRF任意框架、任意歷元的坐標(biāo)轉(zhuǎn)換。同時(shí)，提出快速獲取CGCS2000坐標(biāo)的新方法。

1 坐標(biāo)框架轉(zhuǎn)換模型

同一歷元下，不同的ITRF框架版本之間存在嚴(yán)密的數(shù)學(xué)模型轉(zhuǎn)換，遵循空間直角坐標(biāo)系的布爾莎七參數(shù)模型，在歷元t時(shí)刻，轉(zhuǎn)公式如下[5-6]：

式中，X,Y,Z表示原ITRF框架下的坐標(biāo)；XS,YS,ZS表示目標(biāo)框架下坐標(biāo)；Tx、Ty、Tz表示兩個(gè)框架版本坐標(biāo)原點(diǎn)的三分量偏差值，單位是mm；D表示尺度縮放因子，單位是PPM（百萬分之一）；Rx、Ry、Rz表示旋轉(zhuǎn)角度因子，單位為mas。具體的數(shù)學(xué)模型各參數(shù)由ITRF官網(wǎng)（https://itrf.ign.fr/trans_para.php）提供。官方提供的參數(shù)對應(yīng)有參考時(shí)刻和參數(shù)變化率，即15個(gè)參數(shù)（7參數(shù)+7參數(shù)變化率+1參考時(shí)刻）。表1展示ITRF2014轉(zhuǎn)ITRF2008框架的參數(shù)。

表1 ITRF2014轉(zhuǎn)換到ITRF2008的參數(shù)

根據(jù)官方提供的每個(gè)轉(zhuǎn)換參數(shù)、參考?xì)v元和變化率，在任意歷元t時(shí)刻下，轉(zhuǎn)換參數(shù)求解公式如下：

式中，EPOCH為表1中的參考時(shí)刻2 010.0；P、P.分別表示t時(shí)刻的轉(zhuǎn)換參數(shù)和參數(shù)速率。

測試我國大陸不同地殼活動(dòng)塊體區(qū)域的4個(gè)IGS核心站（圖1）在ITRF2014，2020歷元下坐標(biāo)值轉(zhuǎn)換到ITRF2008，2020歷元下（表2）。UBUM站位于西域活動(dòng)地塊，LHAZ站位于青藏活動(dòng)地塊，WUHN站位于華南活動(dòng)地塊，SHAO站位于華北活動(dòng)地塊。根據(jù)坐標(biāo)框架轉(zhuǎn)換模型計(jì)算結(jié)果發(fā)現(xiàn)相同歷元，不同活 動(dòng)塊體的站點(diǎn)框架坐標(biāo)偏差均很?。? mm以內(nèi)）。

表2 ITRF2008與ITRF2014框架轉(zhuǎn)換

圖1 中國大陸活動(dòng)地塊及4個(gè)IGS核心站點(diǎn)分布圖

2 坐標(biāo)歷元轉(zhuǎn)換模型

在同一坐標(biāo)框架下，不同歷元的坐標(biāo)轉(zhuǎn)換公式（3），其中Xt1、Yt1、Zt1和Xt2、Yt2、Zt2分別是同框架下歷元t1和t2的坐標(biāo)。根據(jù)公式可知，坐標(biāo)歷元轉(zhuǎn)換最核心的問題是獲取準(zhǔn)確的站點(diǎn)速度值（或者塊體運(yùn)動(dòng)模型）Vx、Vy、Vz。

以往獲得速度值（場）模型的研究有很多[7-9]，例如全域歐拉矢量法、局域歐拉矢量法、格網(wǎng)平均值法和塊體歐拉矢量法等。魏子卿研究結(jié)果證明，局域歐拉矢量法獲得的速度精度最高，全域歐拉矢量法速度精度最低，格網(wǎng)平均值法使用方便，且獲得的速度精度較高，適合一般用戶使用。

圖2展示我國大陸地區(qū)137個(gè)3 ×3 的格網(wǎng)[7]，本文收集公開發(fā)表的1998－2014年間全國GNSS站點(diǎn)地表運(yùn)動(dòng)背景速度場計(jì)算出每個(gè)格網(wǎng)在CGCS2000坐標(biāo)框架下的平均速度值，此值代表格網(wǎng)內(nèi)所有點(diǎn)的速度。

圖2 中國大陸地區(qū)3 ×3 地殼運(yùn)動(dòng)速度格網(wǎng)

3 任意框架和歷元轉(zhuǎn)換實(shí)驗(yàn)

為了驗(yàn)證本文歷元轉(zhuǎn)換模型正確性和有效性，設(shè)計(jì)實(shí)驗(yàn)一，相同框架不同歷元轉(zhuǎn)換。將圖2中的 4個(gè)IGS站在ITRF2008框架歷元2020.0轉(zhuǎn)換到歷元2016.0，結(jié)果顯示表3，對比三分量坐標(biāo)偏差，可以發(fā)現(xiàn)精度均在1 cm左右。

表3 ITRF2008框架下2020.0歷元轉(zhuǎn)換至2016.0歷元

實(shí)驗(yàn)二：不同框架和歷元轉(zhuǎn)換；針對不同坐標(biāo)系統(tǒng)之間的轉(zhuǎn)換，操作流程通過ITRF97框架過渡，即：A坐標(biāo)框架轉(zhuǎn)換到ITRF97下，然后利用圖2中的格網(wǎng)速度模型進(jìn)行歷元轉(zhuǎn)換，再經(jīng)過框架轉(zhuǎn)換到B坐標(biāo)。

筆者在日常工作中，經(jīng)常遇到工程項(xiàng)目中需獲取測區(qū)CGCS2000坐標(biāo)[9-10]；通過GNSS靜態(tài)觀測以及GAMIT/GLOBK軟件解算獲得毫米級高精度的ITRF2008框架2020.0歷元坐標(biāo)，通過框架轉(zhuǎn)換到ITRF97框架2000.0歷元，即獲得CGCS2000坐標(biāo)系統(tǒng)。

為了驗(yàn)證方法可行，筆者收集了部分國家B級控制點(diǎn)，按照本文介紹的方法進(jìn)行坐標(biāo)轉(zhuǎn)換計(jì)算，表4統(tǒng)計(jì)轉(zhuǎn)換結(jié)果，對比B等級控制點(diǎn)真實(shí)坐標(biāo)，發(fā)現(xiàn)本文坐標(biāo)系統(tǒng)轉(zhuǎn)換精度較高，三分量較差可以達(dá)到1cm左右精度。

表4 ITRF2008坐標(biāo)轉(zhuǎn)換CGCS2000坐標(biāo)

4 結(jié) 語

本文研究了坐標(biāo)系統(tǒng)轉(zhuǎn)換，涉及框架轉(zhuǎn)換和歷元轉(zhuǎn)換，提出了一種實(shí)際可行的歷元轉(zhuǎn)換方法。通過總結(jié)以往研究中的各類方法，提供我國大陸地區(qū)3 ×3 地殼運(yùn)動(dòng)速度格網(wǎng)，能夠在工程應(yīng)用中快速獲得CGCS2000坐標(biāo)。

本文研究結(jié)果表明，相同歷元下，不同框架之間的坐標(biāo)差異較小，不同活動(dòng)塊體下ITRF2008和ITRF2014在2020.0歷元時(shí)各分量坐標(biāo)差小于3 mm；受不同區(qū)域地殼活動(dòng)程度影響，相同框架下，不同歷元之間的坐標(biāo)差異明顯。本文提供的歷元轉(zhuǎn)換方法，實(shí)際測試精度可達(dá)1 cm左右，能夠滿足工程應(yīng)用。