時域最小二乘配置求定點形變特征曲線與異常識別的應用研究*

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

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

0 引言

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

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

1 最小二乘配置方法原理

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

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

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

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

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

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

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

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

2 定點形變曲線與異常識別應用實例

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

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

圖1 四川姑咱臺鉆孔應變東西分量日均值觀測曲線Fig.1 Daily mean observation cruve of E-W borehole strain at Guzan Station in Sichuan

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

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

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

圖2 姑咱臺鉆孔應變東西分量在不同相關范圍下日均值與擬合值曲線對比，絕對差值曲線及二倍均方差(參考區(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 姑咱臺鉆孔應變東西分量在不同參考區(qū)間內(nèi)的日均值與擬合值曲線對比、絕對值差異曲線及二倍均方差(相關范圍取缺省值)(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 趨勢加速或轉變明顯的實例

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

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

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

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

圖5 永勝臺垂直擺傾斜北南分量在不同參考區(qū)間內(nèi)日均值與擬合值曲線對比、絕對值差異曲線及二倍均方差(相關范圍取缺省值)(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 結論與討論

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

江在森，張希.2000.華北地區(qū)近期地殼水平運動與應力應變場特征［J］.地球物理學報，43(5):657 －665.

孔祥元，郭標明，劉宗泉，等.2001.大地測量學基礎［M］.武漢:武漢大學出版社.

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

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

武艷強，黃立人.2004.時間序列處理的新插值方法［J］.大地測量與地球動力學，24(4):43－47.

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

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

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

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