曾安敏，劉光明，秦顯平，唐穎哲

(1．信息工程大學地理空間信息學院，河南鄭州450052;2．地理信息工程國家重點實驗室，陜西西安710054;3．西安測繪研究所，陜西 西安710054)

一、引 言

中國大陸除遭受印度板塊的強烈碰撞以及菲律賓板塊、太平洋板塊的俯沖影響外，板塊內(nèi)形變也非常復雜和不均勻[1-3]。在具有如此復雜的地殼運動背景的大陸構建大地測量坐標框架，必須顧及板塊的整體運動和局部運動。我國已建立了2000中國坐標系CGCS2000[4-6]，2008年8月正式啟用相應 的 參 考 框 架 CTRF2000[6-7]。 但 由 于 建 立CTRF2000時的資料截止到2002年，重復觀測時間較短，絕大部分框架點未賦予點位速度場信息，無法反映我國坐標框架的動態(tài)變化，不能完整地表述我國坐標系統(tǒng)的現(xiàn)勢性[8-9]。由于板塊運動引起的坐標變化是不容忽視的，速度場數(shù)據(jù)成為現(xiàn)代大地測量不可或缺的基礎性數(shù)據(jù)。地殼運動的描述往往采用歐拉矢量法和統(tǒng)計擬合的方法。基于實測速度信息可以采用最小二乘準則求得歐拉矢量的估值，常用的統(tǒng)計擬合方法有多面函數(shù)法、有限元插值法、擬合推估法等。近10年來，國內(nèi)不少學者在2002年前利用“網(wǎng)絡工程”數(shù)據(jù)陸續(xù)開展我國大陸板塊運動相關研究工作[10-12]，但由于使用數(shù)據(jù)少，其成果的實用價值有一定的局限性。2009年前利用“網(wǎng)絡工程”數(shù)據(jù)研究了最新的速度場，有學者利用局部歐拉矢量法建立了3度格網(wǎng)平均值[13]，但用于建立局部歐拉矢量的點及其分布是其關鍵。也有學者利用有限元法建立了速度模型[14]，但面元的確定是關鍵。多面函數(shù)法的關鍵是節(jié)點的選擇、核函數(shù)的確定以及平滑因子的優(yōu)化，在實踐中很難實踐，需要依靠大量的試驗。最小二乘配置法具有理論上的嚴密性，可綜合顧及非隨機參數(shù)和隨機信號，在重力場的擬合與推估有成功應用[15]。在速度場的建立中，為顧及各局部形變場的變化信息，也可采用擬合推估法，擬合推估法的關鍵是信號向量協(xié)方差函數(shù)的確定[16]。也有研究團隊利用自適應擬合推估模型[17]建立了中國大陸地殼水平運動[18]，從隨機模型誤差的角度通過方差因子部分調整協(xié)方差函數(shù)的不可靠問題，需要迭代計算。

本文嘗試將歐拉矢量作為函數(shù)模型的趨勢項，將局部形變作為函數(shù)模型的信號項，從隨機信號協(xié)方差函數(shù)的建立著手，以建立與觀測向量的隨機誤差相匹配的隨機信號向量協(xié)方差函數(shù)，通過建立相對合理的函數(shù)模型和隨機模型來構建中國大陸坐標速度場模型。

二、板塊運動的擬合推估模型

設測站的n×1維速度觀測向量為L，相應的誤差向量為e，其期望為E(e)=0，協(xié)方差矩陣為;X為3×1維的板塊運動歐拉矢量，S為具有先驗統(tǒng)計特性的u×1維信號向量，協(xié)方差矩陣為。描述研究區(qū)域速度觀測向量的空間分布性與離散度和分別為e和S的權矩陣。塊體運動的擬合推估方程為

式中，A為n×3階列滿秩設計矩陣;B為n×u維系數(shù)矩陣。

相應的誤差方程為

假如事先去除趨向性參數(shù)部分，僅保留隨機信號部分，則式(1)改寫為

其相應的誤差方程也改寫為

構造如下極值條件

如果同時考慮未測點信號^S'，保留 ^S'與 ^S相關性，則式(6)的擬合推估解為

通常情況下B=I，如果不考慮觀測噪聲，即Σe=0，則有

三、板塊運動的趨勢項

理論上，點位坐標變化速度與板塊運動、地心運動等地球動力學因素有關。點位坐標的變化可分解為水平變化和高程變化，但是由于目前用于監(jiān)測點位變化的空間技術在垂直方向的精度較低，在擬合推估時若采用三維地心坐標，可能會因為垂向站速度誤差的影響點位變化的擬合推估精度，因此實踐中僅考慮站心速度的水平分量。若僅考慮板塊運動對點位坐標的影響，則站坐標的水平運動與板塊運動的歐拉矢量的關系為：

式中，R表示地球半徑;λ、φ為點的經(jīng)緯度。

選擇不同的板塊運動物理模型可計算得到點位的不同板塊運動速度，即測站的趨向性參數(shù)。從點位實測速度扣除板塊運動物理模型獲得的速度，則點位剩余速度為

四、協(xié)方差函數(shù)的擬合

在擬合推估中，關鍵的問題是在數(shù)據(jù)處理前已知協(xié)方差矩陣ΣS、ΣSS'，而協(xié)方差矩陣中的元素是根據(jù)協(xié)方差函數(shù)計算的。基于地殼連續(xù)變形的假設，在空間上很靠近的兩個站點的速度值必然相關性強，反之則弱?？梢栽O地殼運動場兩點間的協(xié)方差為距離的連續(xù)函數(shù)，隨機信號協(xié)方差函數(shù)可采用Gauss指數(shù)函數(shù)、Hirvonen函數(shù)等。協(xié)方差函數(shù)的Hirvonen函數(shù)為

協(xié)方差函數(shù)的Gauss指數(shù)函數(shù)為

式中，C0、k為待定常數(shù);C(d)為 i、j兩點之間的協(xié)方差;d為i、j兩點之間的距離。

一般來講，信號的初始方差可按下式計算

式中，n為距離間隔為d的點對個數(shù);li、lj為距離間隔為d的點對的去中心化觀測值，即由式(10)獲得的點位剩余速度?？紤]觀測噪聲Ce(0)，同時假設噪聲e、信號S之間獨立，有

由于li、lj點對間距離dij嚴格等于距離間隔d的點對較少，因此實際計算時，可選擇落在距離間隔d一定誤差范圍Δd內(nèi)的點對即可，即li、lj點對間距離dij滿足|dij-d|=Δd的點對認為是距離間隔d的點對。然后按一定的準則估計式(11)、式(12)中的待定常數(shù) C0、k，由式(14)可計算觀測噪聲方差Ce(0)。

采用“地殼運動觀測網(wǎng)絡工程”直到2008年的觀測數(shù)據(jù)獲得的扣除板塊運動的實測速度作為信號協(xié)方差函數(shù)的觀測數(shù)據(jù)。在實際應用中，通常不考慮觀測噪聲Ce(0)，表1給出了不顧及觀測噪聲的不同協(xié)方差函數(shù)系數(shù);圖1、圖2為不顧及噪聲的距離-協(xié)方差分布。表2給出了顧及觀測噪聲的不同協(xié)方差函數(shù)系數(shù);圖3、圖4為顧及噪聲的距離-協(xié)方差分布。可以看出，顧及觀測噪聲時擬合的協(xié)方差函數(shù)與實測速度的擬合程度高。另外，由于按站間距離分檔計算中短距離的站點分布較少，噪聲偏大，使協(xié)方差分布在中短距離的區(qū)間吻合程度降低。

表1 不顧及噪聲的協(xié)方差函數(shù)系數(shù)

表2 顧及噪聲的協(xié)方差函數(shù)系數(shù)

圖1 不顧及噪聲的Gauss函數(shù)

圖2 不顧及噪聲的Hirvonen函數(shù)

圖3 顧及噪聲的Gauss函數(shù)

圖4 顧及噪聲的Hirvonen函數(shù)

五、點位坐標速度的最小二乘擬合推估解

本文采用“地殼運動觀測網(wǎng)絡工程”直到2008年的觀測數(shù)據(jù)進行點位坐標速度擬合，包括29個連續(xù)運行基準站和站速度優(yōu)于3mm/a的1041個重復觀測站。由于29個連續(xù)運行基準站大都建設在基巖上，其站速度是經(jīng)過長期觀測數(shù)據(jù)計算得到的，具有較高的精度。因此，把這29個連續(xù)運行基準站作為外部檢查點，只用來評定所構建的速度場模型的精度。采用1041個重復觀測站構建速度場模型，共設計了6種計算方案：

1)不考慮隨機信號的板塊運動模型最小二乘解;

2)利用實測速度不考慮觀測噪聲建立信號協(xié)方差函數(shù)，不考慮板塊運動模型的擬合推估解;

3)利用實測速度考慮觀測噪聲建立信號協(xié)方差函數(shù)，不考慮板塊運動模型的擬合推估解;

4)利用實測速度經(jīng)過板塊運動模型中心化后的剩余速度，在不考慮觀測噪聲情況下建立信號協(xié)方差函數(shù)(即表1中的高斯函數(shù)參數(shù))，無系統(tǒng)性參數(shù)的擬合推估解，即采用(4)式;

5)利用實測速度經(jīng)過板塊運動模型中心化后的剩余速度，在考慮觀測噪聲情況下建立信號協(xié)方差函數(shù)(即表2中的高斯函數(shù)參數(shù))，無系統(tǒng)性參數(shù)的擬合推估解，即采用(4)式;

6)在方案5基礎上，考慮板塊運動模型的擬合推估解，即利用表2中的高斯函數(shù)，按式(2)進行計算。

各種方案的殘差統(tǒng)計見表3，外部檢查結果見表4。分析上述結果，可以得出以下結論。

表3 殘差統(tǒng)計 mm/a

表4 29個基準站外部檢查精度 mm/a

1)從殘差看，基于板塊運動模型的最小二乘解的殘差非常大，最大達到30 mm/a;統(tǒng)計精度也差，東方向4．63為mm/a，北方向6．19為 mm/a。從外部檢查點來看，歐拉矢量最小二乘解的外部精度中，東方向為3．92mm/a，北方向6．59 為mm/a。這是因為歐拉矢量最小二乘解只確定了傾向性部分即參數(shù)部分，并未顧及局部的系統(tǒng)性變化。

2)利用實測速度原始觀測數(shù)據(jù)(即不經(jīng)過去中心化處理)，在建立協(xié)方差函數(shù)中無論考慮觀測噪聲與否，是否考慮板塊運動模型，擬合推估解效果都很差(見方案2、3。方案2外部檢查RMS在北方向為5．07mm/a，東方向為6．33mm/a;最大值在北方向為22．90mm/a，東方向為 18．68mm/a。方案 3 外部檢查 RMS在北方向為 3．85 mm/a，東方向為7．21mm/a;最大值在北方向為 9．70 mm/a，東方向為18．79 mm/a。這是因為在擬合推估模型中要求對數(shù)據(jù)進行去中心化處理，即去掉傾向性參數(shù)的影響，否則不消去傾向，就會產(chǎn)生畸變。

3)經(jīng)過去中心化數(shù)據(jù)處理后得到剩余速度，在建立協(xié)方差函數(shù)中考慮觀測噪聲與否，對擬合推估解影響很大(見方案4、5)。不考慮觀測噪聲，方案4外部檢查RMS在北方向為4．45 mm/a，東方向為4．25mm/a;而考慮觀測噪聲的方案5外部檢查RMS 在北方向為 1．93mm/a，東方向為 1．76mm/a。綜合方案2、3結果，說明傾向性參數(shù)是存在的，與板塊整體運動這一物理現(xiàn)象是一致的，同時也說明觀測噪聲是存在的。

4)考慮觀測噪聲，利用去中心化數(shù)據(jù)處理后得到剩余速度建立協(xié)方差函數(shù)，無論是先計算趨向性參數(shù)部分(方案5)，還是在模型中考慮趨向性參數(shù)部分(方案6)，計算結果精度一致，精度都非常高，內(nèi)外部精度一致，北方向和東方向都優(yōu)于2mm/a。這是因為建立的函數(shù)模型和隨機模型都相對合理，在函數(shù)模型方面，同時考慮了傾向性參數(shù)部分和隨機信號。在隨機模型方面，在建立隨機信號協(xié)方差函數(shù)時經(jīng)過了中心化數(shù)據(jù)處理，并顧及了觀測噪聲，這使得信號協(xié)方差函數(shù)與觀測噪聲協(xié)方差具有較好的協(xié)調性。

[1]李延興，張靜華，李智，等．由GPS網(wǎng)融合得到的中國大陸及周邊地區(qū)的地殼水平運動[J]．測繪學報，2003，32(4)：301-307．

[2]朱文耀，王小丫，符養(yǎng)，等．基于ITRF2000的全球板塊運動模型和中國的地殼形變[J]．地球物理學報，2002，45(S0)：197-204．

[3]王敏，沈正康，牛之俊，等．現(xiàn)今中國大陸地殼運動與活動塊體模型[J]．中國科學(D 輯)，2003，33(B04)：21-32．

[4]陳俊勇．中國現(xiàn)代大地基準：中國大地坐標系統(tǒng)2000(CGCS2000)及其框架[J]．測繪學報，2008，37(3)：269-271．

[5]魏子卿．2000中國大地坐標系[J]．大地測量與地球動力學，2008，28(6)：1-5．

[6]YANG Yuanxi．Chinese Geodetic Coordinate System 2000[J]． Chinese Science Bulletin， 2009，54(15)：2714-2721．

[7]陳俊勇，楊元喜，王敏，等．2000國家大地控制網(wǎng)的構建和它的技術進步[J]．測繪學報，2007，36(1)：1-8．

[8]王敏，張祖勝，許明元，等．2000國家GPS大地控制網(wǎng)的數(shù)據(jù)處理和精度評估[J]．地球物理學報，2005，48(4)：817-823．

[9]YANG Yuanxi，TANG Yungzhe，CHEN Chuanlu，et al．Integrated Adjustment of Chinese 2000’GPSControl Network[J]．Survey Review，2009，41(313)：226-237．

[10]柴洪洲，催岳，明鋒．最小二乘配置方法確定中國大陸主要塊體運動模型[J]．測繪學報，2009，38(1)：61-65．

[11]江在森，劉經(jīng)南．應用最小二乘配置建立地殼運動速度場與應變場的方法[J]．地球物理學報，2010，53(5)：1109-1117．

[12]劉經(jīng)南，姚宜斌，施闖．中國地殼運動整體速度場模型的建立方法研究[J]．武漢大學學報：信息科學版，2002，27(4)：331-336．

[13]魏子卿，劉光明，吳富梅．2000中國大地坐標系：中國大陸速度場[J]．測繪學報，2010，39(4)：403-410．

[14]蔣志浩，張鵬，秘金鐘，等．基于CGCS2000的中國地殼水平運動速度場模型研究[J]．測繪學報，2009，38(6)：471-476．

[15]TSCHERNING C C．Collocation and Least Squares Methods as a Tool for Handling Gravity Field Dependent Data Obtained Through Space Research Techniques[R]．Paris：[s．n．]，1978：141-149．

[16]YANG Yuanxi．Robustfying Collocation[J]．Manusctipta Geodactica．1992，17(1)：21-28．

[17]Yang Y，Zeng A，Zhang J．Adaptive Collocation with Application in Height System Transformation[J]．Journal of Geodesy，2009，83(5)：403-410．

[18]楊元喜，曾安敏，吳富梅．基于歐拉矢量的中國大陸地殼水平運動自適應擬合推估模型[J]．中國科學：地球科學，2011，41(8)：1116-1125．