杜仲進(jìn)

（福建省測繪院，福州 350003）

0 引言

由于地殼運(yùn)動(dòng)和形變等因素引起的站點(diǎn)坐標(biāo)變化對(duì)現(xiàn)代毫米級(jí)大地測量技術(shù)具有很大的影響，測站高精度坐標(biāo)必須通過速度場歸算到用戶所需的歷元時(shí)刻才具有工程實(shí)用價(jià)值。速度場的擬合目前除了使用NNR-NUVEL1A[1]模型，主流方案是利用近幾十年來采集的中國大陸空間大地測量數(shù) 據(jù)求取速度場。

自2011 年12 月中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)即陸態(tài)網(wǎng)（crustal movement observation network of China， CMONOC）完成全部建設(shè)任務(wù)后，CMONOC便成為推算中國大陸或局部區(qū)域速度場的主要數(shù)據(jù)源[2-5]。常用的速度場模型建立方法主要分為歐拉矢量法和數(shù)值擬合法：歐拉矢量法分為整體、塊體和局域歐拉矢量法[6-8]，其中局部歐拉矢量法精度最高，可達(dá)到1～2 mm/a 的精度[9]。數(shù)值擬合方法有多種，常用的有臨近點(diǎn)插值、反距離插值、克里金插值、最小曲率法、三角網(wǎng)插值和多面函數(shù)插值等方法，其中臨近點(diǎn)和反距離插值法原理最為簡單，但僅考慮距離作為權(quán)重標(biāo)準(zhǔn)不太合理[10-11]；最小曲率法是構(gòu)造1 個(gè)具有最小曲率的非精確曲面，需要經(jīng)過迭代處理[12-13]；克里金法考慮了隨機(jī)數(shù)據(jù)的期望和方差等信息，適用性較廣且精度通常較高，但原理復(fù)雜，特別是半變異函數(shù)的選取直接影響到結(jié)果的精度[14]；多面函數(shù)法是采用一系列單值數(shù)學(xué)面疊加逼近某一曲面，精度較高，但核函數(shù)和光滑系數(shù)的選擇還是要具體問題具體分析[15-16]；三角網(wǎng)法原理簡單，在考慮了速度和方位信息的同時(shí)，可以刻畫更為精細(xì)的局部特征，因此插值精度高，適用性廣[12]；另外，采用數(shù)值擬合方法得到格網(wǎng)速度場，使用更為簡便。

考慮到沿海區(qū)域速度場模型的數(shù)據(jù)現(xiàn)勢性、高精度性和簡單適用性原則，本文基于福建省及周邊區(qū)域測站2011—2017 年近7 a 的觀測數(shù)據(jù)進(jìn)行處理分析，利用局部優(yōu)化過程（local optimization procedure， LOP）三角網(wǎng)反距離加權(quán)格網(wǎng)模型構(gòu)建福建省區(qū)域格網(wǎng)模型，并討論了區(qū)域歐拉矢量模型的應(yīng)用。

1 基于GAMIT/GLOBK 的福建CORS 數(shù)據(jù)解算

CMONOC 主要由均勻分布在全國的260 余個(gè)連續(xù)運(yùn)行參考站（continuously operating reference stations， CORS）和2 000 余個(gè)不定期觀測區(qū)域站構(gòu)成。為了研究福建省及周邊地區(qū)水平速度場模型，本文選取CMONOC 中位于福建省及其周邊地區(qū)的7 個(gè)CORS 和68 個(gè)區(qū)域站進(jìn)行速度場解算和分析。其中，CORS 觀測數(shù)據(jù)收集大多起始于2010—2011 年，總體時(shí)間跨度約7～8 a，區(qū)域站觀測數(shù)據(jù)為2009 年以后的數(shù)據(jù)，復(fù)測頻率定為 每2 年1 次。研究區(qū)域的站點(diǎn)分布如圖1 所示，圓形代表區(qū)域站，三角形代表CORS。

本文采用 GAMIT/GLOBK 軟件處理福建省及周邊區(qū)域的站點(diǎn)觀測數(shù)據(jù)，為了將同全球衛(wèi)星導(dǎo) 航 系 統(tǒng)(global navigation satellite system， GNSS）測站同國際GNSS 服務(wù)組織(International GNSS Service， IGS）站進(jìn)行聯(lián)合解算，還需選取國內(nèi)及周邊的17 個(gè)IGS 站（AIRA、BJFS、DAEJ、IISC、IRKT、KIT3、KUNM、LHAZ、PIMO、POL2、SHAO、URUM、WUHN、YSSK、TIXI、ARTU和TCMS）進(jìn)行約束。受計(jì)算機(jī)軟件和硬件性能的限制，采用GAMIT 軟件進(jìn)行基線解算時(shí)，一般要求測站數(shù)量不超過60 個(gè)，否則會(huì)大大降低基線解算效率，因此首先將福建省及周邊的測站按照區(qū)域劃分原則分為2 個(gè)分區(qū)，加上17 個(gè)IGS約束站，每個(gè)分區(qū)測站數(shù)不超過55 個(gè)；設(shè)置好基線解算參數(shù)進(jìn)行解算，得到每個(gè)分區(qū)的單日松弛解文件，然后采用GLRED 軟件將全球IGS 單日解與各個(gè)分區(qū)的單日解文件進(jìn)行合并，合并時(shí)僅保留穩(wěn)定的IGS 核心站，需剔除出現(xiàn)粗差的測站信息，合并后得到包含全球IGS 核心站和福建省及周邊區(qū)域測站的坐標(biāo)、衛(wèi)星軌道和極移等參數(shù)的單日整體松弛解；最后將所有的單日整體松弛解法方程文件進(jìn)行網(wǎng)平差，解算福建省及周邊區(qū)域測站在國際地球參考框架2014（international terrestrial reference frame， ITRF2014）下的坐標(biāo)和速度。繪制研究區(qū)域測站的速度場，如圖 2所示。

圖2 福建省及周邊區(qū)域速度場矢量圖

2 福建省速度場模型的構(gòu)建原理

2.1 局部地區(qū)歐拉矢量模型

由于福建省所占區(qū)域不大，且地表內(nèi)部構(gòu)造活動(dòng)不劇烈，因此可以把整個(gè)福建及周邊區(qū)域地表當(dāng)作1 個(gè)剛性塊體，根據(jù)歐拉定理[13]可得

式中：V為塊體內(nèi)某站點(diǎn)的速度；Ω為該站點(diǎn)所處塊體的歐拉矢量；R為站點(diǎn)的位置矢量。式（1）改寫成便于計(jì)算的公式[13]為

式中：VE、VN為某測站的站心坐標(biāo)速度，假設(shè)塊體運(yùn)動(dòng)沿水平方向與VU無關(guān)；r 為地球平均半徑，約為6 371 393 m；λ、φ為該測站的經(jīng)緯度，單位為弧度；Ωx、Ωy、Ωz為歐拉矢量的3 分量。

設(shè)歐拉矢量的估值為參數(shù)向量 ?X ，站點(diǎn)的水平速度觀測向量L，觀測誤差向量為Δ，誤差矩陣為V，則可以組成誤差方程為

上述方程的最小二乘解為

以福建省及周邊區(qū)域測站的坐標(biāo)和速度代入上述方程進(jìn)行平差解算，為避免粗差對(duì)參數(shù)解算的影響，需進(jìn)行迭代最小二乘計(jì)算，每次解算后以3 倍中誤差準(zhǔn)則判斷數(shù)據(jù)集中是否含有粗差，直至不再含有粗差時(shí)結(jié)束迭代，得到研究區(qū)域的歐拉矢量參數(shù)。

2.2 Delaunay 三角網(wǎng)反距離加權(quán)法構(gòu)建區(qū)域格網(wǎng)速度場

2.2.1 格網(wǎng)點(diǎn)生成

以115～121°E、23～28°N 為矩形邊界框，以0.5°為格網(wǎng)間距進(jìn)行格網(wǎng)劃分，得到福建省及周邊區(qū)域的經(jīng)緯度格網(wǎng)點(diǎn)。

2.2.2 Delaunay 三角網(wǎng)生成

進(jìn)行德洛奈（Delaunay）三角剖分，必須要滿足2 個(gè)重要準(zhǔn)則：唯一性和最大化最小角特性。根據(jù)Delaunay 三角網(wǎng)的實(shí)現(xiàn)過程的不同，通常分為逐點(diǎn)插入法、三角網(wǎng)生成法、分治算法和合成算法等[17]。其中逐點(diǎn)插入法原理簡單，易于編程實(shí)現(xiàn)，基本步驟為：1）建立包含所有離散點(diǎn)集的多邊形凸包；2）在構(gòu)成凸包的所有點(diǎn)中，初始化凸包三角網(wǎng)；3）根據(jù)局部優(yōu)化過程局部優(yōu)化過程（local optimization procedure， LOP）算法優(yōu)化初始凸包三角網(wǎng)；4）在已生成的三角網(wǎng)逐點(diǎn)插入其余散點(diǎn)，每次插入1 個(gè)點(diǎn)后，用LOP 算法重新生成三角網(wǎng)；5）重復(fù)步驟4），直到所有點(diǎn)插入完畢。

由福建省及周邊區(qū)域的測站生成的Delaunay三角網(wǎng)如圖3 所示。

圖3 福建省及周邊區(qū)域測站三角網(wǎng)

2.2.3 反距離加權(quán)法求取格網(wǎng)點(diǎn)的水平速度值

1）判斷格網(wǎng)點(diǎn)位于某一三角形中。采取射線法判斷1 點(diǎn)是否在區(qū)域內(nèi)，具體步驟為：從該點(diǎn)向左引水平掃描線，遍歷計(jì)算該射線與每1 個(gè)三角網(wǎng)邊界的相較次數(shù)，如果為奇數(shù)，則認(rèn)為在三角形內(nèi)，若為偶數(shù)，則在三角形外部。如果格網(wǎng)點(diǎn)不位于三角網(wǎng)內(nèi)部，則尋找離其最近的3 個(gè)測站點(diǎn)。

2）反距離加權(quán)法計(jì)算格網(wǎng)點(diǎn)水平速度場。設(shè)空 間 待 插 值點(diǎn) 為 P ( xp， yp， zp)，P 點(diǎn) 所在 鄰 域 內(nèi) 有3 個(gè)已知離散點(diǎn) Qi( xi， yi， zi)（ i =1， 2， 3， xi、 yi，為位置信息， zi為屬性信息）。根據(jù)離散點(diǎn)的值，通過反距離加權(quán)法對(duì)P 點(diǎn)的屬性zp進(jìn)行插值，即

3 福建省速度場模型的精度分析

3.1 歐拉矢量模型內(nèi)外符合精度分析

將福建省及周邊區(qū)域當(dāng)作1 個(gè)剛性塊體，使用68 個(gè)區(qū)域站的坐標(biāo)和速度，根據(jù)式（4）用迭代最小二乘法求取該區(qū)域的歐拉矢量，與文獻(xiàn)[6]利用“地殼網(wǎng)絡(luò)工程”1999—2008 年的觀測數(shù)據(jù)解算的速度場模型（CHINA-2008）進(jìn)行比較，具體參數(shù)見表1。根據(jù)解算得到的歐拉矢量估值計(jì)算福建省區(qū)域的歐拉矢量模型速度場殘差值，如圖4 和圖5 所示。

表1 福建省區(qū)域歐拉矢量參數(shù) 單位： rad/Ma

分析表1 可得：2 套參數(shù)的整體旋轉(zhuǎn)趨勢有較大的相似性，但是CHINA-2008 模型代表的是中國整個(gè)大陸地區(qū)的整體歐拉矢量，反映的是整個(gè)大陸運(yùn)動(dòng)趨勢；而福建省歐拉矢量是針對(duì)局部區(qū)域構(gòu)建的，因此更適合表達(dá)福建省自身塊體近年來的運(yùn)動(dòng)趨勢。因此2 套參數(shù)必然會(huì)產(chǎn)生一定的偏差。

圖4 福建省歐拉矢量模型區(qū)域站速度場殘差矢量圖

圖5 福建省歐拉矢量模型區(qū)域站速度場殘差散點(diǎn)圖

分析圖4 及圖5 可得，采用福建省區(qū)域歐拉 矢量模型可以明顯減弱該區(qū)域塊體的運(yùn)動(dòng)趨勢，除個(gè)別測站的速度殘差值較大，出現(xiàn)該現(xiàn)象可能是該測站的實(shí)測速度不準(zhǔn)造成的。經(jīng)歐拉轉(zhuǎn)換后，福建省區(qū)域測站點(diǎn)的速度殘差在2 mm/a 以內(nèi)，大部分在1 mm/a 以內(nèi)?？鄢＝ㄊ^(qū)域歐拉矢量模型后，該區(qū)域仍具有由西向東的微小運(yùn)動(dòng)趨勢。

為了進(jìn)一步驗(yàn)證福建省區(qū)域歐拉矢量模型的精度和可靠性，以該區(qū)域的8 個(gè)CORS 作為外部檢核點(diǎn)，使用解算得到的歐拉矢量參數(shù)計(jì)算檢核點(diǎn)的水平速度，并與原始實(shí)測速度作差得到殘差值，如圖6 所示：外部檢核點(diǎn)殘差絕對(duì)值不超過1 mm/a，其中大部分測站的速度殘差表現(xiàn)為正向偏差，說明解算得到的歐拉矢量參數(shù)存在一定的系統(tǒng)性誤差，該誤差可能是由于部分測站的實(shí)測速度不準(zhǔn)造成的。

圖6 福建省歐拉矢量模型CORS 速度場殘差柱狀圖

統(tǒng)計(jì)福建省歐拉矢量模型的內(nèi)外符合精度指標(biāo)，如表2 所示。

表2 福建省歐拉矢量模型速度場殘差統(tǒng)計(jì) 單位: mm/a

分析表2 可得，東、北方向的內(nèi)外符合速度殘差絕對(duì)值在0.5 mm/a 左右，東（E）、北（N）方向的內(nèi)外符合速度殘差中誤差均不超過1 mm/a，區(qū)域站殘差絕對(duì)值的最大值超過3 mm/a，而CORS殘差絕對(duì)值的最大值不超過1 mm/a，說明區(qū)域站速度場精度低于CORS 速度場精度?？傮w來看，區(qū)域站內(nèi)符合精度與CORS 站外符合精度保持在同一精度水平上，福建省歐拉矢量模型東、北方向均能達(dá)到1 mm/a 的精度。

3.2 Delaunay 三角網(wǎng)反距離加權(quán)格網(wǎng)速度場模型內(nèi)外符合精度分析

在福建省及周邊區(qū)域內(nèi)，使用68 個(gè)區(qū)域站的經(jīng)緯度坐標(biāo)構(gòu)建三角網(wǎng)，并基于反距離加權(quán)原理計(jì)算每個(gè)格網(wǎng)點(diǎn)的水平速度值，如圖7 所示。根據(jù)解算得到的歐拉矢量估值計(jì)算福建省區(qū)域的歐拉矢量模型速度場殘差值，如圖8 所示。

圖7 福建省區(qū)域格網(wǎng)速度場矢量圖

圖8 福建省三角網(wǎng)模型區(qū)域站速度場殘差散點(diǎn)圖

分析圖7、圖8 可得，采用三角網(wǎng)反距離加權(quán)模型得到的格網(wǎng)速度場整體具有西北至東南的運(yùn)動(dòng)方向，整體測站平均運(yùn)動(dòng)速度為：E 方向32.8 mm/a，N 方向-11.7 mm/a。根據(jù)格網(wǎng)速度場反推區(qū)域站的速度值并與實(shí)測值作差得到速度場殘差值，福建省區(qū)域測站點(diǎn)的速度殘差在2 mm/a 以內(nèi)，大部分在1 mm/a 以內(nèi)。同樣，為了進(jìn)一步驗(yàn)證福建省區(qū)域三角網(wǎng)反距離加權(quán)模型的精度和可靠性，以該區(qū)域的8 個(gè)CORS 作為外部檢核點(diǎn)，使用區(qū)域格網(wǎng)速度場計(jì)算檢核點(diǎn)的水平速度，并與原始實(shí)測速度作差得到殘差值，如圖9 所示，外部檢核點(diǎn)殘差絕對(duì)值不超過1 mm/a。

統(tǒng)計(jì)福建省區(qū)域格網(wǎng)模型的內(nèi)外符合精度指標(biāo)，如表3 所示。分析可得與歐拉矢量模型較為一致的結(jié)果，東、北方向的內(nèi)外符合速度殘差絕對(duì)值在0.5 mm/a 左右，東、北方向的內(nèi)外符合速度殘差中誤差均不超過1 mm/a，區(qū)域站殘差絕對(duì)值的最大值超過2.5 mm/a，而CORS 殘差絕對(duì)值的最大值不超過1 mm/a?？傮w來看，區(qū)域站內(nèi)符合精度與CORS外符合精度保持一致，福建省三角網(wǎng)反距離加權(quán)格網(wǎng)模型東、北方向均能達(dá)到1 mm/a 的精度。

圖9 福建省三角網(wǎng)模型CORS 速度場殘差柱狀圖

表3 福建省三角網(wǎng)模型速度場殘差統(tǒng)計(jì) 單位: mm/a

4 結(jié)束語

本文基于福建省及周邊區(qū)域陸態(tài)網(wǎng)測站2011—2017 年長期全球定位系統(tǒng)(global positioning system， GPS）觀測數(shù)據(jù)解算得到福建省區(qū)域高精度水平速度場，分別使用68 個(gè)區(qū)域站和7 個(gè)CORS構(gòu)建和評(píng)估福建省區(qū)域歐拉矢量塊體運(yùn)動(dòng)模型和Delaunay 三角網(wǎng)反距離加權(quán)格網(wǎng)模型，統(tǒng)計(jì)并分析2 種模型的內(nèi)外符合精度指標(biāo)，結(jié)果表明：

1）福建省及周邊區(qū)域整體呈現(xiàn)由西北至東南方向的運(yùn)動(dòng)趨勢，東、北方向整體平均速度分別為32.8、-11.7 mm/a。由于其塊體內(nèi)部構(gòu)造活動(dòng)并不劇烈，因此本文使用的2 種模型的精度相當(dāng)，東、北方向速度場的精度均能達(dá)到1 mm/a 左右。

2）相較于歐拉矢量塊體運(yùn)動(dòng)模型，Delaunay三角網(wǎng)反距離加權(quán)格網(wǎng)模型雖然具體物理意義不強(qiáng)，但是在考慮了距離和方位信息（最大化最小角特性）的同時(shí)，Delaunay 三角剖分充分利用了每個(gè)點(diǎn)的信息，該特性使得采用少量的信息便可以刻畫更為精細(xì)的局部特征，且使用更為簡便，根據(jù)格網(wǎng)點(diǎn)速度場可以快速得到該區(qū)域內(nèi)任意未知點(diǎn)的速度值。

與中國大陸整體速度場模型相比，區(qū)域速度場模型更加符合局部區(qū)域地殼運(yùn)動(dòng)特征，因此針對(duì)每個(gè)地區(qū)建立適用自身的高精度水平速度場模型對(duì)于現(xiàn)代高精度大地測量工作非常有必要。文中提出的基于Delaunay 三角網(wǎng)的反距離加權(quán)格網(wǎng)模型精度高，且計(jì)算相對(duì)簡單，使用較為方便，值得進(jìn)一步借鑒和推廣。