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

郭樹松 祝意青

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

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

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

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

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

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

1 前言

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

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

2 理論模型

2.1 基本原理

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

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

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

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

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

2.2 建立模型

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

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

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

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

其中

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

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

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

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

3 實(shí)例分析

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

我們選取該觀測網(wǎng)中的134、57和28號(hào)測點(diǎn)(圖2)，用本文方法研究其重力時(shí)序變化的穩(wěn)定性。1994—2009年該測段和測點(diǎn)的重力時(shí)序變化和時(shí)序系統(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)以上幾個(gè)例子較好地反映了該方法對重力時(shí)序系統(tǒng)是否穩(wěn)定的判斷的正確性，大量計(jì)算也表明，雖然各測點(diǎn)所受的干擾因素不同，構(gòu)造活動(dòng)對每一個(gè)測點(diǎn)的影響方式不同，分析某一區(qū)域內(nèi)各測點(diǎn)的時(shí)序穩(wěn)定性也不盡相同，但多數(shù)測點(diǎn)所反映的總體穩(wěn)定性趨勢是一致的。

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

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

1 祝意青，等.河西地區(qū)重力場及其動(dòng)態(tài)演化特征［J］.大地測量與地球動(dòng)力學(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 祝意青.青藏高原東北緣強(qiáng)震前兆特征研究［J］.國際地震動(dòng)態(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é)報(bào)，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ū)重力場動(dòng)態(tài)變化計(jì)算［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 夏德鈴.自動(dòng)控制理論［M］.北京：機(jī)械工業(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］.大地測量與地球動(dòng)力學(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è)科研專項(xiàng)(20090829)

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

P315.72+6

A