應(yīng)用擬合推估法建立我國(guó)大陸地殼水平運(yùn)動(dòng)速度場(chǎng)的研究

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

(1．信息工程大學(xué)地理空間信息學(xué)院，河南鄭州450052;2．地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西西安710054;3．西安測(cè)繪研究所，陜西 西安710054)

一、引 言

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

本文嘗試將歐拉矢量作為函數(shù)模型的趨勢(shì)項(xiàng)，將局部形變作為函數(shù)模型的信號(hào)項(xiàng)，從隨機(jī)信號(hào)協(xié)方差函數(shù)的建立著手，以建立與觀測(cè)向量的隨機(jī)誤差相匹配的隨機(jī)信號(hào)向量協(xié)方差函數(shù)，通過(guò)建立相對(duì)合理的函數(shù)模型和隨機(jī)模型來(lái)構(gòu)建中國(guó)大陸坐標(biāo)速度場(chǎng)模型。

二、板塊運(yùn)動(dòng)的擬合推估模型

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

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

相應(yīng)的誤差方程為

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

其相應(yīng)的誤差方程也改寫(xiě)為

構(gòu)造如下極值條件

如果同時(shí)考慮未測(cè)點(diǎn)信號(hào)^S'，保留 ^S'與 ^S相關(guān)性，則式(6)的擬合推估解為

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

三、板塊運(yùn)動(dòng)的趨勢(shì)項(xiàng)

理論上，點(diǎn)位坐標(biāo)變化速度與板塊運(yùn)動(dòng)、地心運(yùn)動(dòng)等地球動(dòng)力學(xué)因素有關(guān)。點(diǎn)位坐標(biāo)的變化可分解為水平變化和高程變化，但是由于目前用于監(jiān)測(cè)點(diǎn)位變化的空間技術(shù)在垂直方向的精度較低，在擬合推估時(shí)若采用三維地心坐標(biāo)，可能會(huì)因?yàn)榇瓜蛘舅俣日`差的影響點(diǎn)位變化的擬合推估精度，因此實(shí)踐中僅考慮站心速度的水平分量。若僅考慮板塊運(yùn)動(dòng)對(duì)點(diǎn)位坐標(biāo)的影響，則站坐標(biāo)的水平運(yùn)動(dòng)與板塊運(yùn)動(dòng)的歐拉矢量的關(guān)系為：

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

選擇不同的板塊運(yùn)動(dòng)物理模型可計(jì)算得到點(diǎn)位的不同板塊運(yùn)動(dòng)速度，即測(cè)站的趨向性參數(shù)。從點(diǎn)位實(shí)測(cè)速度扣除板塊運(yùn)動(dòng)物理模型獲得的速度，則點(diǎn)位剩余速度為

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

在擬合推估中，關(guān)鍵的問(wèn)題是在數(shù)據(jù)處理前已知協(xié)方差矩陣ΣS、ΣSS'，而協(xié)方差矩陣中的元素是根據(jù)協(xié)方差函數(shù)計(jì)算的?；诘貧みB續(xù)變形的假設(shè)，在空間上很靠近的兩個(gè)站點(diǎn)的速度值必然相關(guān)性強(qiáng)，反之則弱?？梢栽O(shè)地殼運(yùn)動(dòng)場(chǎng)兩點(diǎn)間的協(xié)方差為距離的連續(xù)函數(shù)，隨機(jī)信號(hào)協(xié)方差函數(shù)可采用Gauss指數(shù)函數(shù)、Hirvonen函數(shù)等。協(xié)方差函數(shù)的Hirvonen函數(shù)為

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

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

一般來(lái)講，信號(hào)的初始方差可按下式計(jì)算

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

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

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

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

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

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

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

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

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

五、點(diǎn)位坐標(biāo)速度的最小二乘擬合推估解

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

1)不考慮隨機(jī)信號(hào)的板塊運(yùn)動(dòng)模型最小二乘解;

2)利用實(shí)測(cè)速度不考慮觀測(cè)噪聲建立信號(hào)協(xié)方差函數(shù)，不考慮板塊運(yùn)動(dòng)模型的擬合推估解;

3)利用實(shí)測(cè)速度考慮觀測(cè)噪聲建立信號(hào)協(xié)方差函數(shù)，不考慮板塊運(yùn)動(dòng)模型的擬合推估解;

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

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

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

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

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

表4 29個(gè)基準(zhǔn)站外部檢查精度 mm/a

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

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

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

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

[1]李延興，張靜華，李智，等．由GPS網(wǎng)融合得到的中國(guó)大陸及周邊地區(qū)的地殼水平運(yùn)動(dòng)[J]．測(cè)繪學(xué)報(bào)，2003，32(4)：301-307．

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

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

[4]陳俊勇．中國(guó)現(xiàn)代大地基準(zhǔn)：中國(guó)大地坐標(biāo)系統(tǒng)2000(CGCS2000)及其框架[J]．測(cè)繪學(xué)報(bào)，2008，37(3)：269-271．

[5]魏子卿．2000中國(guó)大地坐標(biāo)系[J]．大地測(cè)量與地球動(dòng)力學(xué)，2008，28(6)：1-5．

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

[7]陳俊勇，楊元喜，王敏，等．2000國(guó)家大地控制網(wǎng)的構(gòu)建和它的技術(shù)進(jìn)步[J]．測(cè)繪學(xué)報(bào)，2007，36(1)：1-8．

[8]王敏，張祖勝，許明元，等．2000國(guó)家GPS大地控制網(wǎng)的數(shù)據(jù)處理和精度評(píng)估[J]．地球物理學(xué)報(bào)，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]柴洪洲，催岳，明鋒．最小二乘配置方法確定中國(guó)大陸主要塊體運(yùn)動(dòng)模型[J]．測(cè)繪學(xué)報(bào)，2009，38(1)：61-65．

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

[12]劉經(jīng)南，姚宜斌，施闖．中國(guó)地殼運(yùn)動(dòng)整體速度場(chǎng)模型的建立方法研究[J]．武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2002，27(4)：331-336．

[13]魏子卿，劉光明，吳富梅．2000中國(guó)大地坐標(biāo)系：中國(guó)大陸速度場(chǎng)[J]．測(cè)繪學(xué)報(bào)，2010，39(4)：403-410．

[14]蔣志浩，張鵬，秘金鐘，等．基于CGCS2000的中國(guó)地殼水平運(yùn)動(dòng)速度場(chǎng)模型研究[J]．測(cè)繪學(xué)報(bào)，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]楊元喜，曾安敏，吳富梅．基于歐拉矢量的中國(guó)大陸地殼水平運(yùn)動(dòng)自適應(yīng)擬合推估模型[J]．中國(guó)科學(xué)：地球科學(xué)，2011，41(8)：1116-1125．