郭樹松 祝意青

(1)中國地震局第二監(jiān)測中心，西安 710054 2)中國科學(xué)院測量與地球物理研究所動力大地測量學(xué)重點實驗室，武漢430077)

重力時序變化系統(tǒng)穩(wěn)定性的研究*

郭樹松1)祝意青1，2)

(1)中國地震局第二監(jiān)測中心，西安 710054 2)中國科學(xué)院測量與地球物理研究所動力大地測量學(xué)重點實驗室，武漢430077)

把重力時序系統(tǒng)視為慢時變系統(tǒng)，根據(jù)線性系統(tǒng)分析的理論方法，用差分方程描述該系統(tǒng)的等價參數(shù)模型，通過判斷參數(shù)序列的穩(wěn)定性來探索前兆異常。對河西流動重力觀測網(wǎng)1994—2009年3個測點的觀測數(shù)據(jù)的計算分析驗證了該方法的有效性，但失穩(wěn)與地震前兆的關(guān)系還需進一步研究。

重力時序變化;慢時變系統(tǒng);差分方程;參數(shù)序列;河西重力觀測網(wǎng)

1 前言

地震前兆系統(tǒng)分析是當(dāng)今科學(xué)研究領(lǐng)域的一大難題，對該系統(tǒng)而言，其結(jié)構(gòu)機制與內(nèi)部機理均不甚明了，無法通過動力系統(tǒng)的模型方法來探索其復(fù)雜特性和過程。研究表明，地震的孕育和發(fā)生過程伴隨著由于斷層構(gòu)造活動和局部地殼應(yīng)力集中而引起地殼形變(高程)和密度(質(zhì)量)的變化，從而引起包括重力場在內(nèi)的多種地球物理場的變化。流動重力測量反映了構(gòu)造活動區(qū)重力場隨時間的非潮汐變化，地殼內(nèi)部的密度異常、質(zhì)量遷移和地震的形成過程等都可在流動重力復(fù)測結(jié)果中反映出來［1，2］。因此，研究區(qū)域重力場的時空動態(tài)演化特征是探索地震前兆異常的一條重要途徑。

目前對流動重力觀測資料的研究方法可以分為兩類，一類通過網(wǎng)格化對數(shù)據(jù)進行內(nèi)插濾波，得出二維的重力場連續(xù)變化圖像，從整體上分析測網(wǎng)內(nèi)應(yīng)力-應(yīng)變場微動態(tài)活動;另一類分析特殊區(qū)域(一般為重力變化高梯度帶和正、負異常變化的過渡區(qū)域)內(nèi)測點或剖面重力的時序變化，時序變化能較好地突出異常動態(tài)，從而探索斷裂的構(gòu)造活動。在研究中一般把兩者綜合起來進行地震趨勢預(yù)報［3，4］。但是，對許多地震來說不能直接得到重力變化的范圍，地震也不是發(fā)生在重力變化的最大點上，短期重力變化的大小也可能包含著測量粗差及淺表局部干擾。針對這樣的問題，本文探索一種判斷重力時序變化穩(wěn)定性的方法，把長期的時序變化看作一個慢時變系統(tǒng)，通過分析時序系統(tǒng)結(jié)構(gòu)變化和穩(wěn)定性的過程探索前兆異常。

2 理論模型

2.1 基本原理

流動重力觀測網(wǎng)所構(gòu)成的蘊震系統(tǒng)是一個非線性、非平衡的慢時變系統(tǒng)，其演變過程是一系列的離散數(shù)據(jù)。對于離散系統(tǒng)一般用差分方程建立參數(shù)模型，系統(tǒng)辨識就是確定這些參數(shù)。模型參數(shù)是對系統(tǒng)行為結(jié)構(gòu)的描述，定常系統(tǒng)的結(jié)構(gòu)不變，參數(shù)亦不變，時變系統(tǒng)的結(jié)構(gòu)在改變，相應(yīng)的參數(shù)亦改變。因此時變參數(shù)序列就反映了系統(tǒng)結(jié)構(gòu)的改變，這是我們探索地震前兆的基礎(chǔ)。

對于單輸入單輸出系統(tǒng)的動態(tài)過程可描述為［5］：

其中u(k)、y(k)分別為系統(tǒng)的輸入、輸出，a1(k)、…、an(k)為模型參數(shù)，n為系統(tǒng)的階數(shù)，時變系統(tǒng)中參數(shù)是時間的函數(shù)。前兆系統(tǒng)的建模就是確定方程的階數(shù)n，求出n個時變參數(shù)函數(shù)，即n個參數(shù)的序列值。

由線性理論分析得知，方程(1)的特征方程為

這是復(fù)平面Z上的方程，其根稱為系統(tǒng)的閉環(huán)極點。閉環(huán)系統(tǒng)的動態(tài)性能與閉環(huán)極點在Z平面上的位置密切相關(guān)。參數(shù)變化時，特征方程的根在Z平面上運動的軌跡稱為根軌跡。在自控工程設(shè)計中，正是利用閉環(huán)極點位置的改變來改善系統(tǒng)穩(wěn)定性及其動態(tài)性能。而在前兆系統(tǒng)分析中，正好相反，是利用根軌跡來分析系統(tǒng)的失穩(wěn)過程與地震前兆異常的關(guān)系。

2.2 建立模型

如圖1所示，我們可以把輸入函數(shù)g(t)合理地劃分為若干個時段，每個時段的輸入函數(shù)為u(k)。假定每個時段的輸入函數(shù)可以用一個階躍函數(shù)近似表示，例如u(k)=b(k)u(1)，其中u(1)為單位階躍函數(shù)，b(k)為該時段的躍度，是一個未知值。故式(1)可表示為

從式(3)可以看出，時變系統(tǒng)的輸出中，既有輸入的因素，又有參數(shù)變化的因素;也就是說，系統(tǒng)輸出中也包含了輸入的信息。因此，在求時變參數(shù)時，可同時求出該時段輸入函數(shù)的躍度b(k)值。因此，式(3)就是我們采用的模型方程。

圖1 輸入函數(shù)的曲線模擬Fig.1 Simulation curve of input function

實際上，時變參數(shù)估計可采用實時遞推估計法［6］的公式計算。該方法的遞推公式為：

其中

式中，θT(N)是由輸入、輸出的模型參數(shù)組成的矩陣，對于遞推初值(0)為任意值;P(0)=α2I，α為無限大，I為單位矩陣。λ為遺忘因子，取值在0到1之間，它對老數(shù)據(jù)不斷進行截尾，從而使實時算法一直保持著對參數(shù)估計的校正能力。

這里需要注意的是確定參數(shù)的階數(shù)n是建模成敗的關(guān)鍵，遺忘因子λ取值的大小直接影響到對參數(shù)的跟蹤能力。對每一個定點觀測系統(tǒng)，二者都要經(jīng)過反復(fù)試算比較才能確定最佳取值。

2.3 系統(tǒng)失穩(wěn)判定

系統(tǒng)失穩(wěn)判定是分析前兆異常的重點，首先要解決特征方程求根問題。在線性系統(tǒng)分析中，一般是將差分方程描述通過變量代換化為狀態(tài)變量描述，與狀態(tài)變量相乘的矩陣A稱為系統(tǒng)矩陣。對單輸入、單輸出系統(tǒng)而言，矩陣A的特征值就是特征方程(2)的根。由此推導(dǎo)出矩陣A的組成為：

3 實例分析

北祁連河西流動觀測網(wǎng)位于青藏高原東北緣(圖2)。測網(wǎng)所在區(qū)域是中國大陸地殼運動最強烈、地震活動頻度最高、強度最大的地區(qū)之一。從1989年開始中國地震局第二地形變監(jiān)測中心在該地區(qū)初步建立了流動重力監(jiān)測網(wǎng)，后經(jīng)多次擴建和加密，形成目前的北祁連河西流動觀測網(wǎng)。1989年開始在該地區(qū)每年進行一次重力測量(1993年停測1期，2009年開始每年觀測2期)，至今已有21期測量成果。各期觀測精度都在12×10-8ms-2以上，觀測精度較高，觀測資料可靠。

我們選取該觀測網(wǎng)中的134、57和28號測點(圖2)，用本文方法研究其重力時序變化的穩(wěn)定性。1994—2009年該測段和測點的重力時序變化和時序系統(tǒng)矩陣根的模如圖3～5所示。其中參數(shù)模型的階數(shù)n=3，遺忘因子λ=0.95。

圖2 河西地區(qū)重力觀測網(wǎng)及構(gòu)造分布Fig.2 Distribution of gravity survey network and tectonic diagram in Hexi area

4 結(jié)語

1)以上幾個例子較好地反映了該方法對重力時序系統(tǒng)是否穩(wěn)定的判斷的正確性，大量計算也表明，雖然各測點所受的干擾因素不同，構(gòu)造活動對每一個測點的影響方式不同，分析某一區(qū)域內(nèi)各測點的時序穩(wěn)定性也不盡相同，但多數(shù)測點所反映的總體穩(wěn)定性趨勢是一致的。

2)從數(shù)值分解的角度分析，我們通過建模將一列觀測值分解為n+1個序列，和觀測序列相比，引起每個參數(shù)變化的因素相對減少了，曲線中包含的成份簡單了。但是，數(shù)學(xué)模型只是分析系統(tǒng)的工具，并不代表形成系統(tǒng)的機理。因此對系統(tǒng)結(jié)構(gòu)穩(wěn)定性的判斷，仍然要從系統(tǒng)環(huán)境變化方面去推測，仍然要從分析參數(shù)曲線的正?；稻€去理解。

圖3 134號測點1994—2009年的重力時序變化圖(a)及其系統(tǒng)矩陣根的模(b)Fig.3 Gravity time-variation at 134 observation point from 1994 to 2009(a)and the modules of characteristic root of its coefficient matrix(b)

圖4 57號測點1994—2009年的重力時序變化圖(a)及其系統(tǒng)矩陣根的模(b)Fig.4 Gravity time-variation at 57 observation point from 1994 to 2009(a)and the modules of characteristic root of its coefficient matrix(b)

圖5 28號測點1994—2009年的的重力時序變化圖(a)及其系統(tǒng)矩陣根的模(b)Fig.5 Gravity time-variation at 28 observation point from 1994 to 2009(a)and the modules of characteristic root of its coefficient matrix(b)

3)本文只是初步研究，針對重力時序系統(tǒng)失穩(wěn)要研究的問題很多，如失穩(wěn)是否發(fā)生在最大的一個或兩個特征值序列上，最大根植序列是否有遷移，失穩(wěn)的形態(tài)是持續(xù)型還是振蕩型，地震前兆異常的特征是什么…等，失穩(wěn)并不都意味著地震前兆異常，這些問題還需要進行大量計算，結(jié)合構(gòu)造活動尤其是震例在根軌跡圖上仔細分析，才能最后確定。

1 祝意青，等.河西地區(qū)重力場及其動態(tài)演化特征［J］.大地測量與地球動力學(xué)，2003，(4)：44-48.(Zhu Yiqing，et al.Study on gravity field and its dynamic evolutional characteristics in Hexi area［J］.Journal of Geodesy and Geodynamics，2003，(4)：44-48)

2 祝意青.青藏高原東北緣強震前兆特征研究［J］.國際地震動態(tài)，2007，(5)：16-21.(Zhu Yiqing.Precursory characteristics of stronger earthquakes innortheastern edge of Qinghai-Tibet plateau by using mobile gravity technique［J］.Recent Developments in World Seismology，2007，(5)：16-21)

3 祝意青，等.龍門山斷裂帶重力變化與汶川8.0級地震關(guān)系研究［J］.地球物理學(xué)報，2009，52(10)：2 538-2 546. (Zhu Yiqing，et al.Relations between gravity variation of Longmenshan fault zone and Wenchuan Ms8.0 earthquake［J］.Chinese Journal of Geophysics，2009，52(10)：2 538 -2 546)

4 李輝，等.滇西地區(qū)重力場動態(tài)變化計算［J］.地殼形變與地震，2000，20(1)：60-66.(Li Hui，et al.Computation of dynamic gravity changes in Westen Yunnan［J］.Crustal Deformation and Earthquake，2000，20(1)：60-66)

5 夏德鈴.自動控制理論［M］.北京：機械工業(yè)出版社，2000.(Xia Deling.Automatic control theory［M］.Bejing：China Machine Press，2000)

6 謝新民，丁鋒.自適應(yīng)控制系統(tǒng)［M］.北京：清華大學(xué)出版社，2002.(Xie Xinmin and Ding Feng.Adaptive control system［M］.Bejing：Qinghua University Press，2002)

7 祝意青，等.景泰5.9級地震前后的重力變化研究［J］.中國地震，2001，17(4)：356-363.(Zhu Yiqing，et al.Research on gravity variation before and after Jingtai Ms5.9 earthquake［J］.Earthquake Researth，2001，17(4)：356-363)

8 祝意青，等.民樂6.1、岷縣5.2級地震前區(qū)域重力場變化研究［J］.大地測量與地球動力學(xué)，2005，(1)：24-29.(Zhu Yiqing，et al.Research on the variation of gravity field before Minle Ms6.1 and Minxian Ms5.2 earthquakes［J］.Journal of Geodesy and Geodynamics，2005，(1)：24-29)

9 王宏華.現(xiàn)代控制理論［M］.北京：電子工業(yè)出版社，2006.(Wang Honghua.Mordern control theory［M］.Beijng：Publishing House of Electronics Industry，2006)

RESEARCH ON STABLILITY OF GRAVITY TIME-VARYING SYSTEM

Guo Shusong1)and Zhu Yiqing1，2)

(1)Second Crust Monitoring and Application Center，CEA，Xi‘a(chǎn)n 710054 2)Institute of Geodesy and Geophysics，CAS，Wuhan430077)

We regard gravity time sequence system as slow time-varying system and use difference equations for describing its equivalent parameter model based on the theory of linear systems analysis.By judging the stability of parametric sequence，we are going to explore the precursors.On the basis of gravity time-varying series from 1994 to 2009 of three observation points in Hexi gravity survey network，we prove that this method is effective，but more research are required to confirm the relations between unstability and seismic precursor anomalies.

gravity time-variation;slow time-varying system;difference equations;parametric sequence;Hexi gravity survey network

1671-5942(2011)05-0061-05

2011-03-09

國家自然科學(xué)基金(40874035);地震行業(yè)科研專項(20090829)

郭樹松，男，1980年生，碩士，主要從事重力數(shù)據(jù)處理及應(yīng)用研究.E-mail：98212541@sina.com

P315.72+6

A