重力時序變化系統穩定性的研究*

郭樹松 祝意青

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

重力時序變化系統穩定性的研究*

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

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

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

重力時序變化;慢時變系統;差分方程;參數序列;河西重力觀測網

1 前言

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

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

2 理論模型

2.1 基本原理

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

對于單輸入單輸出系統的動態過程可描述為［5］：

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

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

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

2.2 建立模型

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

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

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

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

其中

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

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

2.3 系統失穩判定

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

3 實例分析

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

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

圖2 河西地區重力觀測網及構造分布Fig.2 Distribution of gravity survey network and tectonic diagram in Hexi area

4 結語

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

2)從數值分解的角度分析，我們通過建模將一列觀測值分解為n+1個序列，和觀測序列相比，引起每個參數變化的因素相對減少了，曲線中包含的成份簡單了。但是，數學模型只是分析系統的工具，并不代表形成系統的機理。因此對系統結構穩定性的判斷，仍然要從系統環境變化方面去推測，仍然要從分析參數曲線的正常基值線去理解。

圖3 134號測點1994—2009年的重力時序變化圖(a)及其系統矩陣根的模(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)及其系統矩陣根的模(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)及其系統矩陣根的模(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)本文只是初步研究，針對重力時序系統失穩要研究的問題很多，如失穩是否發生在最大的一個或兩個特征值序列上，最大根植序列是否有遷移，失穩的形態是持續型還是振蕩型，地震前兆異常的特征是什么…等，失穩并不都意味著地震前兆異常，這些問題還需要進行大量計算，結合構造活動尤其是震例在根軌跡圖上仔細分析，才能最后確定。

1 祝意青，等.河西地區重力場及其動態演化特征［J］.大地測量與地球動力學，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］.國際地震動態，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級地震關系研究［J］.地球物理學報，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 李輝，等.滇西地區重力場動態變化計算［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］.北京：機械工業出版社，2000.(Xia Deling.Automatic control theory［M］.Bejing：China Machine Press，2000)

6 謝新民，丁鋒.自適應控制系統［M］.北京：清華大學出版社，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級地震前區域重力場變化研究［J］.大地測量與地球動力學，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 王宏華.現代控制理論［M］.北京：電子工業出版社，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‘an 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

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

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

P315.72+6

A