管貽亮 張 玲 李希亮 董曉娜

1 山東省地震局，濟(jì)南市港西路2066號，250100

地磁轉(zhuǎn)換函數(shù)可以表達(dá)地磁水平分量和垂直分量之間的線性轉(zhuǎn)換關(guān)系，能夠反映地下電性結(jié)構(gòu)的橫向不均勻性。地震地磁研究中轉(zhuǎn)換函數(shù)的時序變化可以反映孕震區(qū)的電性結(jié)構(gòu)變化[1]，用于震前微觀異常信息的提取。諸多學(xué)者對轉(zhuǎn)換函數(shù)進(jìn)行深入研究[2-5]，獲取了豐富的震例經(jīng)驗。感應(yīng)矢量最早用于分析海岸效應(yīng)和低阻異常體的影響[6]，Schmucker[7]最先確立感應(yīng)矢量的計算公式，明確地磁水平分量和垂直分量之間的關(guān)系。由于不受局部異常體的影響，感應(yīng)矢量可更加直觀地展現(xiàn)地下結(jié)構(gòu)的電性差異，被廣泛應(yīng)用于在電磁測深中對宏觀電性結(jié)構(gòu)的研究[8-9]。但在地震地磁觀測中，感應(yīng)矢量和區(qū)域電性結(jié)構(gòu)的研究較少，王橋等[10]和龔紹京等[11]對部分地磁臺站資料進(jìn)行感應(yīng)矢量計算，并對其時序特征和構(gòu)造指示意義進(jìn)行研究。

本文計算山東省地磁臺陣15個臺站的轉(zhuǎn)換函數(shù)和感應(yīng)矢量，并對其時序變化特征進(jìn)行多臺站對比分析，對區(qū)域電性結(jié)構(gòu)和構(gòu)造指示關(guān)系等進(jìn)行研究。所選臺站均為固定連續(xù)觀測且數(shù)量較多，其數(shù)據(jù)質(zhì)量和穩(wěn)定性明顯優(yōu)于傳統(tǒng)的大地電磁測深(MT)數(shù)據(jù)，適合開展感應(yīng)矢量的背景特征分析。本文研究可以為區(qū)域磁場變化、構(gòu)造活動性和地震危險性研究提供背景數(shù)據(jù)，后續(xù)可在該基礎(chǔ)上開展時移動態(tài)監(jiān)測，為地震活動趨勢判定提供地球物理依據(jù)。

1 數(shù)據(jù)來源與計算方法

1.1 數(shù)據(jù)來源

2016年山東省建成地磁臺陣，采用GM4磁通門磁力儀進(jìn)行秒值觀測，臺站分布較為均勻，2017年開始多數(shù)臺站數(shù)據(jù)趨于穩(wěn)定。本文基于數(shù)據(jù)有效率、背景噪聲等評價指標(biāo)選取數(shù)據(jù)質(zhì)量較好的15個臺站，臺站分布如圖1所示。其中沂沭斷裂帶(F2)附近分布有萊州、安丘、留山、玉皇山、陵陽、郯城6個臺站，濟(jì)南、長清、泰安、嘉祥、鄒城5個臺站分布于魯西隆起西側(cè)，萊州、北長山、文登、平度4個臺站分布于膠東隆起附近且距離海域較近。為保證結(jié)果的穩(wěn)定性，選取分鐘值預(yù)處理數(shù)據(jù)進(jìn)行轉(zhuǎn)換函數(shù)和感應(yīng)矢量的計算。

F1:五蓮-即墨-牟平斷裂帶；F2:沂沭斷裂帶(F2-1:昌邑-大店斷裂；F2-2:安丘-莒縣斷裂；F2-3:沂水-湯頭斷裂；F2-4:鄌郚-葛溝斷裂)；F3:聊城-蘭考斷裂；F4:齊河-廣饒斷裂；Ⅰ1:蘇魯造山帶；Ⅱ1:膠遼陸塊；Ⅱ2:魯西隆起；Ⅱ3:華北坳陷圖1 地磁觀測臺站分布

1.2 計算方法

地磁場垂直分量與2個水平分量之間的關(guān)系可表示為[7]：

Z(ω)=A(ω)X(ω)+B(ω)Y(ω)

(1)

式中，所有量均為復(fù)數(shù)；Z(ω)、X(ω)、Y(ω)分別為地磁垂直分量、南北分量和東西分量在某頻率ω上的頻譜；A、B為地磁轉(zhuǎn)換函數(shù)，其表現(xiàn)形式與MT中阻抗張量的估算方式一致?？墒褂米钚《朔ㄊ箽埐钰呌谧钚≈担錃埐瞀趴杀硎緸椋?/p>

(2)

式中，k為事件序號；N為事件總數(shù)；Z、X、Y為各分量的復(fù)數(shù)頻譜；A、B為復(fù)數(shù)轉(zhuǎn)換函數(shù)，其中A=Ar+iAi,B=Br+iBi。由于地磁連續(xù)觀測可能會受到較多干擾，導(dǎo)致計算結(jié)果出現(xiàn)較大偏差，因此可采用加權(quán)最小二乘法對強干擾數(shù)據(jù)進(jìn)行加權(quán)，降低其參與計算的比重，減小對計算結(jié)果的影響。

基于轉(zhuǎn)換函數(shù)可計算得到實感應(yīng)矢量Gr和虛感應(yīng)矢量Gi，具體表示為：

Gr=Ari+Brj,Gi=Aii+Bij

(3)

式中，i和j分別為指向正北和正東方向的單位矢量。在實際應(yīng)用時，為使實感應(yīng)矢量指向電流聚集方向，應(yīng)對實感應(yīng)矢量取反向。文中感應(yīng)矢量均指帕金森矢量，其大小和方向為：

α=arctan(Br/Ar)，

(4)

式中α、L分別為帕金森矢量的方位角和大小，I為地磁變化矢量優(yōu)勢平面傾角。

首先對數(shù)據(jù)進(jìn)行預(yù)處理，剔除有效率低于20%的數(shù)據(jù)，對超過5倍均方差的數(shù)據(jù)作刪除處理；然后將每天的預(yù)處理分鐘數(shù)據(jù)分4段進(jìn)行離散傅里葉變換；最后對每個頻點數(shù)據(jù)進(jìn)行整合，單獨計算各頻點的轉(zhuǎn)換函數(shù)。其中，頻譜基于地磁分析預(yù)報軟件進(jìn)行計算，其他操作均用Python語言編程實現(xiàn)。

2 結(jié)果分析

2.1 轉(zhuǎn)換函數(shù)時頻分析

基于加權(quán)最小二乘法對山東省地磁臺陣數(shù)據(jù)進(jìn)行轉(zhuǎn)換函數(shù)計算。為保證結(jié)果的穩(wěn)定性、減小計算誤差、便于開展背景特征研究，選取2019-05～2020-05共1 a的數(shù)據(jù)參與計算。分別以Ar和Br表示轉(zhuǎn)換函數(shù)A和B的實部，15個臺站的計算結(jié)果如圖2所示，圖中大部分臺站的轉(zhuǎn)換函數(shù)曲線比較光滑、計算誤差較小，各頻點結(jié)果的連續(xù)性較好，說明計算結(jié)果具有穩(wěn)定性。

從轉(zhuǎn)換函數(shù)曲線可以看出，20 min周期內(nèi)轉(zhuǎn)換函數(shù)的數(shù)值波動較大，誤差也相對較大，Ar和Br之間沒有明顯的趨勢性變化，說明轉(zhuǎn)換函數(shù)對于淺層電性結(jié)構(gòu)的表現(xiàn)力較弱。原因主要有2點：一是地磁觀測在淺層受到的外部干擾較多，導(dǎo)致高頻數(shù)據(jù)計算誤差較大；二是由于山東省區(qū)域地表多為第四系覆蓋，沒有明確的構(gòu)造特征，轉(zhuǎn)換函數(shù)特征不明顯。沂沭斷裂帶附近的安丘、留山、玉皇山、陵陽4個臺站的轉(zhuǎn)換函數(shù)隨周期呈波動性變化，沒有趨勢性特征。相比于魯西地區(qū)，上述4個臺站的變化更加復(fù)雜，其中玉皇山臺和陵陽臺相距19 km，安丘臺和留山臺相距23 km，2組臺站轉(zhuǎn)換函數(shù)的相關(guān)系數(shù)分別為-0.25和0.45，相關(guān)性較低，說明2組臺站的地下電性結(jié)構(gòu)和構(gòu)造特征存在較大差異。沂沭斷裂帶是由4條分支斷裂組成的深大斷裂，其內(nèi)部可能為深度破碎帶，結(jié)構(gòu)復(fù)雜。而魯西地區(qū)的泰安、長清、濟(jì)南等臺站的轉(zhuǎn)換函數(shù)隨周期變化較為平穩(wěn)，說明該地區(qū)地下電性結(jié)構(gòu)隨深度變化較小。

圖2 轉(zhuǎn)換函數(shù)Ar和Br周期變化曲線

圖3 泰安臺轉(zhuǎn)換函數(shù)時間變化曲線及頻譜分析

基于轉(zhuǎn)換函數(shù)研究地下電性結(jié)構(gòu)的長期緩慢變化，從而提取與孕震有關(guān)的磁異常信息，是地震地磁轉(zhuǎn)換函數(shù)研究的主要目的。為對長期變化趨勢和背景特征進(jìn)行分析，本文基于逐日滑動算法對泰安臺和郯城臺2011～2020年共10 a的數(shù)據(jù)進(jìn)行逐日滑動計算。以N為窗長(windows length, WL)計算得到第1 d的值，其后窗長不變，逐日向后滑動計算每天的轉(zhuǎn)換函數(shù)。為得到最優(yōu)的計算參數(shù)，分別以10 d、30 d、60 d、180 d為窗長對40 min周期的轉(zhuǎn)換函數(shù)及頻譜進(jìn)行對比分析，計算結(jié)果如圖3、4所示?？梢钥闯觯?0 d窗長的結(jié)果波動較大，難以進(jìn)行長趨勢分析；30 d窗長的轉(zhuǎn)換函數(shù)具有一定的年變趨勢，數(shù)據(jù)波動相對較小。窗長為10 d和30 d的頻譜曲線中含有較多的短周期成分，說明滑動計算并沒有將短周期信號濾掉。窗長為60 d和180 d的頻譜曲線以年周期變化為主，短周期成分較少，60 d窗長顯示出0.5 a周期變化。由于長期背景成分主要來自空間磁場變化，與地下電性結(jié)構(gòu)關(guān)系不大，因此本文分別計算泰安臺和郯城臺Ar、Br的相關(guān)系數(shù)(表1)。由表可知，窗長為30 d和60 d時2個臺站的相關(guān)系數(shù)最高，說明此時計算結(jié)果中長周期背景最突出。王橋等[10]選取10 d窗長計算得到泰安臺和郯城臺長期變化趨勢具有較大差異，與本文結(jié)論不一致，可能與窗長的選取有關(guān)。綜合時頻曲線變化特征和相關(guān)性分析結(jié)果可知，選取60 d窗長計算山東地區(qū)逐日滑動轉(zhuǎn)換函數(shù)更加合適。不同窗長的時頻分析均證實了轉(zhuǎn)換函數(shù)的年變特征，因此如果要從轉(zhuǎn)換函數(shù)時間序列中提取震磁異常信息或構(gòu)造活動信息，需對長周期背景成分進(jìn)行分析、提取和濾波，才能從強大的背景場中提取出微弱、可靠的異常信息。

分析60 d窗長的時序曲線發(fā)現(xiàn)，轉(zhuǎn)換函數(shù)的峰谷值出現(xiàn)的時間大約為1月和7月，呈明顯的季節(jié)性變化特征，Ar和Br呈反向變化，即Ar的峰值對應(yīng)Br的谷值。對比Ar和Br的曲線可以看出，Ar的年變特征較Br更加明顯，數(shù)據(jù)穩(wěn)定性更好；Br隨時間變化波動較大，年變特征被壓制，數(shù)據(jù)毛刺較多。由此可知，選取Ar結(jié)果參與異常分析更加可靠。該結(jié)論對后續(xù)時序分析和震前轉(zhuǎn)換函數(shù)異常信息的提取具有一定的指導(dǎo)意義。

圖4 郯城臺轉(zhuǎn)換函數(shù)時間變化曲線及頻譜分析

表1 不同窗長計算的泰安臺和郯城臺Ar、Br的相關(guān)系數(shù)

2.2 感應(yīng)矢量特征分析

基于轉(zhuǎn)換函數(shù)計算1 a尺度下(2019-05～2020-05)的感應(yīng)矢量。為保證曲線的連續(xù)性，矢量方向角度范圍設(shè)置為-360°～360°，15個臺站的計算結(jié)果如圖5所示。各頻點的感應(yīng)矢量指向電流聚集方向(高導(dǎo)體)，其大小可反映該深度下的橫向不均勻程度，矢量越大說明其電性結(jié)構(gòu)差異越大。不同周期的感應(yīng)矢量可反映不同深度下的電性結(jié)構(gòu)差異，感應(yīng)矢量隨周期的變化能夠反映出地下電性結(jié)構(gòu)的垂向變化。從各臺站感應(yīng)矢量曲線可以看出，矢量大小隨周期變化穩(wěn)定，矢量方向隨周期變化波動。大部分臺站矢量方位角在高頻(80 min以下)部分存在波動變化，在低頻部分主要呈趨勢性變化，整個頻段的方位角變化小于30°，說明山東省區(qū)域地下電性結(jié)構(gòu)的橫向相對差異隨深度變化不明顯，沒有較大的電性反轉(zhuǎn)或構(gòu)造轉(zhuǎn)折。各臺站大部分頻點矢量絕對值小于0.2，且隨頻率降低呈趨勢性減小，說明整體電性橫向不均勻性較小，橫向差異隨深度增加而減小，而淺層電性差異相對更加明顯。

圖5 帕金森矢量周期變化曲線

圖6 帕金森矢量方向統(tǒng)計

為了更加直觀地體現(xiàn)地下電性結(jié)構(gòu)的變化情況，以15°為區(qū)間對矢量方向進(jìn)行統(tǒng)計并繪制玫瑰圖(圖6)。圖中玫瑰花瓣的長度代表該方向區(qū)間的頻點在所有頻點中的占比，花瓣數(shù)則表示方向區(qū)間數(shù)量，花瓣越多說明該觀測點的電性方向越分散，主花瓣越長且其他花瓣越少則說明該點的垂向電性越均勻。結(jié)合感應(yīng)矢量曲線和玫瑰圖綜合分析可知，文登、長清、郯城、萊州、北長山、陵陽6個臺站的矢量方位角最集中，其主要方向占比均超過50%，說明其地下電性結(jié)構(gòu)隨深度變化較小，構(gòu)造指示單一且明確。其中北長山、萊州、文登3個臺站沿海，其感應(yīng)矢量均指向海洋，呈現(xiàn)出明顯的海岸效應(yīng)。位于沂沭斷裂帶東側(cè)的郯城、安丘、平度3個臺站的矢量主要指向NW向，位于斷層內(nèi)部的留山臺矢量主要指向NEE向。上述4個臺站矢量均近似垂直指向斷裂帶方向，說明沂沭斷裂帶是一個顯著的電性分界面，其兩側(cè)的膠東隆起和魯西隆起電性結(jié)構(gòu)差異較大。結(jié)合矢量大小可知，淺層矢量越大，說明斷層兩側(cè)淺層的電性結(jié)構(gòu)差異越明顯，距離深部越近電性結(jié)構(gòu)差異越小。泰安臺感應(yīng)矢量主要指向SEE向，大小在0.05以下且變化較小，該結(jié)果與龔紹京等[11]的結(jié)論基本一致，說明該區(qū)域電性結(jié)構(gòu)相對比較簡單，電性橫向差異較小。長清臺和濟(jì)南臺相距38 km，二者矢量指向一致且各頻段明確指向NNW向的齊廣斷裂，周期曲線高度相似，相關(guān)系數(shù)為0.9，說明其地下電性結(jié)構(gòu)基本一致，均處于同一構(gòu)造活動區(qū)。嘉祥臺和鄒城臺相距60 km，二者相關(guān)系數(shù)為0.68，2個臺站感應(yīng)矢量在40～80 min周期頻段內(nèi)具有明顯的方向轉(zhuǎn)折，尤其是嘉祥臺的矢量方向變化較大，高導(dǎo)體由SW向往NE向轉(zhuǎn)換，60 min周期后又向W轉(zhuǎn)換，說明該區(qū)域在整個研究頻率所對應(yīng)深度下的SW向為高導(dǎo)區(qū)。40～80 min周期對應(yīng)深度下的NE向存在較大范圍的高導(dǎo)異常體，使電流聚集方向發(fā)生反轉(zhuǎn)，可能對應(yīng)較大的高導(dǎo)構(gòu)造區(qū)，但由于感應(yīng)矢量較小，該高導(dǎo)區(qū)與鄰區(qū)的電性差異并不明顯。

3 結(jié) 語

1)山東地區(qū)的地磁轉(zhuǎn)換函數(shù)對淺層的表現(xiàn)特征不明顯，在斷裂帶等構(gòu)造復(fù)雜地區(qū)的變化較大。

2)利用逐日滑動算法進(jìn)行轉(zhuǎn)換函數(shù)時序分析，并確立最優(yōu)計算參數(shù)。分析顯示，轉(zhuǎn)換函數(shù)具有長周期變化背景，呈明顯的季節(jié)性波動，周期為1 a，轉(zhuǎn)換函數(shù)A的穩(wěn)定性優(yōu)于B，更適合進(jìn)行時移分析。

3)感應(yīng)矢量結(jié)果顯示，沂沭斷裂帶是一個顯著的電性分界面，其兩側(cè)的膠東隆起和魯西隆起具有明顯的電性結(jié)構(gòu)差異，魯西隆起區(qū)內(nèi)部的橫向電性結(jié)構(gòu)差異較小。山東省地下電性結(jié)構(gòu)整體橫向不均勻性未隨深度發(fā)生變化，淺層不均勻性較強。處于魯西隆起的濟(jì)寧地區(qū)在40～80 min周期內(nèi)電流聚集方向發(fā)生近180°轉(zhuǎn)換，說明在該深度下的NE向存在較大范圍的高導(dǎo)構(gòu)造區(qū)。

本文對山東省地磁轉(zhuǎn)換函數(shù)和感應(yīng)矢量開展背景分析和時序特征分析，由于山東地區(qū)震例較少，因此未開展震前異常信息的提取分析。后續(xù)可結(jié)合川滇等地區(qū)的資料進(jìn)行震前異常信息提取，從逐日滑動轉(zhuǎn)換函數(shù)的時移變化中挖掘孕震區(qū)電性結(jié)構(gòu)變化信息，加強對震前區(qū)域電性結(jié)構(gòu)變化的研究，為地震危險性評價提供依據(jù)。

致謝:本文頻譜計算使用的是馮志生研究員、朱培育高工編寫的軟件，在此表示感謝！