基于自參考方法求取地磁測(cè)深響應(yīng)函數(shù)

張艷輝 楊 悅

1）中國(guó)河北050043 石家莊鐵道大學(xué)安全工程與應(yīng)急管理學(xué)院

2）中國(guó)青島266061 自然資源部第一海洋研究所

0 引言

地磁臺(tái)站所觀測(cè)到的磁場(chǎng)變化數(shù)據(jù)包含許多關(guān)于地幔過(guò)渡帶和下地幔電性結(jié)構(gòu)的信息。而通過(guò)對(duì)地球內(nèi)部電性結(jié)構(gòu)的研究可以探索地球的演化歷史、動(dòng)力學(xué)特征等問(wèn)題。

Banks（1969）首先提出將地磁數(shù)據(jù)轉(zhuǎn)換為C-響應(yīng)函數(shù)，進(jìn)而通過(guò)反演C-響應(yīng)函數(shù)來(lái)研究地球內(nèi)部電性結(jié)構(gòu)的分布。目前被廣泛應(yīng)用的C-響應(yīng)函數(shù)求取方法有2 種：①Z/H方法（H為水平磁場(chǎng)，由單一臺(tái)站提供）。該方法基于球諧函數(shù)分析求解，并被應(yīng)用在了研究美國(guó)西南地區(qū)的導(dǎo)電結(jié)構(gòu)上；②Z/Y方法，垂直分量由單一臺(tái)站提供，而水平分量由全球臺(tái)站數(shù)據(jù)的球諧系數(shù)表示（Olsen，1999）。

國(guó)內(nèi)學(xué)者對(duì)求取地磁測(cè)深C-響應(yīng)做了諸多研究，并取得一系列成果。如：杜興信等（1995）使用單臺(tái)Z/H方法，分析了陜西地區(qū)深部電導(dǎo)率；范國(guó)華等（1997）利用梯度法，分別求取了25 個(gè)臺(tái)站周期為24 h、12 h、8 h、6 h 的C-響應(yīng)值；張貴賓（1998）詳細(xì)介紹了地磁測(cè)深方法的原理，并采用Z/Y方法求取我國(guó)東北地區(qū)地磁臺(tái)站數(shù)據(jù)的響應(yīng)函數(shù)。因?yàn)榈卮判盘?hào)微弱，易受干擾（李雪華等，2012），所以地磁響應(yīng)曲線(xiàn)質(zhì)量較差。

BIRRP（Bounded Influence，Remote Reference Processing）是由Chave 等（2004）在Robust法基礎(chǔ)上開(kāi)發(fā)的數(shù)據(jù)處理軟件。Bounded Influence 方法將標(biāo)準(zhǔn)的M 估計(jì)方法和Hat Matrix diagonal 統(tǒng)計(jì)分析相結(jié)合，可以由時(shí)間域的電磁場(chǎng)信號(hào)通過(guò)遠(yuǎn)參考方法得到穩(wěn)健的頻率域響應(yīng)。該方法可以消除地磁場(chǎng)信號(hào)中的相關(guān)噪聲，同時(shí)具有較高效率。

為了更好地處理中國(guó)地區(qū)地磁臺(tái)站數(shù)據(jù)，將BIRRP 中使用的遠(yuǎn)參考方式改為基于本身的自參考方式，并以此求取地磁臺(tái)站C-響應(yīng)函數(shù)。分析驗(yàn)證結(jié)果顯示，本方法可以得到較為可靠的地磁測(cè)深C-響應(yīng)數(shù)據(jù)，為以后的科研工作奠定了基礎(chǔ)。

1 C-響應(yīng)求取的基本原理

在磁層源可以近似表示為的條件下，C-響應(yīng)計(jì)算公式（Banks，1969）可以表示為

其中，α0為地球的平均半徑，Hr、Hθ分別為地表處垂直指向地心和水平南向的磁場(chǎng)，θ為地磁余緯度，ω為角頻率。關(guān)于Hr和Hθ的求取可參考Semenov 等（2012）的研究。

2 數(shù)據(jù)處理基本技術(shù)

常規(guī)地磁C-響應(yīng)求取流程為：①原始時(shí)間序列中長(zhǎng)期變化場(chǎng)的去除；②C-響應(yīng)的求取需在地磁坐標(biāo)下進(jìn)行，因此進(jìn)行地磁分量時(shí)間序列的坐標(biāo)系旋轉(zhuǎn)；③時(shí)間序列轉(zhuǎn)到頻率域分析，一般使用離散傅里葉變換；④按照公式（1）進(jìn)行不同周期的C-響應(yīng)求??；⑤海洋效應(yīng)影響去除。

為了得到準(zhǔn)確、穩(wěn)定的C-響應(yīng)曲線(xiàn)，在每一個(gè)步驟均采取不同技術(shù)手段，以保證處理過(guò)程的正確性。

此外，高緯度地磁臺(tái)站數(shù)據(jù)易受極光效應(yīng)影響。目前的研究并不能較好地消除極光效應(yīng)影響（Semenov et al，2012）。然而，由于使用的地磁臺(tái)站均位于10°—50°N 范圍，受極光效應(yīng)影響較小，故極光效應(yīng)不予考慮。

2.1 地磁分量求取

國(guó)家地磁臺(tái)網(wǎng)中心地磁數(shù)據(jù)存儲(chǔ)采用IAGA（The International Association of Geomagnetism and Aeronomy，國(guó)際地磁學(xué)和高層大氣物理協(xié)會(huì)）格式，包含地磁偏角D、水平分量H和垂直分量Z的絕對(duì)記錄值或相對(duì)記錄值，在實(shí)際使用時(shí)一般需要將數(shù)據(jù)轉(zhuǎn)換為水平分量X、Y和垂直分量Z（圖1）。而在C-響應(yīng)求取中，需要使用水平南向磁場(chǎng)Hθ和垂直地心的磁場(chǎng)分量Hr（即Z）。將水平分量H分解得到水平南向分量Hθ，即

圖1 蘭州臺(tái)站地磁三分量時(shí)間序列Fig.1 Time series of the three components of the geomagnetic data of Lanzhou station

2.2 長(zhǎng)期變化場(chǎng)去除

地磁測(cè)深研究的是磁層電流（外部場(chǎng)源）產(chǎn)生的變化磁場(chǎng)。地磁場(chǎng)的主要組成部分是起源于地球內(nèi)部的長(zhǎng)期變化場(chǎng)。在進(jìn)行較長(zhǎng)周期的電磁感應(yīng)研究中，需要將長(zhǎng)期變化場(chǎng)去除（Olsen，1999）。國(guó)際地磁參考場(chǎng)（IGRF）是描述地球主磁場(chǎng)和長(zhǎng)期變化場(chǎng)的數(shù)學(xué)模型。在地球表層和上部區(qū)域，主磁場(chǎng)（地球內(nèi)部場(chǎng)源引起的）各分量可以由地磁位函數(shù)V的梯度來(lái)表示

其中，φ和θ分別是經(jīng)度和余緯度，α=6 371.2 km，r是地心距，N表示最高階數(shù)，和是地磁位的球諧系數(shù)，是Schmidt 歸一化締合Legendre 函數(shù)。

IAGA 提供了計(jì)算地磁場(chǎng)不同分量長(zhǎng)期變化的程序代碼。該程序可用于去除各臺(tái)站不同分量對(duì)應(yīng)的長(zhǎng)期變化場(chǎng)（https://www.ngdc.noaa.gov/IAGA/vmod/igrf12.f）。

為了驗(yàn)證長(zhǎng)期變化場(chǎng)去除與否對(duì)C-響應(yīng)求取結(jié)果的影響，從國(guó)家地磁臺(tái)網(wǎng)中心獲取部分臺(tái)站的地磁數(shù)據(jù)，并將長(zhǎng)春臺(tái)（CNH）和滿(mǎn)洲里臺(tái)（MZL）去除與未去除長(zhǎng)期變化場(chǎng)的C-響應(yīng)計(jì)算結(jié)果進(jìn)行對(duì)比（圖2）。由圖2 可見(jiàn)，長(zhǎng)春臺(tái)站數(shù)據(jù)長(zhǎng)期場(chǎng)去除與否對(duì)短周期的C-響應(yīng)影響不大，二者幾乎重合，但在長(zhǎng)周期段，二者存在明顯差別；滿(mǎn)洲里臺(tái)站數(shù)據(jù)的計(jì)算結(jié)果顯示，二者的響應(yīng)曲線(xiàn)趨勢(shì)類(lèi)似，然而具體數(shù)值存在明顯區(qū)別。因此，在C-響應(yīng)求取過(guò)程中，有必要去除長(zhǎng)期變化場(chǎng)。

圖2 長(zhǎng)春臺(tái)和滿(mǎn)洲里臺(tái)長(zhǎng)期場(chǎng)去除與否的C-響應(yīng)處理結(jié)果對(duì)比Fig.2 Comparison of C-response results of Changchun (CNH) and Manzhouli (MZL)stations with or without secular variation removal

2.3 坐標(biāo)系轉(zhuǎn)換

對(duì)已知地球模型進(jìn)行數(shù)值模擬時(shí)，通常假設(shè)源為地球磁層中的環(huán)形電流，并以球諧函數(shù)來(lái)表示，該電流環(huán)與地球的磁赤道同心。地磁測(cè)深數(shù)值模擬和反演均在地磁坐標(biāo)系中進(jìn)行。因此，數(shù)據(jù)處理時(shí)需要將地理坐標(biāo)系中的磁場(chǎng)數(shù)據(jù)旋轉(zhuǎn)到地磁坐標(biāo)系中，以得到地磁坐標(biāo)系中的C-響應(yīng)。

坐標(biāo)系旋轉(zhuǎn)的基本思路為，通過(guò)國(guó)際地磁參考場(chǎng)（IGRF）移除長(zhǎng)期變化場(chǎng)，利用IGRF 計(jì)算程序得到地磁坐標(biāo)系中某點(diǎn)的磁偏角D，由去除長(zhǎng)期場(chǎng)的Hx和Hy可得到去除長(zhǎng)期場(chǎng)的H及其與地理北的夾角，通過(guò)簡(jiǎn)單加減法即可得到該夾角，也就得到了地磁北和東分量。

圖3 為坐標(biāo)系轉(zhuǎn)換示意圖，O為臺(tái)站位置，H為臺(tái)站觀測(cè)水平磁場(chǎng)強(qiáng)度，Dobs為臺(tái)站觀測(cè)磁偏角（并不一定等于該點(diǎn)的磁坐標(biāo)系與地理坐標(biāo)系的夾角，因?yàn)樵擖c(diǎn)的水平磁場(chǎng)不一定指向磁北，但在一維地球和源的條件下，取等是成立的），Hx和Hy為臺(tái)站觀測(cè)的地理北向和地理東向分量，Hx-sv、Hy-sv是減去長(zhǎng)期變化場(chǎng)（由國(guó)際地磁參考場(chǎng)IGRF 計(jì)算得到）的地理北向和地理東向分量，H-sv是減去長(zhǎng)期場(chǎng)之后的水平磁場(chǎng)強(qiáng)度，DIGRF是由IGRF 計(jì)算得到的磁偏角，可以認(rèn)為是地磁北和地理北方向的夾角。

圖3 地理坐標(biāo)系旋轉(zhuǎn)到地磁坐標(biāo)系示意Fig.3 Schematic diagram of rotation from geographical coordinate system to geomagnetic coordinate system

2.4 時(shí)頻轉(zhuǎn)換

在進(jìn)行地磁測(cè)深C-響應(yīng)求取時(shí)，使用一個(gè)短的自回歸濾波器對(duì)每一個(gè)時(shí)間序列進(jìn)行預(yù)白化，可以有效減小頻譜偏差。將數(shù)據(jù)分段，截取每個(gè)頻率對(duì)應(yīng)的時(shí)間段。由于常規(guī)錐形截取時(shí)間序列會(huì)造成頻譜泄漏和影響分辨率，故采用Slepian 序列窗來(lái)進(jìn)行數(shù)據(jù)截?。–have et al，2004），可以通過(guò)調(diào)節(jié)時(shí)間帶寬來(lái)滿(mǎn)足不同的周期要求。

由于需要在頻率域分析地磁響應(yīng)，故使用離散傅里葉變換將截取的時(shí)間序列轉(zhuǎn)到頻率域，并采用重疊平均法來(lái)增加計(jì)算的可靠性（Chave et al，2004）。

2.5 基于Bounded Influence 方法的C-響應(yīng)處理

采用Bounded Influence（BI）方法計(jì)算地磁測(cè)深的C-響應(yīng)函數(shù)。Semenov 等（2012）通過(guò)對(duì)比，驗(yàn)證了在地磁測(cè)深數(shù)據(jù)處理中，采用臺(tái)站自參考方法與遠(yuǎn)參考具有良好的一致性。因此，采用自參考方法進(jìn)行C-響應(yīng)的求取。

其中：H*為H的復(fù)共軛；*表示矩陣的復(fù)共軛；α=6 371.2 km，為地球半徑；為疊加自（互）功率譜（Huber，1981），表達(dá)式為

式中，A和B代表不同的磁場(chǎng)分量或者參考分量。

為了判斷不同質(zhì)量數(shù)據(jù)的臺(tái)站C-響應(yīng)的可靠性，通過(guò)求取平方相關(guān)系數(shù)和數(shù)據(jù)誤差來(lái)進(jìn)行C-響應(yīng)質(zhì)量評(píng)價(jià)。在實(shí)際數(shù)據(jù)處理中，可以去除相關(guān)系數(shù)較低或數(shù)據(jù)誤差較大的響應(yīng)數(shù)據(jù)。

在計(jì)算得到C-響應(yīng)之后，本文按照Semenov 等（2012）的方法計(jì)算了C-響應(yīng)不同周期的平方相關(guān)系數(shù)。

其中，H*、Z*分別表示H和Z的復(fù)共軛。C-響應(yīng)的數(shù)據(jù)誤差計(jì)算式（Schmucker，1999）為

其中，1-β是置信水平，本文的置信水平設(shè)置為0.9（Chave et al，2004）；L是時(shí)間序列的分段數(shù)量，周期越長(zhǎng)，L越小。

3 臺(tái)站數(shù)據(jù)處理結(jié)果

我們利用上述方法對(duì)ASP（Alice Springs）和HON（Honolulu）臺(tái)站數(shù)據(jù)（數(shù)據(jù)下載地址：http://geomag.org/）求取C-響應(yīng)并與已有結(jié)果進(jìn)行了對(duì)比，見(jiàn)圖4。可以看到本文方法處理得到的響應(yīng)曲線(xiàn)在短周期段與Semenov 等（2012）的結(jié)果保持著較好的一致性；在長(zhǎng)周期段更加穩(wěn)定，“飛點(diǎn)”情況明顯減少。從平方相關(guān)系數(shù)和誤差棒也能明顯看到本方法求得的C-響應(yīng)相關(guān)性系數(shù)更高，誤差更小，這說(shuō)明了該方法的正確性和穩(wěn)定性。

圖4 本方法求取的C-響應(yīng)與Semenov 等（2012）結(jié)果對(duì)比三角形為Semenov 的結(jié)果；圓圈為本文的結(jié)果Fig.3 Comparison of the C-response obtained by this method with the result of Semenov et al (2012)

為了進(jìn)一步檢驗(yàn)本方法求取的地磁臺(tái)站C-響應(yīng)的正確性，利用從中國(guó)國(guó)家地磁臺(tái)網(wǎng)中心獲取的北京和蘭州地磁臺(tái)站記錄數(shù)據(jù)求取C-響應(yīng)，所得結(jié)果見(jiàn)圖5。由圖5 可見(jiàn)，2個(gè)臺(tái)站數(shù)據(jù)的響應(yīng)曲線(xiàn)穩(wěn)定，且在短周期段具有較低的數(shù)據(jù)誤差和較高的一致性；而在長(zhǎng)周期段，受限于數(shù)據(jù)記錄長(zhǎng)度，無(wú)論是平方一致性曲線(xiàn)還是數(shù)據(jù)誤差以及C-響應(yīng)曲線(xiàn)的連續(xù)性均相對(duì)較差。整體而言，利用該方法可以為中國(guó)地磁臺(tái)站C-響應(yīng)求取提供良好的技術(shù)支持。

圖5 使用本方法求取的北京（BJI）臺(tái)站和蘭州（LZH）臺(tái)站C-響應(yīng)和平方一致性曲線(xiàn)Fig.5 C-responses and square coherences curve of Beijing (BJI) and Lanzhou(LZH) stations obtained by the method in this paper

4 結(jié)論

利用自參考Bounded Influence 方法求取地磁測(cè)深C-響應(yīng)值可有效減少實(shí)測(cè)地磁場(chǎng)數(shù)據(jù)的環(huán)境干擾和噪聲影響。C-響應(yīng)曲線(xiàn)的“ 飛點(diǎn)”明顯減少，同時(shí)，響應(yīng)曲線(xiàn)的平方相關(guān)系數(shù)更高、誤差更小。本方法為高分辨率三維地磁測(cè)深反演研究奠定了基礎(chǔ)。