山東省地磁轉(zhuǎn)換函數(shù)和感應(yīng)矢量特征研究

管貽亮 張 玲 李希亮 董曉娜

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

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

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

1 數(shù)據(jù)來(lái)源與計(jì)算方法

1.1 數(shù)據(jù)來(lái)源

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

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

1.2 計(jì)算方法

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

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

(1)

式中，所有量均為復(fù)數(shù)；Z(ω)、X(ω)、Y(ω)分別為地磁垂直分量、南北分量和東西分量在某頻率ω上的頻譜；A、B為地磁轉(zhuǎn)換函數(shù)，其表現(xiàn)形式與MT中阻抗張量的估算方式一致。可使用最小二乘法使殘差趨于最小值，其殘差ε可表示為：

(2)

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

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

Gr=Ari+Brj,Gi=Aii+Bij

(3)

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

α=arctan(Br/Ar)，

(4)

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

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

2 結(jié)果分析

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

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

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

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

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

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

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

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

表1 不同窗長(zhǎng)計(jì)算的泰安臺(tái)和郯城臺(tái)Ar、Br的相關(guān)系數(shù)

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

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

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

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

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

3 結(jié) 語(yǔ)

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

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

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

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

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