龔紹京 劉雙慶 張明東

（中國(guó)天津300201天津市地震局）

?

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

（中國(guó)天津300201天津市地震局）

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

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

引言

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

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

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

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

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

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

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

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

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

（1）

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

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

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

（2）

（3）

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

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

（4）

其中

（5）

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

（6）

同理也可推導(dǎo)出B實(shí)部與虛部誤差公式為

（7）

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

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

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

（8）

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

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

圖1 成都臺(tái)2007年9月22—28日（a）與2009年7月11—21日（b）的周期響應(yīng)曲線

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

1.4 穩(wěn)健估計(jì)方法效果評(píng)價(jià)

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

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

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

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 資料來(lái)源及可靠性分析

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

注： 表中某些周期值實(shí)際代表一定的周期范圍, 即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）， 括號(hào)內(nèi)為其范圍.

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

2.2 事件選取

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

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

注：ΔZ/ΔH用的是前沿時(shí)間， 大體相當(dāng)于半周期．

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

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

3 結(jié)果分析

3.1 成都臺(tái)兩種時(shí)間變化曲線的比較

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

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

3.2 3個(gè)地磁臺(tái)時(shí)間變化曲線的比較與分析

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

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

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

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

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

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

表3 成都臺(tái)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年的平均誤差， 而不是平均值的標(biāo)準(zhǔn)差．

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

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

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

量圖方法得出的a與b沒(méi)有嚴(yán)格的周期概念， 因此這種對(duì)比是很粗糙的． 為了進(jìn)一步驗(yàn)證轉(zhuǎn)換函數(shù)和帕金森矢量的長(zhǎng)期變化， 作者在計(jì)算機(jī)上處理了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年在磁照?qǐng)D上的量圖結(jié)果相比， b， L和F都有很大變化， 甚至比Fr的變化還大（表4）．

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

4 討論與結(jié)論

綜合上述分析， 本文得出以下幾點(diǎn)看法：

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

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

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

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

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

4） 西昌和重慶臺(tái)在這兩次大地震前后沒(méi)有可察覺(jué)的異常變化， 這從一個(gè)側(cè)面說(shuō)明成都臺(tái)在地震前后的可察覺(jué)變化也許與地震有關(guān)．同時(shí)， 唐山地震和本文結(jié)果均說(shuō)明地磁轉(zhuǎn)換函數(shù)的異常范圍是不大的． 其異常的范圍還與構(gòu)造帶的產(chǎn)狀或余震區(qū)的形狀有關(guān)． 例如離唐山地震震中180km的白家疃臺(tái)無(wú)異常， 但震中距110km的天津青光臺(tái)卻有異常． 這是因?yàn)樘粕降卣鸬臄嗔褞Ш陀嗾饏^(qū)呈北東—西南長(zhǎng)條分布． 成都臺(tái)對(duì)相距180km的松潘地震有反應(yīng)， 也是因?yàn)樵摰卣痣x龍門(mén)山斷裂帶較近． 但西昌和重慶臺(tái)則離龍門(mén)山斷裂帶較遠(yuǎn)．

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

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

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

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

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

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

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

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

龔紹京. 1987. 廣東省地磁臺(tái)的帕金森矢量及廣州臺(tái)系數(shù)在河源地震前后的時(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]∥松潘地震預(yù)報(bào)學(xué)術(shù)討論會(huì)文集. 北京： 地震出版社: 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級(jí)地震前后菏澤臺(tái)轉(zhuǎn)換函數(shù)隨時(shí)間變化的研究: 兼與王锜同志商榷[J]. 地震學(xué)報(bào)， 13（1）: 113--120.

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

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

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

龔紹京， 田真麗， 戚成柱， 何淑敏， 閻嘵梅， 陳化然， 栗連弟. 2001. 地磁水平場(chǎng)轉(zhuǎn)換函數(shù)的短期前兆[J]. 地震學(xué)報(bào)， 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]∥王子昌先生誕辰百年紀(jì)念文集. 北京： 北京大學(xué)出版社： 45--51.

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

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

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

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

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

袁寶珠， 陳化然， 張素琴， 李琪， 楊冬梅， 朱榮， 劉曉燦. 2009. 地磁轉(zhuǎn)換函數(shù)與汶川大地震關(guān)系的初步研究[J]． 國(guó)際地震動(dòng)態(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.

天津市地震局局長(zhǎng)科研基金資助.

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月成都、 西昌和重慶臺(tái)地磁轉(zhuǎn)換函數(shù)的時(shí)間變化. 地震學(xué)報(bào), 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.