趙麗華，楊元喜，王慶良

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

1 引 言

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

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

2 相對(duì)穩(wěn)定點(diǎn)組的確定

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

（1）將各點(diǎn)速率觀測(cè)值分別取中位數(shù)med（L），求各速率值與相應(yīng)中位數(shù)之差，即初始?xì)埐?/p>

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

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

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

式中，k0可取1.0或1.5。

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

式中，k＝1.5～2。

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

3 塊體運(yùn)動(dòng)模型的精化及相鄰塊體邊界的確定

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

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

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

式中，S′為推估點(diǎn)上信號(hào)，S′～N（0，ΣS′），且基于最小二乘原則可得

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

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

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

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

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

（4）假設(shè)檢驗(yàn)

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

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

迭代的過(guò)程即是對(duì)形變模型精化的過(guò)程，同時(shí)也確定了斷層邊界的分布。

4 計(jì)算與分析

4.1 算例1

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

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

圖1 模擬觀測(cè)點(diǎn)及斷裂分布Fig.1 The distribution of simulated measurement points and fault

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

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

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

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

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

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

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

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

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

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

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

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

圖4 各方案檢核點(diǎn)不符值Fig.4 Residuals of check points

4.2 算例2

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

圖5 各點(diǎn)位置分布及速率Fig.5 Distribution and velocity of observation points

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

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

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

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

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

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

圖6 算例2中不同方案模擬的水平運(yùn)動(dòng)與實(shí)測(cè)值比較Fig.6 Comparison between simulated velocity and observed velocity

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

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

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

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

（4）由于篇幅限制，沒(méi)有一一列出各檢核點(diǎn)不符值（擬合模型估計(jì)的檢核點(diǎn)速率與該點(diǎn)實(shí)測(cè)速率之差），但由檢核點(diǎn)不符值計(jì)算的RMS也反映出分區(qū)域抗差擬合推估法估計(jì)的各點(diǎn)水平運(yùn)動(dòng)速率更接近真實(shí)值（表4）。

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

表4 不同方案檢核點(diǎn)RMSTab.4 RMS of check points mm

5 結(jié) 論

若形變場(chǎng)中有隱伏斷層或存在較顯著形變差異時(shí)，有時(shí)無(wú)法確切地知道差異運(yùn)動(dòng)的分界帶。這時(shí)若采用一般擬合推估方法建立形變場(chǎng)，在斷層附近會(huì)產(chǎn)生較大偏差，無(wú)法可靠地反映斷層分布情況。本文首先通過(guò)地質(zhì)資料進(jìn)行簡(jiǎn)單判斷，對(duì)可能存在斷層的區(qū)域進(jìn)行概略劃分，然后通過(guò)高崩潰污染率抗差估計(jì)法挑選各區(qū)域內(nèi)相對(duì)穩(wěn)定的點(diǎn)，建立各區(qū)域形變場(chǎng)置信區(qū)間。對(duì)分布在斷層附近分界不明點(diǎn)的歸屬進(jìn)行判斷，確定各塊體差異運(yùn)動(dòng)的邊界。邊界明確后，分別建立各區(qū)域能夠描述本塊體運(yùn)動(dòng)趨勢(shì)的擬合推估模型，計(jì)算區(qū)域形變場(chǎng)。試驗(yàn)算例表明，該方法在發(fā)現(xiàn)隱伏斷層并確定其分布邊界的同時(shí)，建立的形變場(chǎng)較好地表現(xiàn)出不同運(yùn)動(dòng)塊體邊界上的差異運(yùn)動(dòng)，具有較高推估精度。當(dāng)研究區(qū)域觀測(cè)點(diǎn)數(shù)量相對(duì)較少時(shí)，利用此方法解得的形變參數(shù)可以推算任意未測(cè)點(diǎn)的形變量，達(dá)到分析形變和應(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é)果揭示中國(guó)大陸水平應(yīng)變場(chǎng)與構(gòu)造變形［J］.地球物理學(xué)報(bào)，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.（楊元喜.抗差估計(jì)理論及其應(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ì)［J］.測(cè)繪學(xué)報(bào)，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.（歐吉坤.粗差的擬準(zhǔn)檢定法（QUAD法）［J］.測(cè)繪學(xué)報(bào)，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］.測(cè)繪學(xué)報(bào)，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)計(jì)［M］.西安：西安交通大學(xué)出版社，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.（虢順民，江在森，張崇立.青藏高原東北緣晚第四紀(jì)塊體劃分與運(yùn)動(dòng)態(tài)勢(shì)研究［J］.地震地質(zhì)，2009，22（3）：219－231.）