考慮區(qū)域構(gòu)造特征的地殼形變分析擬合推估模型

趙麗華，楊元喜，王慶良

1.長安大學地質(zhì)工程與測繪學院，陜西西安710054；2.西安測繪研究所，陜西西安710054；3.中國地震局第二監(jiān)測中心，陜西西安710054

1 引 言

最小二乘擬合推估，又稱最小二乘配置，根據(jù)最小二乘推估來內(nèi)插和外推重力異常［1］。近年來，國內(nèi)外眾多學者基于擬合推估法開展大量的地殼形變分析研究。文獻［2—3］將其用于地殼垂直形變速度場的研究，認為該方法是研究構(gòu)造變形的有利工具；文獻［4］應(yīng)用擬合推估方法研究日本東京地區(qū)水平運動形變場，驗證該方法可以排除觀測噪聲干擾，分離出所需要的形變信號；文獻［5］對華北盆地154個GPS觀測速度運用擬合推估方法得到連續(xù)的速度場，計算該地區(qū)的應(yīng)變分布，并與傳統(tǒng)的離散三角法相比較，得出擬合推估方法能更清楚地反映出地殼的形變狀態(tài)；文獻［6—9］在地殼形變分析時均采用擬合推估方法內(nèi)插位移向量來計算應(yīng)變場；文獻［10］提出一種改進的擬合推估法，根據(jù)趨勢項的梯度迭代計算協(xié)方法差函數(shù)，從而提取非隨機形變場的趨勢項與形變信號。這些文獻將擬合推估法應(yīng)用于地殼形變分析時，大都假設(shè)被擬合參數(shù)在空間是連續(xù)變化的。事實上對于大尺度、地殼運動復(fù)雜的形變場，由于小斷層或隱伏斷層的存在，形變的連續(xù)介質(zhì)中常常存在著若干非連續(xù)面，若仍按統(tǒng)一的擬合推估模型進行形變分析，則會扭曲形變分析結(jié)果，因為統(tǒng)一的擬合推估模型只能給出形變（應(yīng)變）場在空間域的總體變化趨勢，在斷裂附近只能給出一個連續(xù)的光滑結(jié)果。如此，很難合理地描述不同活動塊體間的運動（或變形）差異。

一般認為，由斷層分割的各區(qū)域塊體的運動方式和大小具有整體上一致趨勢。但由于觀測異常誤差或塊體內(nèi)部非均勻變形，使同一塊體內(nèi)有些觀測點上的位移觀測值表現(xiàn)得“不合群”，這些點會導(dǎo)致形變參數(shù)解有偏，一般不宜參與塊體整體運動參數(shù)計算。由于局部非均勻變形和干擾會影響塊體實際的運動和變形參數(shù)精度，采用高崩潰污染率抗差估計方法［11］，將這些顯著變形或受異常誤差干擾的觀測值判別篩選出來。篩選之后的點組即為相對穩(wěn)定點組，如此建立的塊體整體運動模型，可望較合理地反映塊體整體運動趨勢。

2 相對穩(wěn)定點組的確定

如引言所述，在建立塊體的速率模型時，為提高模型的可靠性，同一塊體內(nèi)部具有異常形變的點不宜參與建模，一般將這些異常觀測和非均勻變形看作對所在塊體整體位移和變形的一種干擾，或者說是“粗差”將其剔除［12－13］。這里采用抗差估計的方法來確定塊體內(nèi)運動相對平穩(wěn)的點。但一般抗差估計方法在迭代計算中，參數(shù)的初值一般為相應(yīng)的最小二乘解。當觀測值出現(xiàn)異常時，形變參數(shù)的最小二乘解會有較大偏離，從而使后繼的殘差值偏離，等價權(quán)也隨之失效。為此，尋求一種崩潰污染率較高的參數(shù)估值，即當觀測值中有較多異常或較大異常時，參數(shù)估值不至受到太大影響，從而為后續(xù)的抗差等價權(quán)的確定提供可靠的測值殘差?？刹捎萌缦滤惴ǎ?/p>

（1）將各點速率觀測值分別取中位數(shù)med（L），求各速率值與相應(yīng)中位數(shù)之差，即初始殘差

（2）計算ΔL的均方差因子抗差解［14－15］

基于此式求得的變形參數(shù)抗差解具有50%的高崩潰污染率［11］。

（3）取強淘汰函數(shù)

式中，k0可取1.0或1.5。

（5）進行后繼M估計迭代計算?？紤]到觀測速率的相對穩(wěn)定性，權(quán)函數(shù)IGG方案做了如下調(diào)整［11］

式中，k＝1.5～2。

迭代結(jié)束后等價權(quán)為0的觀測值即為異常觀測，其余點即為相對穩(wěn)定點組參加后繼的分區(qū)域擬合推估模型的建立。

3 塊體運動模型的精化及相鄰塊體邊界的確定

對于大尺度、地殼運動復(fù)雜、小斷層或隱伏斷層比較多的形變場，若按通常的擬合推估方法建模，對于分布在斷層兩側(cè)的點，只能根據(jù)速度場在空間域的總體變化趨勢進行平滑處理，其結(jié)果不能正確地反映斷層附近形變的真實情況。因此，根據(jù)所研究區(qū)域的地質(zhì)構(gòu)造特征，建立區(qū)域運動擬合推估模型。

（1）利用已有地質(zhì)資料按斷層分布將研究區(qū)域劃分為若干塊體。相鄰斷層附近界限不明的點不參與劃分，等候判別。按前述方法找出各區(qū)域內(nèi)相對穩(wěn)定的點組，利用擬合推估原理建立如下區(qū)域運動模型［16－17］

式中，L為沿坐標軸的速度分量；A、B為設(shè)計矩陣；β為非隨機參數(shù)向量；Δ為觀測誤差向量，Δ～N（0，ΣΔ）；S為已測點上隨機信號向量，S～N（0，ΣS）；ΣΔS＝0，即Δ與隨機參數(shù)S不相關(guān)。Aβ為系統(tǒng)形變，反映塊體運動的一般趨勢；BS為隨距離變化的隨機函數(shù)，反映塊體內(nèi)連續(xù)的、不規(guī)則的形變。這里S是以單點信號的形式出現(xiàn)的，所以B為單位陣。為推估未測點上信號，將式（5）改寫為

式中，S′為推估點上信號，S′～N（0，ΣS′），且基于最小二乘原則可得

式中，ΣL＝ΣΔ＋BΣSBT，為觀測值協(xié)方差矩陣。任一未測點上信號則為

（2）根據(jù)測量誤差的性質(zhì)，設(shè)V服從正態(tài)分布，分別對不同分塊中的V進行區(qū)間估計。給定置信概率為1－α＝0.90，則V的置信下限Vlow和置信上限Vup分別為［18］

則（Vlow，Vup）即為模型殘差的置信區(qū)間。

（3）選取邊界附近的點進行判別。將斷層附近分界不明的點作為檢驗點，分別用各區(qū)域的擬合推估模型推估其形變量，計算預(yù)測殘差，即推估值與觀測值較差

式中，LBi表示邊界上點Bi的觀測值；ABi、＾S′Bi分別表示該點對應(yīng)的形變參數(shù)系數(shù)陣和信號。

（4）假設(shè)檢驗

若H0成立，則將Bi點添加到相應(yīng)塊體的穩(wěn)定點組。

（5）用新組成的穩(wěn)定點組重復(fù)（1）～（4）的計算。

迭代的過程即是對形變模型精化的過程，同時也確定了斷層邊界的分布。

4 計算與分析

4.1 算例1

模擬一個100km×100km的垂直位移形變場，其內(nèi)部100個監(jiān)測點坐標由隨機函數(shù)生成；同時假設(shè)一條斷層穿過該形變場，將形變場分割為區(qū)域1和區(qū)域2兩部分（圖1）。給位于區(qū)域1的監(jiān)測點賦予2mm／y的速率，位于區(qū)域2的監(jiān)測點賦予1mm／y的速率，并且加入均方差為1× 10－6的正態(tài)隨機誤差。把隨機誤差和固定速率相加得到監(jiān)測點垂直位移速率的模擬觀測值。若考慮由已有地質(zhì)資料所判斷的斷層的位置可能存在的不確定性，將模擬斷層附近的20個點作為分界不明的點等候判別；另外，為檢核擬合推估的精度，在區(qū)域1、區(qū)域2中隨機均勻選取19個點作為外部檢核點，檢核精度的均方根誤差公式為

式中，m為檢核點個數(shù)；di為檢核點不符值，是擬合模型模擬的檢核點速率與該點實測速率之差。

圖1 模擬觀測點及斷裂分布Fig.1 The distribution of simulated measurement points and fault

為比較分析高崩潰污染率抗差估計法確定相對穩(wěn)定點組的有效性，試驗采用兩種方案進行。

方案1：無異常數(shù)據(jù)影響情形。分別采用一般擬合推估方法；根據(jù)地質(zhì)資料的分區(qū)域擬合推估法；基于穩(wěn)定點組的分區(qū)域抗差擬合推估法。

方案2：有異常數(shù)據(jù)影響情形。各區(qū)域分別隨機選取四個點人為加入2mm～3mm的異常值，進行與方案1完全相同的計算。

圖2、圖3分別表示兩個方案中不同方法計算得到的垂直形變場空間分布圖；不同方法估計的檢核點均方根誤差列于表1；圖4為兩種方案中各擬合模型檢合點不符值分布。

表1 不同方案計算結(jié)果比較Tab.1 Comparison of different methods mm

分析上述計算結(jié)果可發(fā)現(xiàn)：

（1）一般擬合推估方法得到的形變場在斷層附近只是將非連續(xù)面間的點值進行光滑，沒有合理反映形變的真實情況，使推估的形變量有較大偏差（圖2（a）），而分區(qū)域擬合推估法（圖2（b））和分區(qū)域抗差擬合推估法（圖2（c））在方案1沒有異常數(shù)據(jù)的影響下都較好地反映出斷裂活動及塊體間運動差異，較真實地反映區(qū)域的形變特征。

（2）當觀測存在異常時，一般擬合推估法（圖3（a））和分區(qū)域擬合推估法（圖3（b））受異常數(shù)據(jù)的影響非常大，由這兩種方法得到的形變場空間分布圖形與圖2（b）和圖2（c）中相應(yīng)圖形相比，存在顯著偏差；而分區(qū)域抗差擬合推估法（圖3（c））仍然較好地反映研究區(qū)域真實的形變特征，與方案1圖形非常一致，反映出該方法良好地抵制異常數(shù)據(jù)影響的能力。

（3）從外部檢核結(jié)果可以看出（圖4），一般擬合推估方法推估的未測點形變量與其真實形變差異較大；而不論是一般擬合推估還是分區(qū)域擬合推估，在異常數(shù)據(jù)的影響下（方案2），檢核點不符值變化都比較大，反映出未測點推估值的誤差較大。分區(qū)域抗差擬合推估則較好地抵制異常值的影響，檢核點不符值整體上較為均勻地分布于0值附近，均方根誤差明顯降低。

（4）從表1中也可以看出，無論有無異常污染，分區(qū)域抗差擬合推估解都基本一致，外部檢核點均方根誤差都比較小，且大致相等，反映出該方法具有較好的穩(wěn)定性。這說明當采用高崩潰污染率抗差估計方法對觀測值中的異常值進行篩選后，分區(qū)域擬合推估模型能更好地逼近未測點的真實值。

圖2 方案1中不同計算方法得到的垂直形變場空間分布圖Fig.2 Vertical velocity fields by different methods in case 1

圖3 方案2中不同方法得到的垂直形變場空間分布圖Fig.3 Vertical velocity fields by different methods in case 2

圖4 各方案檢核點不符值Fig.4 Residuals of check points

4.2 算例2

實測數(shù)據(jù)為100°E～110°E，33°N～37°N區(qū)域中100個GPS站點的速率及相應(yīng)中誤差，點位分布及速率見圖5。該數(shù)據(jù)由中國地震局第二地形變中心于1999－07月和2001－07月分別采用Z－型GPS接收機觀測所得，該觀測網(wǎng)位于青藏高原東北緣，是國家重大科學工程“中國地殼運動觀測網(wǎng)絡(luò)”的一部分。對該GPS網(wǎng)進行每站96h蛙跳式觀測，相應(yīng)數(shù)據(jù)與中國大陸及周邊約15個IGS站進行同步解算，分別采用GAMIT軟件和GLOBK軟件進行統(tǒng)一處理。

圖5 各點位置分布及速率Fig.5 Distribution and velocity of observation points

根據(jù)晚第四紀活動斷裂帶分布情況，固關(guān)—功縣斷裂帶從該區(qū)域通過［19］。根據(jù)此斷裂帶對該區(qū)域概略劃分為A區(qū)、B區(qū)，對分布于斷裂帶附近的19個點進行檢驗判別。對該區(qū)域形變模型的建立采用如下兩種方案：① 一般擬合推估法；②本文提出的分區(qū)域抗差擬合推估法。

為說明新方法的有效性，將兩種方案得到的模型擬合誤差列于表2，同時將邊界點判別歸區(qū)前后各區(qū)的擬合誤差列于表3。擬合誤差由擬合殘差的均方根誤差表示

式中，V為模型擬合殘差向量；m為擬合點的個數(shù)。

外符合精度檢核則在區(qū)域中均勻選取27個點作為外部檢核點不參與計算，通過以上兩種方案分別模擬其速率值。由這些點的模擬值與已知值（觀測速率）之差，根據(jù)式（14）計算檢核點的均方根誤差并列于表4；將不同計算方案模擬的各檢核點水平運動速度與相應(yīng)實測值進行比較，比較結(jié)果見圖6。

表2 不同方案模型擬合誤差MTab.2 Mof different cases mm

表3 邊界點判別前后各區(qū)擬合誤差MTab.3 Min area Aand B mm

圖6 算例2中不同方案模擬的水平運動與實測值比較Fig.6 Comparison between simulated velocity and observed velocity

從計算結(jié)果中可以看出：

（1）分區(qū)域抗差擬合推估模型的擬合誤差明顯小于一般擬合推估模型（表2），說明該模型能更合理的描述此類區(qū)域的變形特征。

（2）受不同區(qū)域運動差異影響，邊界上的點在判別歸區(qū)之前，各區(qū)的擬合RMS都比較大，而經(jīng)過合理的歸區(qū)后，擬合均方根誤差明顯減?。ū?）。實際上，在判別歸區(qū)前，邊界點在各區(qū)的殘差表現(xiàn)出明顯的區(qū)域特征，如有7個點在A區(qū)中具有較大殘差而在B區(qū)中殘差較小，同時A區(qū)中殘差較小的8個點在B區(qū)中具有較大誤差，體現(xiàn)顯著檢驗的效果。

（3）一般擬合推估方法，由于要反映形變場在空間域的整體變化趨勢，對不同塊體運動差異的描述顯得較為粗略，由此求得的點位形變的精度相對較低（圖6）；分區(qū)域抗差擬合推估法由于顧及不同塊體各自的運動特征，擬合效果明顯提高，在一定程度上反映斷裂活動和塊體運動差異，求得的點位形變與各點的實測值具有較好的一致性（圖6）。

（4）由于篇幅限制，沒有一一列出各檢核點不符值（擬合模型估計的檢核點速率與該點實測速率之差），但由檢核點不符值計算的RMS也反映出分區(qū)域抗差擬合推估法估計的各點水平運動速率更接近真實值（表4）。

（5）分區(qū)域抗差擬合推估法在確定塊體運動模型的同時，還確定塊體邊界分布（圖5中用測站點的不同顏色來反映區(qū)域劃分邊界）。這與該地區(qū)地質(zhì)情況和塊體劃分一致［19］。

表4 不同方案檢核點RMSTab.4 RMS of check points mm

5 結(jié) 論

若形變場中有隱伏斷層或存在較顯著形變差異時，有時無法確切地知道差異運動的分界帶。這時若采用一般擬合推估方法建立形變場，在斷層附近會產(chǎn)生較大偏差，無法可靠地反映斷層分布情況。本文首先通過地質(zhì)資料進行簡單判斷，對可能存在斷層的區(qū)域進行概略劃分，然后通過高崩潰污染率抗差估計法挑選各區(qū)域內(nèi)相對穩(wěn)定的點，建立各區(qū)域形變場置信區(qū)間。對分布在斷層附近分界不明點的歸屬進行判斷，確定各塊體差異運動的邊界。邊界明確后，分別建立各區(qū)域能夠描述本塊體運動趨勢的擬合推估模型，計算區(qū)域形變場。試驗算例表明，該方法在發(fā)現(xiàn)隱伏斷層并確定其分布邊界的同時，建立的形變場較好地表現(xiàn)出不同運動塊體邊界上的差異運動，具有較高推估精度。當研究區(qū)域觀測點數(shù)量相對較少時，利用此方法解得的形變參數(shù)可以推算任意未測點的形變量，達到分析形變和應(yīng)變的目的。

［1］ KRARUP T.Mathematical Foundation of Geodesy：Selected Papers of Torben Krarup［M］.Berlin：Springer，2006：29－90.

［2］ FUJII Y，XIA S.Estimation of Distribution of the Rates of Vertical Crustal Movements in the Tokai District with the Aid of Least－squares Prediction［J］.Journal of Physics of the Earth，1993，41（4）：239－256.

［3］ EL－FIKY G S，KATO T，F(xiàn)UJII Y.Distribution of the Vertical Crustal Movement Rates in the Tohoku District，Japan，Predicted by Least－squares Collocation［J］.Journal of Geodesy，1997，71（7）：432－442.

［4］ EL－FIKY G S，KATO T.Continuous Distribution of the Horizontal Strain in the Tohoku District，Japan，Predicted by Least－squares Collocation［J］.Journal of Geodynamics，1999，27（2）：213－236.

［5］ WU J C，TANG H W，CHEN Y Q.The Current Strain Distribution in the North China Basin of Eastern China by Least－square Collocation［J］.Journal of Geodynamics，2006，41（5）：462－470.

［6］ HOLLENSTEIN C，MULLER M D，GEIGER A，et al.Crustal Motion and Deformation in Greece from a Decade of GPS Measurements，1993－2003［J］.Tectonophysics，2008，449（1）：17－40.

［7］ KAHLE H G，MULLER M V，GEIGER A，et al.The Strain Field in Northwestern Greece and the Ionian Islands：Results Inferred from GPS Measurements［J］.Tectonophysics，1995，249（1）：41－52.

［8］ KAHLE H G，COCARD M，PETER Y，et al.GPS－derived Strain Rate Field within the Boundary Zones of the Eurasian，African and Arabian Plates［J］.Journal Geophysical Research，2000，105（10）：23353－23370.

［9］ JIANG Zaisen，MA Zongjin，ZHANG Xi，et al.Horizontal Strain Field and Tectonic Deformation of China Mainland Revealed by Preliminary GPS Result［J］.Journal of Geodesy and Geodynamics，2003，46（3）：353－358.（江在森，馬宗晉，張希，等.GPS初步結(jié)果揭示中國大陸水平應(yīng)變場與構(gòu)造變形［J］.地球物理學報，2003，46（3）：353－358.）

［10］ EGLI R，GEIGER A，WIGET A，et al.A Modified Least－squares Collocation Method for the Determination of Crustal Deformation：First Results in the Swiss Alps［J］.Geophysical Journal International，2007，168（1）：1－12.

［11］ YANG Yuanxi.Robust Estimation Theory and Its Applications［M］.Beijing：Bayi Publishing House，1993.（楊元喜.抗差估計理論及其應(yīng)用［M］.北京：八一出版社，1993.）

［12］ YANG Yuanxi.Adaptively Robust Least Squares Estimation［J］.Acta Geodaetica et Cartographica Sinica，1996，25（3）：206－211.（楊元喜.自適應(yīng)抗差最小二乘估計［J］.測繪學報，1996，25（3）：206－211.）

［13］ OU Jikun.Quasi－accurate Detection of Gross Errors（QUAD）［J］.Acta Geodaetica et Cartographica Sinica，1999，28（1）：15－20.（歐吉坤.粗差的擬準檢定法（QUAD法）［J］.測繪學報，1999，28（1）：15－20.）

［14］ YANG Y X.Robust Estimation of Geodetic Datum Transformation［J］.Journal of Geodesy，1999，73（5）：268－274.

［15］ YANG Y，CHENG M K，SHUM C K，et al.Robust Estimation of Systematic Errors of Satellite Laser Range［J］.Journal of Geodesy，1999，76（6）：345－349.

［16］ ZHOU Jiangwen.Collocation by Two Variant Methods［J］.Acta Geodaetica et Cartographica Sinica，1981，10（1）：9－12.（周江文.擬合推估的兩種解法［J］.測繪學報，1981，10（1）：9－12.）

［17］ YANG Y X.Robustfying Collocation［J］.Manuscripta Geodaetica，1992，17（1）：21－28.

［18］ WANG Rongxin.Mathematical Statistics［M］.Xi’an：Xi’an Jiaotong University Press，1986.（汪榮鑫.數(shù)理統(tǒng)計［M］.西安：西安交通大學出版社，1986.）

［19］ GUO Shunmin，JIANG Zaisen，ZHANG Chongli.Division and Motion Status of Blocks for the Northeastern Tibetan Plateau in Late Quaternary［J］.Seismology and Geology，2009，22（3）：219－231.（虢順民，江在森，張崇立.青藏高原東北緣晚第四紀塊體劃分與運動態(tài)勢研究［J］.地震地質(zhì)，2009，22（3）：219－231.）