2007年5月—2013年12月成都、 西昌和重慶臺地磁轉(zhuǎn)換函數(shù)的時間變化

龔紹京 劉雙慶 張明東

（中國天津300201天津市地震局）

?

2007年5月—2013年12月成都、 西昌和重慶臺地磁轉(zhuǎn)換函數(shù)的時間變化

（中國天津300201天津市地震局）

本文較詳細地討論了地磁轉(zhuǎn)換函數(shù)的計算方法， 介紹了復轉(zhuǎn)換函數(shù)實部和虛部的誤差公式．利用成都、 西昌和重慶3個地磁臺2007年5月—2013年12月的地磁數(shù)字化資料， 計算出轉(zhuǎn)換函數(shù)A，B及帕金森矢量．結(jié)果表明， 在2008年汶川MS8.0地震前后成都臺存在可察覺的小鼓包， 而西昌、 重慶兩臺則沒有可察覺的異常．另外與20世紀80年代的結(jié)果比較， 汶川地震前后成都臺的b值和帕金森矢量的長度及方向都有明顯的長期變化.

地磁轉(zhuǎn)換函數(shù) 復數(shù)最小二乘法 多元回歸分析 穩(wěn)健估計 帕金森矢量 地震異常

引言

1940年賈普曼最早提出地殼的不規(guī)則性可能會影響地磁場尤其是較短周期的瞬變場（Chapman, Bartels, 1940）, 這一論點已被20世紀50—60年代發(fā)現(xiàn)的大量事實所證實． Parkinson （1959）通過研究快速地磁變化矢量的方向， 發(fā)現(xiàn)地磁變化矢量有限定在一個平面上的趨勢, 此平面被稱為“地磁變化矢量優(yōu)勢面”. 在此基礎上提出了帕金森矢量的概念， 推導出著名的帕金森矢量系數(shù)經(jīng)驗公式． 大量研究表明， 帕金森矢量有3種反映地下電性結(jié)構(gòu)的效應： 內(nèi)陸異常、 海岸效應和電流通道（Parkinson, 1959, 1962； Untiedt, 1970； 龔紹京, 1987）．

此后引入具有嚴格周期概念的轉(zhuǎn)換函數(shù)， 它表達輸出與輸入之間的函數(shù)關系， 用以描述復頻率域內(nèi)線性系統(tǒng)的動態(tài)特性． 在地磁領域， 輸入是外源施感場， 輸出是內(nèi)源感應場． 由于分解地磁內(nèi)外源場很困難且需較多臺站資料， Schmucker（1970）引入正常場和異常場的概念， 假設外源場準均勻， 經(jīng)近似簡化處理推導出較簡便且實用的表達形式， 使單臺和雙臺地磁復轉(zhuǎn)換函數(shù)的實際應用成為可能并取得許多成果． 按Schmucker的概念， 當?shù)叵码娦越Y(jié)構(gòu)為水平分層且橫向比較均勻時， 其轉(zhuǎn)換函數(shù)值很??； 垂直場轉(zhuǎn)換函數(shù)及由它們組成的帕金森矢量的長度反映地下電性的橫向不均勻程度， 矢量指向?qū)щ娐矢叩囊环剑?/p>

20世紀70年代以來， 國外許多研究人員利用研究地下電性結(jié)構(gòu)的參數(shù)來探尋地震引起的地下電性變化， 出現(xiàn)明顯異常的有關東地震（Yanagihara， 1972）、 塔什干地震（Miyakoshi， 1975）、 錫特卡地震（Rikitake， 1979）和卡萊爾地震（Beamish， 1982）． 從20世紀80年代起， 作者利用磁變儀資料， 采用量圖和磁照圖數(shù)字化方法， 計算帕金森矢量系數(shù)a和b（Gong， 1985； 龔紹京等， 1989； 林美， 龔紹京， 1991）、 單臺垂直場轉(zhuǎn)換函數(shù)A和B（龔紹京等， 1991）以及水平場臺際轉(zhuǎn)換函數(shù)C、G、E、F（龔紹京等， 1997， 2001）， 發(fā)現(xiàn)在唐山、 松潘、 菏澤、 花蓮地震前后均存在或大或小、 或長或短的異常． 特別對唐山地震， 作者研究了上述所有參量， 發(fā)現(xiàn)a（Gong， 1985）和Au、Cu、Cv、Fv在唐山地震前后有持續(xù)約5—7年的異常， 震后一段時間才恢復（龔紹京等， 1997）； 同時水平場轉(zhuǎn)換函數(shù)C和F在震前的1976年2—4月出現(xiàn)十分明顯的短期前兆， 這是1972—1997年26年間出現(xiàn)的唯一一次短期前兆（龔紹京等， 2001）． 陳伯舫（1998）也研究了關島等地震， 同樣發(fā)現(xiàn)了這些參數(shù)的異常．

對復轉(zhuǎn)換函數(shù)的計算， Beamish （1982）采用功率譜互譜方法， Everett和Hyndman （1967）以及作者采用復數(shù)最小二乘法（龔紹京等， 1991， 1997， 2001）． 國內(nèi)一些研究人員對轉(zhuǎn)換函數(shù)也提出了各自的計算方法和公式， 對其中某些問題龔紹京等（2012）已經(jīng)作過分析. 本文擬對新發(fā)現(xiàn)的問題進行分析探討．

2007—2008年間， 我國的大部分地磁相對觀測已改用磁通門磁力儀， 實現(xiàn)了數(shù)字化記錄． 隨著計算科學的發(fā)展出現(xiàn)了許多新的計算語言， 為適應這些新的變化， 作者將過去龔紹京編的Fortran語言程序改編成Matlab語言程序*在改編程序的課題中， Matlab程序的大部分內(nèi)容如提取數(shù)據(jù)、 顯示圖形、 消除干擾、 顯示菱形截取數(shù)據(jù)等由劉雙慶編寫. 復數(shù)最小二乘法求轉(zhuǎn)換函數(shù)及實、 虛部誤差的計算， 穩(wěn)健估計的應用以及求帕金森矢量部分由龔紹京編寫.， 并用新程序計算了成都、 西昌、 重慶3個地磁臺的轉(zhuǎn)換函數(shù)A，B和帕金森矢量， 擬探討汶川MS8.0地震前后是否存在異?，F(xiàn)象．

1 轉(zhuǎn)換函數(shù)計算方法， 公式及討論

1.1 帕金森矢量系數(shù)

帕金森矢量系數(shù)經(jīng)驗公式為

ΔZ=a·ΔH+b·ΔD，

（1）

1.2 復轉(zhuǎn)換函數(shù)計算方法及實部和虛部的誤差公式

經(jīng)簡化處理后的地磁復轉(zhuǎn)換函數(shù)表達式為

Z=A·H+B·D+εz，

（2）

（3）

式中：A，B為單臺垂直場轉(zhuǎn)換函數(shù)；C，G，E，F(xiàn)為水平場臺際轉(zhuǎn)換函數(shù)；ε為殘差. 這些參數(shù)皆為復數(shù)且是頻率的函數(shù)． 下標“a”代表異常場， “n”代表正常場， “k”代表異常區(qū)臺站的磁場． 對于式（2）、 （3）的常規(guī)計算方法請參閱相關文獻（Everett, Hyndman, 1967; 龔紹京等, 1991, 1997, 2001）， 但其中未給出A，B實部和虛部的誤差公式.

用復數(shù)最小二乘法推導的計算復轉(zhuǎn)換函數(shù)的公式為

（4）

其中

（5）

A、B實部和虛部的誤差公式由本文第一作者采用單位矢量法和拉格朗日乘數(shù)法則推導， 并曾寄給陳伯舫博士驗證．A實部與虛部的誤差公式為

（6）

同理也可推導出B實部與虛部誤差公式為

（7）

1.3 轉(zhuǎn)換函數(shù)計算方法討論

由于Matlab語言為諸多數(shù)學問題提供了成熟的子函數(shù)程序， 有人提出了與復數(shù)最小二乘不同的算法， 即將等式兩邊除同一變量， 將聯(lián)立方程組中的復系數(shù)和復變量分解成實部和虛部， 再采用多元回歸分析處理方程數(shù)倍增的聯(lián)立方程組． 例如對式（2）進行如下變換得到式（8）.

令Z/D=Tr+iTi，S=H/D=Sr+iSi， 則T=A·S+B+ε. 由此可以得到

（8）

式中所有元素都為實數(shù)， 有4個待定參數(shù)， 變成2m個方程． 調(diào)用Matlab中的多元回歸（regress）子函數(shù)可求出4個系數(shù)值．

然而， 依據(jù)式（8）調(diào)用regress子函數(shù)得到的結(jié)果與依據(jù)復數(shù)最小二乘法（式（4））得到的結(jié)果有一定差別， 有時差別還相當大， 如圖1所示． 顯然圖中“方法1”（姑且稱之為‘多元回歸分析法’）的離散度較采用復數(shù)最小二乘的“龔算法1”大， 效果差一些． 起初找不到導致這種差別的原因， 嘗試了多方面的分析對比研究． 當確認式（4）推導無誤后， 曾懷疑是否多元回歸子函數(shù)計算的結(jié)果與復數(shù)最小二乘的結(jié)果不同， 但最終找出了出現(xiàn)差別的原因： 變換使式（8）與式（2）的殘差表達式不同， 因而對殘差平方和偏微商的結(jié)果也不相同． 為此用數(shù)字試驗的辦法進行了驗證， 即用復數(shù)最小二乘法的式（4）去計算經(jīng)變換（除D）后的數(shù)據(jù)， 用T代替式（5）中的Z， 用S代替H，D則為1． 所得到的結(jié)果與調(diào)用多元回歸子函數(shù)的結(jié)果完全一樣．

圖1 成都臺2007年9月22—28日（a）與2009年7月11—21日（b）的周期響應曲線

同時對最簡單的情形進行了研究. 設z=a·x+b·y， 式中的元素均為實數(shù)． 用y（或x）去除等式的兩邊． 如果是確定性問題， 即如果是解方程， 那么， 未進行變換得到的a、b值與進行變換后得到的a、b值相同． 但如果用最小二乘法進行“估算”， 由于存在殘差， 數(shù)字試驗表明， 進行變換與未變換的兩種估計值是有差別的． 估計值的殘差越小， 越接近確定性問題的解， 兩種算法的差別也越??； 殘差越大， 兩種方法的差別也越大． 這從圖1中也可得到證實． 圖1中的周期響應曲線僅用了式（4）和（8）的計算結(jié)果， 未引入穩(wěn)健估計．

1.4 穩(wěn)健估計方法效果評價

穩(wěn)健統(tǒng)計方法是在加權統(tǒng)計方法基礎上發(fā)展起來的， 其中心思想是選擇某種統(tǒng)計量， 使理想假設不受某些偏離數(shù)據(jù)的影響， 并引入“極端點”（outlier points）的概念， 以便自動識別偏離的數(shù)據(jù)和自動修正所加之權重， 通過反復迭代使估計值趨于穩(wěn)定． 關于穩(wěn)?。≧obust）估計的原理和做法詳見Egbert和Booker（1986）以及龔紹京等（1997， 2001）文章．

具體進行穩(wěn)健估算時， 需設定兩個控制參量， 即加權時用到的損失函數(shù)以及迭代終止的判據(jù)． 作者共設計了3種穩(wěn)健估計， 其中一種調(diào)用了Matlab中的子函數(shù). 圖2給出了3種穩(wěn)健估計的周期響應曲線. 可以看出， 用方法1（式（8））得到的實線與本文用虛線表示的3種穩(wěn)健估計曲線存在差異， 而3種虛線的差別在圖中幾乎看不出來， 僅數(shù)值有小的差別． 這也進一步驗證了作者過去和現(xiàn)在所采用方法的可靠性．

圖2 3種穩(wěn)健估計方法的周期響應曲線

Fig.2 Period respond curves of three robust estimations. The solid line denotes the result by method 1, and the dotted line denote three kinds of robust estimate methods

2 資料處理

2.1 資料來源及可靠性分析

表1 用Z2與實測Z得到的轉(zhuǎn)換函數(shù)之比較Table 1 Comparison between transfer functions calculated by using Z2 and measured Z value

注： 表中某些周期值實際代表一定的周期范圍, 即34.4（32.0—36.57）， 25.9（23.27—28.44）， 19.2（17.07—21.33）， 14.1（12.19—16.0）， 10.1（8.53—11.64）， 括號內(nèi)為其范圍.

本文使用經(jīng)高斯濾波預處理的分鐘值， 計算轉(zhuǎn)換函數(shù)時沒有進行高通濾波． 由于一些小的干擾可能沒有被發(fā)覺， 或沒有辦法識別、 消除， 因此短周期的結(jié)果并不穩(wěn)定． 同時由于采取端點去傾， 數(shù)據(jù)沒有零均值化， 使得存在快速傅里葉變換的非零直流項（第一項）． 由于本文快速傅里葉變換的樣本長度不算很長， 所以還會有譜線擴散問題． 當3個分量的時間曲線大多集中在起、 終點連線一側(cè)時， 非零直流項有時可能會較大， 以致于影響長周期的結(jié)果． 因此本文在研究轉(zhuǎn)換函數(shù)的時間變化時， 只取周期為10—64min的值.

2.2 事件選取

用數(shù)字化資料計算復轉(zhuǎn)換函數(shù)， 需利用地球變化磁場的短周期成分． 鑒于要進行譜分析， 因此需選用有連續(xù)擾動的時段， 也可以包含急始、 灣擾和脈動等． 而量圖求a, b的方法僅選典型的地磁短周期變化如急始、 類急始、 灣擾和孤立的擾動． 根據(jù)作者的經(jīng)驗和數(shù)據(jù)的預處理方法以及當前儀器的狀態(tài)， 本文提出如下事件選取原則： ① 選取目標周期成分， 根據(jù)過去的工作總結(jié)出一些地震出現(xiàn)異常的周期（表2）， 可供參考． ② 避開3個分量的日變化， 盡量在夜間的時段選?。?③ 避開和識別包括儀器、 環(huán)境在內(nèi)的各種干擾, 有些干擾持續(xù)時間較短， 可以用線性內(nèi)插的辦法改正． ④ 一般不選特別激烈擾動， 因激烈的擾動可能不能滿足源場準均勻的假設． 由于源場不均勻產(chǎn)生的源場效應會造成估計值的漲落, 為此盡量選全球同時發(fā)生的事件， 不選地方性擾動如鉤擾． ⑤ 遇到缺數(shù)出現(xiàn)99999時， 如時間不長可以用線性內(nèi)插的辦法補數(shù)． ⑥ 起止點要選在比較平坦的位置而不是快速變化的途中， 要兼顧H和D的變化， 還要兩頭兼顧． 由于是端點線性去傾， 如果起止位置處在一個相近水平則最好, 如不能滿足上述條件， 則起止點也要選在拐點（停頓）或極值點的位置（圖3）． ⑦ 需要注意識別Z的較長周期成分是否由于H和D的變化所引起， 如果不是， 則會影響長周期的轉(zhuǎn)換函數(shù)結(jié)果， 這種情況我們曾遇到過.

表2 幾個震例中異常的最大周期Table 2 The maximum periods of anomalies for several examples of earthquakes

注：ΔZ/ΔH用的是前沿時間， 大體相當于半周期．

本文采取3個臺同時挑選事件的辦法． 為便于識別Z的變化， 作圖3時3個臺的Z都放大了1倍. 可以看出, 3個臺的H和D變化大體相同， 但Z變化不同； 還可看出重慶臺的Z受H的影響較明顯， 成都臺的Z受D的影響較明顯， 而西昌臺受H和D的共同影響較明顯． 這反映了地下電性結(jié)構(gòu)的差異． 根據(jù)龔紹京（1983）對地磁短周期事件時間序列的分析， 選m=12—16計算一組轉(zhuǎn)換函數(shù)， 多數(shù)都是13—15個事件． 由于是分鐘值， 所能挑選的事件數(shù)和所能求出的周期都會受到限制， 無法求出更短周期的轉(zhuǎn)換函數(shù)， 好在以往震例顯示10—48min周期的轉(zhuǎn)換函數(shù)也能很好地檢測到異常（表2）． 每組轉(zhuǎn)換函數(shù)所占的時間跨度視事件多少而變， 一般7—10天求一組， 也有5天可求一組的； 數(shù)據(jù)不好時半個月或更長時間才能求出一組． 以每個時間跨度的中點日期代表該組轉(zhuǎn)換函數(shù)的時間點， 誤差僅有半天.

圖3 事件選取示例（事件起止時間： 2009年7月20日14:12—18:28．

3 結(jié)果分析

3.1 成都臺兩種時間變化曲線的比較

本文首先研究了最初得到的數(shù)據(jù)， 繪制出Au， Av， Bu和Bv的時間變化曲線． 10—64min的周期范圍可以畫出8條時間變化曲線， 但圖形顯得“很擠”. 為此選擇了6條曲線， 其對應的周期從下至上分別是： 51.2， 42.7， 34.4， 29.5， 19.2和10.1min（圖4）． 除Au在地震前后有所變化外（圖4a）， 其它幾條曲線都未顯示出明顯變化， 且Bu和Bv的漲落較大， 限于篇幅僅列出Au圖． 圖4a中， 2007—2008年用M15預處理分鐘值， 2009年開始用GM4預處理分鐘值， 其中2008年1月26日至年底的Z數(shù)據(jù)因受潮不可靠． 由圖4a可看出， 42.7， 34.4和29.5min這3個周期成分在地震發(fā)生前后有明顯的“鼓包”. 只是由于儀器受潮， 不能確信上述“鼓包”就是異常， 需用Z2數(shù)據(jù)繪出圖4b以相互佐證.

圖4b中， 2007年用M15預處理數(shù)據(jù)， 2008年1月—5月19日用計算的Z2值， 5月20日09時07分開始用GM4預處理數(shù)據(jù)． 可以看出, 圖4a與圖4b在2008年的形態(tài)是不同的． 但仔細查看， 會發(fā)現(xiàn)在地震發(fā)生和緊隨其后不長的時間內(nèi)（即圖4a出現(xiàn)“鼓包”的時段內(nèi)）， Au曲線有可察覺的上升， 即Au絕對值減?。?而從圖4b還可看到震前有一不太深的“凹陷”（42.7， 34.4， 29.5和10.1min）， 即Au的絕對值變大． 有意思的是， 蘆山地震前似乎也有一點“凹陷”．

3.2 3個地磁臺時間變化曲線的比較與分析

圖5給出了西昌臺和重慶臺的Au時間曲線， 圖中數(shù)據(jù)均來自GM4． 可以看出， 2007—2009年上半年兩個臺站的Au值都很平穩(wěn)， 沒有異常． 西昌臺自2009年下半年以后Au值漲落較大， 說明有來自儀器和環(huán)境的干擾． 對于成都臺的結(jié)果, 作者認為異常應遵從成片原則： ① 要有多個周期同時出現(xiàn)異常； ② 要有連續(xù)多個點同時出現(xiàn)異常． 圖4a的異常是明顯的， 但地磁臺網(wǎng)中心認為那段資料不可靠， 盡管本文在選擇事件時已排除了可以察覺到的不正常情況， 但仍可能有人眼不能察覺的因素影響了數(shù)據(jù)． 圖4b的結(jié)果是對圖4a的佐證， 盡管圖形不太一樣， 但大部分周期的時間曲線均在汶川地震前有一個下降的趨勢， 并在地震時回升， 以10.1， 34.3和42.7min3個周期較明顯． 圖4中的曲線下降表明Au的絕對值增大， 說明地下電性的橫向不均勻性增加． 地震時Au絕對值變小， 曲線回升．

圖5 （a） 西昌臺Au時間曲線（Au為正值）; （b） 重慶臺Au時間曲線（Au為正值， 2011年11月以后為重慶附近的仙女山臺資料）

為了更好地考察成都臺轉(zhuǎn)換函數(shù)的變化， 本文用成都臺GM4資料的Au和Bu求出實帕金森矢量的長度Lr和方位角Fr， 并求出Au， Bu， Lr及Fr的均值（圖6）， 求均值所取的時間跨度參考圖4b． 圖6中， 在汶川和蘆山地震前成都臺均表現(xiàn)為Au下降（絕對值增大）， 地震時及震后（2008年5—6月和2013年4—7月）Au上升． 對于汶川地震， 25.9—64.0min周期的Au是在地震發(fā)生時上升至最大， 而10.1—19.2min周期的Au是在震后達到最大; Bu的大多數(shù)周期在地震時有所增大， 而震前相對較低; Fr表現(xiàn)為在兩個地震時或震后減小， 對汶川地震， 短周期成分的減小要滯后一點， 表明深層Fr的變小先于淺層． 帕金森矢量長度Lr規(guī)律性較差， 這與D的記錄值為“分”， 由“分”換算成“nT”的換算系數(shù)接近“10”， 因而Bu值誤差較大有關． 但還是能找到一些規(guī)律， 如Lr較大的周期是25.9， 34.4和42.7min， 以34.4min最大， 說明這些周期對應的深度電性結(jié)構(gòu)橫向不均勻性較大． 又如10.1， 14.1， 34.4和42.7min周期的Lr在汶川地震前已變大， 19.2， 25.9， 51.2和64.0min則在震時變大， 說明在地震前或地震時電性的橫向不均勻性有所增加． 帕金森矢量的方位角Fr變小， 表明帕金森矢量的指向往地震發(fā)生的方向移動． 也說明上述括號中兩個時段的Au絕對值變小是由于帕金森矢量的方向變化所致． 成都臺的Lr值變大和Fr變小是重要的異常標志． 然而， 這只是一種短期的前兆， 是否有更長期的異常變化還有待研究．

圖6 成都臺轉(zhuǎn)換函數(shù)A的實部Au, B的實部Bu以及實帕金森矢量的長度Lr和方位角Fr的均值變化

3.3 成都臺轉(zhuǎn)換函數(shù)和帕金森矢量的長期變化

成都、 西昌、 重慶臺2007年5月以前無磁通門磁力儀資料． 由于客觀條件的限制， 本文作者目前未能收集磁變儀資料并將其數(shù)字化， 以求出復轉(zhuǎn)換函數(shù)， 只好與以前在成都臺用量圖方法做的a， b對比（龔紹京等， 1989）． 該文量取地磁短周期變化事件（如急始、 類急始、 各類孤立的擾動等）的前沿時間Δt及相應的ΔZ，ΔH，ΔD， 其持續(xù)時間Δt只能與本文14.1和10.1min周期的結(jié)果（對應的周期范圍是8.5—16.0min）對比． 為考察是否存在長期變化， 表3列出了成都臺2009—2011年的均值.

表3 成都臺2009—2011年Au, Bu, Lr和Fr的均值及其誤差Table 3 The average values and corresponing errors of Au, Bu, Lr， Frat the Chengdu station in 2009—2011

注: 誤差dAu,dLr,dFr和dBu是2009—2011年3年的平均誤差， 而不是平均值的標準差．

圖7 成都、 西昌和重慶臺的帕金森矢量（Δt=3—10 min）

成都、 西昌、 重慶3個臺不同周期范圍的帕金森矢量示于圖7． 可以看出， 成都臺的帕金森矢量指向龍門山斷裂帶并有較長期的變化， 黑色箭頭是20世紀80年代的結(jié)果， 黃色箭頭為表3的結(jié)果， 兩箭頭相差約27°． 重慶臺的帕金森矢量比另兩個臺都長， 其指向背離四川盆地． 西昌臺的帕金森矢量指向西南， 不同周期帕金森矢量的指向有差別．

帕金森矢量是一個用地磁方法反映地下電性結(jié)構(gòu)的參量， 它的長度越大， 表示地下電性橫向不均勻的程度越高． 其方向指向?qū)щ娐矢叩囊环剑?例如它會垂直于海岸線總的走向或指向附近的深海． 在北半球， 它的畫法是由磁南右旋Fr角． 在地理坐標中需考慮偏角的年均絕對值Ds， 例如Ds為-1.7°、 Fr=108.3°表示與磁北的夾角是-71.7°， 則帕金森矢量與地理北的夾角是-1.7°-71.7°=-73.4°． 帕金森矢量的物理意義和求法見Parkinson（1959， 1962）與龔紹京（1987）等文章． 用1979—1987年成都臺的a與b均值求出帕金森矢量， 其長度L=0.13， 方位角F=171.5°， 與地理北的夾角為-10.1°． 而2009—2011年成都臺周期為8.5—16.0min的實帕金森矢量長度Lr=0.182±0.021， 方位角Fr=144.8°±8.3°， 與地理北的夾角為-36.9°， 變化超出了2倍誤差． 該結(jié)果與龔紹京和劉雙慶（2012）的初步結(jié)果吻合， 只是該文用的是未經(jīng)預處理的分鐘值， 本文用的是經(jīng)過預處理的分鐘值資料．

量圖方法得出的a與b沒有嚴格的周期概念， 因此這種對比是很粗糙的． 為了進一步驗證轉(zhuǎn)換函數(shù)和帕金森矢量的長期變化， 作者在計算機上處理了2011年4—6月的資料， 量取了Δt，ΔD，ΔH，ΔZ值， 67組數(shù)求出的結(jié)果為 a=-0.1123， b=0.1412， L=0.1776， F=128.5°． 由a和b求得的帕金森矢量與地理北的夾角為-53.2°， 與1979—1987年在磁照圖上的量圖結(jié)果相比， b， L和F都有很大變化， 甚至比Fr的變化還大（表4）．

表4 成都臺地磁短周期變化參量的長期變化Table 4 The long-term changes in geomagnetic short-period variation parameters at the Chengdu station

4 討論與結(jié)論

綜合上述分析， 本文得出以下幾點看法：

1） 無論是求a與b還是求A與B， 都不應該進行變換， 即不能在等式兩邊除同一變量然后作回歸分析． 對求a和b要進行二元回歸分析． 對復轉(zhuǎn)換函數(shù)A和B， 最早和大量使用的都是式（4）． 如用多元回歸分析， 所求系數(shù)和方程數(shù)都將增加一倍．

2） 無論是用量圖方法得到的a， b， L， F， 還是用譜分析方法得到的Au， Bu， Lr， Fr， 與20世紀80年代成都臺的結(jié)果相比較， 都顯示轉(zhuǎn)換函數(shù)和帕金森矢量有明顯的長期變化． b， L和F的變化比較大， 而a的長期變化不大． L增大說明這些年該地區(qū)地下電性橫向不均勻程度增加， F變小說明帕金森矢量的方向朝地震發(fā)生的區(qū)域移動．

轉(zhuǎn)換函數(shù)的長期變化已有柿崗臺（Yanagihara， 1972）和廣州臺的例子（林美， 龔紹京， 1991）． 前者與1923年關東大地震有關， 后者處理了1960—1987年共28年資料， 其長期變化也許與菲律賓板塊和歐亞板塊的相對運動或是大陸架隆起因而海水變淺有關.

為了準確而嚴格地考察成都臺各種周期的長期變化及變化過程， 需將1975—2007年的磁變儀資料數(shù)字化， 計算復轉(zhuǎn)換函數(shù)并求出實、 虛帕金森矢量．

3） 汶川和蘆山地震前后Au有一點略可察覺的短期變化. 從圖6可看出, Lr, Bu和Fr也在同期有所變化， 但這些變化都不顯著． 作者過去的研究表明， 垂直場轉(zhuǎn)換函數(shù)在唐山地震前后有5—7年的中長期變化（龔紹京等， 1985， 1997）， 但并未發(fā)現(xiàn)明顯的短期前兆． 只是水平場轉(zhuǎn)換函數(shù)有十分明顯的短期異常（龔紹京等， 2001）． 因此汶川地震前后Au等的短期變化量不大是可以理解的．

4） 西昌和重慶臺在這兩次大地震前后沒有可察覺的異常變化， 這從一個側(cè)面說明成都臺在地震前后的可察覺變化也許與地震有關．同時， 唐山地震和本文結(jié)果均說明地磁轉(zhuǎn)換函數(shù)的異常范圍是不大的． 其異常的范圍還與構(gòu)造帶的產(chǎn)狀或余震區(qū)的形狀有關． 例如離唐山地震震中180km的白家疃臺無異常， 但震中距110km的天津青光臺卻有異常． 這是因為唐山地震的斷裂帶和余震區(qū)呈北東—西南長條分布． 成都臺對相距180km的松潘地震有反應， 也是因為該地震離龍門山斷裂帶較近． 但西昌和重慶臺則離龍門山斷裂帶較遠．

5） 本文成都、 西昌、 重慶3個地磁臺的Au值及以前的一些研究結(jié)果均表明， 轉(zhuǎn)換函數(shù)值在正常情況下都比較穩(wěn)定， 即使出現(xiàn)‘突跳’等情況， 其變化幅度也不會超出一個量級． 這與袁寶珠等（2009）的結(jié)果形成鮮明對比． 對其所引用的方法和計算公式作者已經(jīng)作過分析（龔紹京等， 2012）． 由于轉(zhuǎn)換函數(shù)反映的是地下的電性結(jié)構(gòu)， 在無異常、 無干擾、 無儀器故障的情況下給出的轉(zhuǎn)換函數(shù)值表現(xiàn)出‘穩(wěn)定’才是正常的．

感謝中國地震局地球物理研究所地磁臺網(wǎng)中心主任楊冬梅研究員對我們工作的大力支持， 特別感謝地磁臺網(wǎng)中心張素琴副研究員不厭其煩地多次為我們拷貝數(shù)據(jù)．

陳伯舫. 1998. 關島8.1級大地震和地磁轉(zhuǎn)換函數(shù)時間變化的關系[J]. 地震學報， 20（2）： 217--219.

ChenPF. 1998.TheGuangreatearthquake（Msz=8.1）andtimechangesingeomagnetictransferfunctions[J]. Acta Seismologica Sinica, 20（2）： 217--219 （inChinese）.

龔紹京. 1983. 青光臺地磁短周期事件的時間序列分析[J]. 地震， （1）: 6--10.

GongSJ. 1983.TheanalysisonthetimeseriesofshortperiodgeomagneticeventatQingguangseismographicstation[J]. Earthquake, （1）: 6--10 （inChinese）.

龔紹京, 吳占峰, 蔣邦本. 1984. 地磁場瞬時擾動ΔZ/ΔH的異常變化[J]. 地震科學研究， （5）: 48--51.

GongSJ,WuZF,JiangBB. 1984.TheanomalouschangesofΔZ/ΔHoninstantaneousdisturbanceofgeomagneticfield[J]. Research of Seismological Science, （5）: 48--51 （inChinese）.

龔紹京. 1987. 廣東省地磁臺的帕金森矢量及廣州臺系數(shù)在河源地震前后的時間變化[J]. 地震研究， 10（5）: 575--582.

GongSJ． 1987.ParkinsonvectorsatthegeomagneticstationsinGuangdongProvinceandthetime-dependentchangesoftheirratioatGuangzhoustationbothbeforeandaftertheHeyuanearthquake[J]. Journal of Seismological Research, 10（5）: 575--582 （inChinese）.

龔紹京， 王伯維， 于彬. 1989. 松潘地震前后地磁轉(zhuǎn)換函數(shù)的變化[C]∥松潘地震預報學術討論會文集. 北京： 地震出版社: 95--98.

GongSJ，WangBW，YuB. 1989.ChangesofgeomagnetictransferfunctionsbeforeandafterSongpanearthquake[C]∥The Collected Works for Symposium About Forecasting Songpan Earthquake.Beijing:SeismologicalPress: 95--98 （inChinese）.

龔紹京， 楊桂君， 田山， 于彬. 1991. 菏澤5.9級地震前后菏澤臺轉(zhuǎn)換函數(shù)隨時間變化的研究: 兼與王锜同志商榷[J]. 地震學報， 13（1）: 113--120.

GongSJ,YangGJ,TianS,YuB. 1992.ResearchonthetimechangesoftransferfunctionsatHezeobservatorybeforeandafterHezeearthquakeofmagnitude5.9:AdiscussionwithWangQi[J]. Acta Seismologica Sinica, 5（1）: 187--195.

龔紹京， 陳化然， 張翠芬， 馬淑芹， 楊桂君. 1997. 地磁水平場轉(zhuǎn)換函數(shù)在唐山地震前的異常反應[J]. 地震學報， 19（1）: 51--58.

GongSJ,ChenHR,ZhangCF，YangGJ，MaSQ. 1997.TheanomalousreactionsofthegeomagnetichorizontalfieldtransferfunctionsbeforeTangshanearthquake[J]. Acta Seismologica Sinica, 10（1）: 61--70.

龔紹京， 田真麗， 戚成柱， 何淑敏， 閻嘵梅， 陳化然， 栗連弟. 2001. 地磁水平場轉(zhuǎn)換函數(shù)的短期前兆[J]. 地震學報， 23（3）: 280--288.

GongSJ,TianZL,QiCZ,HeSM,YanXM,ChenHR,LiLD. 2001.Shorttermprecursorofthegeomagnetichorizontalfieldtransferfunctions[J]. Acta Seismologica Sinica, 14（3）： 293--302.

龔紹京， 劉雙慶. 2012. 汶川M8.0地震前帕金森矢量的變化[G]∥王子昌先生誕辰百年紀念文集. 北京： 北京大學出版社： 45--51.

GongSJ,LiuSQ. 2012.ThechangeofParkinsonvectorbeforetheWenchuanM8.0earthquake[G]∥The Collected Works for Centennial of Wang Zichang Professor.Beijing：PekingUniversityPress： 45--51 （inChinese）.

龔紹京， 馬驥， 劉雙慶， 栗連弟. 2012. 對《轉(zhuǎn)換函數(shù)與汶川大地震關系的初步研究》一文的分析[J]. 國際地震動態(tài)， （8）: 20--26.

GongSJ,MaJ,LiuSQ,LiLD. 2012.Commenton“StudyontherelationshipbetweentheWenchuanstrongearthquakeandthegeomagnetictransferfunction”byBaozhuYuanetal[J]. Recent Developments in World Seismology, （8）: 20--26 （inChinese）.

林美， 龔紹京. 1991. 廣州臺轉(zhuǎn)換函數(shù)的長期變化和季節(jié)變化[J]. 地震學報， 13（4）: 480--488.

LinM,GongSJ. 1992.SecularandseasonalchangesintransferfunctionatGuangzhougeomagneticobservatory[J]. Acta Seismologica Sinica, 5（3）: 587--595.

袁寶珠， 陳化然， 張素琴， 李琪， 楊冬梅， 朱榮， 劉曉燦. 2009. 地磁轉(zhuǎn)換函數(shù)與汶川大地震關系的初步研究[J]． 國際地震動態(tài), （7）: 69--75.

YuanBZ,ChenHR，ZhangSQ，LiQ，YangDM,ZhuR,LiuXC. 2009.StudyontherelationshipbetweentheWenchuanstrongearthquakeandthegeomagnetictransferfunction[J]. Recent Developments in World Seismology, （7）: 69--75 （inChinese）.

BeamishD. 1982.Ageomagneticprecursortothe1979Carlisleearthquake[J]. Geophys J R astr Soc, 68（2）: 531--543.

ChapmanS，BartelsJ. 1940. Geomagnetism[M].Oxford:ClarendonPress: 1049.

EgbertGD,BookerJR． 1986．Robustestimationofgeomagnetictransferfunctions[J]. Geophys J R astr Soc, 87（1）: 173--194.

EverettJE,HyndmanRD. 1967．Geomagneticvariationsandelectricalconductivitystructureinsouth-westernAustralia[J]． Phys Earth Planet Inter, 1（1）: 24--34.

GongSJ． 1985．Anomalouschangesintransferfunctionsandthe1976Tangshanearthquake（MS=7.8）[J]. J Geomag Geoelectr, 37（4）： 503--508.

GongSJ,ChenPF,YangGJ. 1991.Researchonthetimechangesofinter-stationtransferfunctionsforthehorizontalgeomagneticfieldandtheirrelationshipwiththeHualian7.6earthquakeinTaiwanregion[C]∥International Conference of Seismicity in Eastern Asia.HongKong,October1991.

MiyakoshiJ. 1975.SecularvariationofParkinsonvectorsinaseismicallyactiveregionofMiddleAsia[J]. J Fac General Education, Tottori Univ, 8: 209--218.

ParkinsonWD． 1959．Directionsofrapidgeomagneticfluctuations[J]． Geophys J R astr Soc， 2（1）： 1--14.

ParkinsonWD． 1962．Theinfluenceofcontinentsandoceansongeomagneticvariations[J]. Geophys J R astr Soc, 6（4）: 441--449.

RikitakeT. 1979.Changesinthedirectionofmagneticvectorofshort-periodgeomagneticbeforethe1972Sitka,Alaska,earthquake[J]. J Geomg Geoelectr, 31（4）: 441--448.

SchmuckerU. 1970.AnomaliesofgeomagneticvariationsinthesouthwesternUnitedStates[J]． Bull Scripps Inst Oceanogr, 13： 1--165.

UntiedtJ. 1970.ConductivityanomaliesincentralandsouthernEurope（Appendix）[J]. J Geomag Geoelectr, 22（1/2）: 131--149.

YanagiharaK. 1972.SecularvariationoftheelectricalconductivityanomalyinthecentralpartofJapan[J]. Memo Kakioka Mag Obs, 15: 1--11.

Time changes in the geomagnetic transfer functions at Chengdu, Xichang and Chongqing stations from May 2007 to December 2013

（EarthquakeAdministrationofTianjinMunicipality,Tianjin300201,China）

This paper discusses the method for computing geomagnetic transfer functions in detail and introduces the error formulae of the real and imaginary parts of complex transfer functions. And then the transfer functionsAandBas well as Parkinson vectors are calculated by using the digital geomagnetic data recorded at the stations Chengdu, Xichang and Chongqing from May 2007 to December 2013. The results show that there are perceptible small bulges at Chengdu station before and after the 2008 WenchuanMS8.0 earthquake, but the other two stations do not exhibit perceptible anomalies. Comparison with the results in the 1980s indicates that obvious long-term variations of thebvalue as well as the length and direction of Parkinson vectors are observed at the station Chengdu.

geomagnetic transfer function; complex least squares method; multiple regression analysis; robust estimation; Parkinson vector; earthquake anomaly

10.11939/jass.2015.01.013.

天津市地震局局長科研基金資助.

2014-03-04收到初稿， 2014-06-03決定采用修改稿.

e-mail: caogong2003@hotmail.com

10.11939/jass.2015.01.013

P315.72+1

A

龔紹京, 劉雙慶, 張明東. 2015. 2007年5月—2013年12月成都、 西昌和重慶臺地磁轉(zhuǎn)換函數(shù)的時間變化. 地震學報, 37（1）: 144--159.

Gong S J, Liu S Q, Zhang M D. 2015. Time changes in the geomagnetic transfer functions at Chengdu, Xichang and Chongqing stations from May 2007 to December 2013.ActaSeismologicaSinica, 37（1）: 144--159. doi:10.11939/jass.2015.01.013.