賈 鵬，張 希，唐紅濤，李瑞莎

(中國(guó)地震局第二監(jiān)測(cè)中心，陜西西安710054)

0 引言

最小二乘配置是根據(jù)已知點(diǎn)信號(hào)、協(xié)方差及其與待估算點(diǎn)的協(xié)方差關(guān)系而獲得的待估算點(diǎn)的無偏最優(yōu)估計(jì)，它綜合了平差、推估和濾波。張希等(1998，1999a，b)，江在森和張希(2000)對(duì)該方法進(jìn)行過深入的研究和探討，將協(xié)方差經(jīng)驗(yàn)函數(shù)進(jìn)行了簡(jiǎn)化，并且將一維時(shí)間域內(nèi)的推估內(nèi)插進(jìn)行了驗(yàn)證，武艷強(qiáng)等(2007，2004)，馬超和單新建(2005)將其用于GPS連續(xù)站資料分析，驗(yàn)證了反映時(shí)序變化特征的可行性。

本文選取“十五”數(shù)字化臺(tái)站中兩個(gè)定點(diǎn)測(cè)項(xiàng)，分別為四川姑咱鉆孔應(yīng)變和云南永勝垂直擺傾斜，時(shí)間范圍為2008年1月1日至2011年12月31日(4年共1 461個(gè)日均值)。定點(diǎn)測(cè)項(xiàng)日均值很多(3年以上至少數(shù)千)，如果都作為已知觀測(cè)值，求逆矩陣相當(dāng)費(fèi)時(shí)，因此通過間隔一定天數(shù)取均值(本文均取5天)來進(jìn)行計(jì)算。并將張希和江在森(1999b)的水平面最小二乘配置擬合模型應(yīng)用于時(shí)間域內(nèi)插，并進(jìn)一步增加外推，應(yīng)用到定點(diǎn)形變觀測(cè)特征曲線尋找及其異常識(shí)別。

1 最小二乘配置方法原理

張希和江在森(1999b)，武艷強(qiáng)等(2007)，陶本藻等(1996)，孔祥元等(2001)都對(duì)最小二乘配置方法的原理進(jìn)行過相關(guān)的探討和論證，對(duì)所選取高斯型經(jīng)驗(yàn)協(xié)方差函數(shù)模型選擇簡(jiǎn)單進(jìn)行敘述。假設(shè)待內(nèi)插區(qū)域有m0個(gè)已知點(diǎn)觀測(cè)(或計(jì)算)值，設(shè)為L(zhǎng)=(g1，g2，…，gm0)T，其中每個(gè)點(diǎn)值gi的中誤差值為mgi(i=1，2，…，m0)。t為要濾波的已知點(diǎn)信號(hào)，n為觀測(cè)誤差向量，s是待估算點(diǎn)信號(hào)，t和n都是中心分布的。那么最小二乘配置的基本方程為

其中，Ctt為已知點(diǎn)信號(hào)t的先驗(yàn)自協(xié)方差矩陣;Cnn為觀測(cè)誤差向量的自協(xié)方差矩陣;Cst為待估算點(diǎn)信號(hào)與已測(cè)點(diǎn)信號(hào)t之間的協(xié)方差。經(jīng)推算的，的協(xié)方差分別為

將某一坐標(biāo)為(x，y)的待估算點(diǎn)的值表示為g，用c(a，b)表示變量a，b間協(xié)方差，則

上述各式中，Cst，Ctt均由高斯型經(jīng)驗(yàn)協(xié)方差函數(shù)f(d)=f(0)e－k2d2確定。其中，d為兩點(diǎn)間距離，k為待定參數(shù)(張希等，1998，1999b;江在森，張希，2000)。

假設(shè)Cnn為對(duì)角元素相等的對(duì)角陣，此對(duì)角元f(0)=fL(0)－fr(0)。為保證f(0)＞0，定義fr(0)=αfL(0)(0＜α≤0.2)，令α＞0表示為必須濾波。

確定參數(shù)是最小二乘配置實(shí)現(xiàn)的關(guān)鍵，但一般情況下很難得到可靠性較好的協(xié)方差圖形，故根據(jù)具體地區(qū)測(cè)點(diǎn)分布情況來確定參數(shù)k，即首先確定擬合量在整個(gè)區(qū)域的相關(guān)距離S(即超出這一距離，則點(diǎn)間協(xié)方差值接近于零)(張希，江在森，1999b)，則參數(shù)

設(shè) dij(i，j=1，2，…，m0)為任意兩點(diǎn)間距離，majx{mjindij}、Dmax=maix{majxdij}分別定義為最小相鄰點(diǎn)距、平均相鄰點(diǎn)距、最大相鄰點(diǎn)距、最大點(diǎn)距。相關(guān)距離S可取為

事實(shí)上，相關(guān)距離S越大，得到的擬合場(chǎng)會(huì)比較平緩，某些局部劇烈變化可能被平滑掉;而太小則不夠光滑，且對(duì)整個(gè)區(qū)域整體變化趨勢(shì)的反映不夠。k值可人為設(shè)定，2年為2，4年即為4;也可以范圍確定后，根據(jù)式(5)～(6)，程序自動(dòng)計(jì)算缺省值，本文兩個(gè)算例中k缺省計(jì)算值均為0.896(張希等，1998，1999b)。

2 定點(diǎn)形變曲線與異常識(shí)別應(yīng)用實(shí)例

2.1 參數(shù)的選取和確定

選取四川姑咱臺(tái)“十五”數(shù)字前兆觀測(cè)曲線(圖1)作為算例一。四川姑咱臺(tái)鉆孔應(yīng)變東西分量總體顯示準(zhǔn)周期性波動(dòng)特征，在2008年5月12日汶川M8.0地震前該測(cè)項(xiàng)產(chǎn)生了持續(xù)的壓性變化，震時(shí)至2008年8月下旬在下降趨勢(shì)中呈現(xiàn)顯著跳變，可能與攀枝花M6.1地震有一定關(guān)系，期間有人為停電的干擾，但總體呈現(xiàn)大震后非穩(wěn)定狀態(tài)，8月底至9月呈現(xiàn)波動(dòng)上升;2009年7～9月、2010年7～8月、2011年7～8月出現(xiàn)相似的波動(dòng)上升，觀測(cè)日志未說明原因、也無氣象數(shù)據(jù)，可能反映了季節(jié)因素。計(jì)算整個(gè)觀測(cè)時(shí)段相關(guān)范圍分別為缺省值、2和4的特征曲線，所得結(jié)果見圖2。

圖1 四川姑咱臺(tái)鉆孔應(yīng)變東西分量日均值觀測(cè)曲線Fig.1 Daily mean observation cruve of E-W borehole strain at Guzan Station in Sichuan

由圖2a可以看出，當(dāng)k值為缺省值0.896，可較好地模擬該測(cè)項(xiàng)4年來在總體下降的基礎(chǔ)上呈周期性波動(dòng)的特征，只有在剛起測(cè)不久，及2008年5～8月、2009年7～9月、2010年7～8月、2011年7～8月的這些時(shí)段差異較大，可能與汶川M8.0、攀枝花 M6.1地震、干擾或季節(jié)影響等有關(guān)。

當(dāng)選定相關(guān)范圍為2時(shí)，如圖2b所示，該擬合曲線總體也呈現(xiàn)下降且周期性的波動(dòng)，擬合曲線更加平滑，異常程度不顯著(超過二倍均方差的也只是2008年5～8月、2009年7～9月、2010年7～8月、2011年7～8月與地震、干擾或季節(jié)影響可能有關(guān)的時(shí)段);當(dāng)選定相關(guān)范圍為4時(shí)，則擬合曲線如圖2c所示，基本為一條下行的直線，波動(dòng)起伏不明顯，突出反映總體下降趨勢(shì)，看不出任何異常。由圖2可知，最小二乘配置對(duì)不同相關(guān)范圍(即不同頻段，式(5)中參數(shù)隨即改變)能反映觀測(cè)曲線的不同特征信息，與武艷強(qiáng)等(2004)研究結(jié)果類似，而從本例可以看出，選取相關(guān)范圍為缺省值0.896時(shí)更適合于定點(diǎn)測(cè)項(xiàng)的異常分析。如果認(rèn)為以往若干年(長(zhǎng)時(shí)間尺度)變化相對(duì)正?；蛴幸?guī)律，想借此判定近期變化形態(tài)是否與以往正常狀態(tài)一致、從而識(shí)別異常，可根據(jù)不同臺(tái)站資料具體情況選取一定時(shí)間段(小于參數(shù)k為缺省值的相關(guān)范圍，本文取小于0.896年的參考區(qū)間進(jìn)行外推。外推的參考區(qū)間分別選取為2008年1月1日至2011年6月30日(即外推至2011年下半年)和2008年1月1日至2011年9月30日(即外推至2011年第4季度)，外推結(jié)果見圖3。

從圖3a可以看出，2011年7～8月出現(xiàn)的波動(dòng)陡升更加顯著，異常差異明顯，對(duì)外推結(jié)果影響大，但總體年變周期趨勢(shì)良好;如果參考區(qū)間范圍改到2011年9月底結(jié)束(圖3b)，外推效果更佳，總體趨勢(shì)反應(yīng)好，說明在相關(guān)范圍內(nèi)進(jìn)行一定時(shí)間的外推對(duì)判定未知異常存在一定意義。因?yàn)楣迷叟_(tái)周期性較好、并無明顯異常，外推和實(shí)測(cè)值的結(jié)果符合周期性趨勢(shì)變化。

圖2 姑咱臺(tái)鉆孔應(yīng)變東西分量在不同相關(guān)范圍下日均值與擬合值曲線對(duì)比，絕對(duì)差值曲線及二倍均方差(參考區(qū)間為2008年1月1日至2011年12月31日)(a)k=0.896;(b)k=2;(c)k=4Fig.2 Comparison between daily mean and fitting curves of E-W borehole strain at Guzan Station，absolute difference curve and the double mean square error between daily mean and fitting values in different relevance range(time scale is from 2008－01－01 to 2011－12－31)

圖3 姑咱臺(tái)鉆孔應(yīng)變東西分量在不同參考區(qū)間內(nèi)的日均值與擬合值曲線對(duì)比、絕對(duì)值差異曲線及二倍均方差(相關(guān)范圍取缺省值)(a)2008－01－01至2011－06－30;(b)2008－01－01至2011－09－30Fig.3 Comparison between daily mean and fitting curves of E-W borehole strain at Guzan Station，absolute difference curve and the double mean square error between daily mean and fitting values in different time scale(k=0.896)

2.2 趨勢(shì)加速或轉(zhuǎn)變明顯的實(shí)例

永勝臺(tái)垂直擺北南分量從2007年至2010年6月總體呈現(xiàn)N向持續(xù)傾斜，汶川M8.0地震前出現(xiàn)過下行加速和轉(zhuǎn)折，震后上行加速明顯;攀枝花M6.1地震前有小幅折返，2010年6月前整體波動(dòng)并不劇烈;2010年6月出現(xiàn)明顯加速，波動(dòng)至2011年6月，目前南向轉(zhuǎn)折后幅度加大，有一些小的波動(dòng)或毛刺，可能觀測(cè)日志與調(diào)零、雷電干擾等有關(guān)，姚安地震后附近并無中強(qiáng)地震出現(xiàn)，但趨勢(shì)轉(zhuǎn)折變化特征非常明顯。故選云南永勝臺(tái)“十五”數(shù)字前兆觀測(cè)曲線作為算例二，如圖4所示。

計(jì)算整個(gè)觀測(cè)時(shí)段特征曲線(相關(guān)范圍為缺省計(jì)算值0.896)所得結(jié)果如圖5a所示。從圖中可以看到，擬合曲線較好地模擬4年來時(shí)緩、時(shí)陡的趨勢(shì)，擬合曲線總體可以分為兩個(gè)階段:2010年6月前上行速度緩慢;2010年6月后速度加快，轉(zhuǎn)折幅度大。差異比較大的時(shí)段出現(xiàn)在2008年的4～5月，2010年的7～8月及2011年年末，可能受到汶川M8.0地震、雷電和調(diào)零的影響，但可以反映出幾個(gè)時(shí)段內(nèi)的曲線異常變化。

圖4 云南永勝臺(tái)垂直擺傾斜北南分量日均值觀測(cè)曲線Fig.4 Daily mean observational records of N-S vertical pendulum tiltmeter at Yongsheng Station in Yunnan

進(jìn)而通過外推求該測(cè)項(xiàng)是否存在破年變異常，將參考區(qū)間改為2008年至2011年6月底，外推2011年7～12月，如圖5b所示。2010年6月的線性加速趨勢(shì)更加明顯，而2011年下半年，尤其9月份以后，測(cè)項(xiàng)轉(zhuǎn)折下降趨勢(shì)與擬合上行加速截然不同;若參考區(qū)間改到2011年9月底結(jié)束，外推至2011年10～12月，如圖5c所示，2011年7～8月趨勢(shì)并未改變，而外推的預(yù)測(cè)結(jié)果與實(shí)際測(cè)項(xiàng)9月轉(zhuǎn)折下行依然相反，說明存在破年變異常。綜合分析認(rèn)為通過最小二乘配置外推能夠檢驗(yàn)偏離已有數(shù)學(xué)模型所能描述的部分信息，對(duì)輔助判定未知異常有一定意義，但不能作為識(shí)別異常的確切指標(biāo)。

圖5 永勝臺(tái)垂直擺傾斜北南分量在不同參考區(qū)間內(nèi)日均值與擬合值曲線對(duì)比、絕對(duì)值差異曲線及二倍均方差(相關(guān)范圍取缺省值)(a)2008－01－01至2011－12－31;(b)2008－01－01至2011－06－30;(c)2008－01－01至2011－09－30Fig.5 Comparison between daily mean and fitting curves of N-S vertical pendulum tiltmeter at Yongsheng Station，absolute difference curve and the double mean square error between daily mean and fitting values in different time scale(k=0.896)

3 結(jié)論與討論

利用最小二乘配置進(jìn)行定點(diǎn)形變曲線特征模擬和尋找異常特征，可以較好地處理連續(xù)變化的點(diǎn)，能夠考慮到內(nèi)插區(qū)域內(nèi)所有已知點(diǎn)的相關(guān)性，反映其隨時(shí)間變化的趨勢(shì)性。外推時(shí)根據(jù)時(shí)間的相關(guān)范圍內(nèi)取一定時(shí)段外推，能夠輔助識(shí)別出定點(diǎn)形變曲線的短、中期趨勢(shì)變化和異常走勢(shì)，對(duì)于輔助識(shí)別定點(diǎn)形變曲線異常有一定的意義。最小二乘配置方法僅是一種數(shù)學(xué)處理方法，通過濾波和擬合可以處理接近觀測(cè)誤差噪聲相對(duì)中、高頻的部分，可以較好地反應(yīng)曲線周期，但并不能代替其他如小波分析等濾波數(shù)學(xué)工具，外推對(duì)于異常的最終判定還需其它物理手段予以落實(shí)。

江在森，張希.2000.華北地區(qū)近期地殼水平運(yùn)動(dòng)與應(yīng)力應(yīng)變場(chǎng)特征［J］.地球物理學(xué)報(bào)，43(5):657 －665.

孔祥元，郭標(biāo)明，劉宗泉，等.2001.大地測(cè)量學(xué)基礎(chǔ)［M］.武漢:武漢大學(xué)出版社.

馬超，單新建.2005.昆侖山口西MS8.0地震INSAR斜距向同震位錯(cuò)分解［J］.地震研究，28(3):245 －247.

陶本藻，周勇前，高士純，等.1996.測(cè)量平差基礎(chǔ)［M］.武漢:測(cè)繪出版社.

武艷強(qiáng)，黃立人.2004.時(shí)間序列處理的新插值方法［J］.大地測(cè)量與地球動(dòng)力學(xué)，24(4):43－47.

武艷強(qiáng)，江在森，楊國(guó)華.2007.最小二乘配置方法在提取GPS時(shí)間序列信息中的應(yīng)用［J］.國(guó)際地震動(dòng)態(tài)，(7):99－103.

張希，江在森，張四新.1998.借助最小二乘配置整體解算地殼視應(yīng)變場(chǎng)［J］.地殼形變與地震，18(2):57－62.

張希，江在森.1999a.對(duì)華北GPS監(jiān)測(cè)區(qū)近期地殼應(yīng)變連續(xù)分布的估計(jì)［J］.地震學(xué)刊，75(2):47－52.

張希，江在森.1999b.用最小二乘配置獲得地形變應(yīng)變場(chǎng)動(dòng)態(tài)圖像的幾個(gè)問題研究［J］.地殼形變與地震，19(3):32－39.