應用擬合推估法建立我國大陸地殼水平運動速度場的研究

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

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

一、引 言

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

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

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

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

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

相應的誤差方程為

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

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

構造如下極值條件

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

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

三、板塊運動的趨勢項

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

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

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

四、協方差函數的擬合

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

協方差函數的Gauss指數函數為

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

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

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

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

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

表1 不顧及噪聲的協方差函數系數

表2 顧及噪聲的協方差函數系數

圖1 不顧及噪聲的Gauss函數

圖2 不顧及噪聲的Hirvonen函數

圖3 顧及噪聲的Gauss函數

圖4 顧及噪聲的Hirvonen函數

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

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

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

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

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

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

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

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

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

表3 殘差統計 mm/a

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

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

2)利用實測速度原始觀測數據(即不經過去中心化處理)，在建立協方差函數中無論考慮觀測噪聲與否，是否考慮板塊運動模型，擬合推估解效果都很差(見方案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。這是因為在擬合推估模型中要求對數據進行去中心化處理，即去掉傾向性參數的影響，否則不消去傾向，就會產生畸變。

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

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

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

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

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

[4]陳俊勇．中國現代大地基準：中國大地坐標系統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國家大地控制網的構建和它的技術進步[J]．測繪學報，2007，36(1)：1-8．

[8]王敏，張祖勝，許明元，等．2000國家GPS大地控制網的數據處理和精度評估[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]．地球物理學報，2010，53(5)：1109-1117．

[12]劉經南，姚宜斌，施闖．中國地殼運動整體速度場模型的建立方法研究[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．