基于自參考方法求取地磁測深響應函數

張艷輝 楊 悅

1）中國河北050043 石家莊鐵道大學安全工程與應急管理學院

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

0 引言

地磁臺站所觀測到的磁場變化數據包含許多關于地幔過渡帶和下地幔電性結構的信息。而通過對地球內部電性結構的研究可以探索地球的演化歷史、動力學特征等問題。

Banks（1969）首先提出將地磁數據轉換為C-響應函數，進而通過反演C-響應函數來研究地球內部電性結構的分布。目前被廣泛應用的C-響應函數求取方法有2 種：①Z/H方法（H為水平磁場，由單一臺站提供）。該方法基于球諧函數分析求解，并被應用在了研究美國西南地區的導電結構上；②Z/Y方法，垂直分量由單一臺站提供，而水平分量由全球臺站數據的球諧系數表示（Olsen，1999）。

國內學者對求取地磁測深C-響應做了諸多研究，并取得一系列成果。如：杜興信等（1995）使用單臺Z/H方法，分析了陜西地區深部電導率；范國華等（1997）利用梯度法，分別求取了25 個臺站周期為24 h、12 h、8 h、6 h 的C-響應值；張貴賓（1998）詳細介紹了地磁測深方法的原理，并采用Z/Y方法求取我國東北地區地磁臺站數據的響應函數。因為地磁信號微弱，易受干擾（李雪華等，2012），所以地磁響應曲線質量較差。

BIRRP（Bounded Influence，Remote Reference Processing）是由Chave 等（2004）在Robust法基礎上開發的數據處理軟件。Bounded Influence 方法將標準的M 估計方法和Hat Matrix diagonal 統計分析相結合，可以由時間域的電磁場信號通過遠參考方法得到穩健的頻率域響應。該方法可以消除地磁場信號中的相關噪聲，同時具有較高效率。

為了更好地處理中國地區地磁臺站數據，將BIRRP 中使用的遠參考方式改為基于本身的自參考方式，并以此求取地磁臺站C-響應函數。分析驗證結果顯示，本方法可以得到較為可靠的地磁測深C-響應數據，為以后的科研工作奠定了基礎。

1 C-響應求取的基本原理

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

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

2 數據處理基本技術

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

為了得到準確、穩定的C-響應曲線，在每一個步驟均采取不同技術手段，以保證處理過程的正確性。

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

2.1 地磁分量求取

國家地磁臺網中心地磁數據存儲采用IAGA（The International Association of Geomagnetism and Aeronomy，國際地磁學和高層大氣物理協會）格式，包含地磁偏角D、水平分量H和垂直分量Z的絕對記錄值或相對記錄值，在實際使用時一般需要將數據轉換為水平分量X、Y和垂直分量Z（圖1）。而在C-響應求取中，需要使用水平南向磁場Hθ和垂直地心的磁場分量Hr（即Z）。將水平分量H分解得到水平南向分量Hθ，即

圖1 蘭州臺站地磁三分量時間序列Fig.1 Time series of the three components of the geomagnetic data of Lanzhou station

2.2 長期變化場去除

地磁測深研究的是磁層電流（外部場源）產生的變化磁場。地磁場的主要組成部分是起源于地球內部的長期變化場。在進行較長周期的電磁感應研究中，需要將長期變化場去除（Olsen，1999）。國際地磁參考場（IGRF）是描述地球主磁場和長期變化場的數學模型。在地球表層和上部區域，主磁場（地球內部場源引起的）各分量可以由地磁位函數V的梯度來表示

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

IAGA 提供了計算地磁場不同分量長期變化的程序代碼。該程序可用于去除各臺站不同分量對應的長期變化場（https://www.ngdc.noaa.gov/IAGA/vmod/igrf12.f）。

為了驗證長期變化場去除與否對C-響應求取結果的影響，從國家地磁臺網中心獲取部分臺站的地磁數據，并將長春臺（CNH）和滿洲里臺（MZL）去除與未去除長期變化場的C-響應計算結果進行對比（圖2）。由圖2 可見，長春臺站數據長期場去除與否對短周期的C-響應影響不大，二者幾乎重合，但在長周期段，二者存在明顯差別；滿洲里臺站數據的計算結果顯示，二者的響應曲線趨勢類似，然而具體數值存在明顯區別。因此，在C-響應求取過程中，有必要去除長期變化場。

圖2 長春臺和滿洲里臺長期場去除與否的C-響應處理結果對比Fig.2 Comparison of C-response results of Changchun (CNH) and Manzhouli (MZL)stations with or without secular variation removal

2.3 坐標系轉換

對已知地球模型進行數值模擬時，通常假設源為地球磁層中的環形電流，并以球諧函數來表示，該電流環與地球的磁赤道同心。地磁測深數值模擬和反演均在地磁坐標系中進行。因此，數據處理時需要將地理坐標系中的磁場數據旋轉到地磁坐標系中，以得到地磁坐標系中的C-響應。

坐標系旋轉的基本思路為，通過國際地磁參考場（IGRF）移除長期變化場，利用IGRF 計算程序得到地磁坐標系中某點的磁偏角D，由去除長期場的Hx和Hy可得到去除長期場的H及其與地理北的夾角，通過簡單加減法即可得到該夾角，也就得到了地磁北和東分量。

圖3 為坐標系轉換示意圖，O為臺站位置，H為臺站觀測水平磁場強度，Dobs為臺站觀測磁偏角（并不一定等于該點的磁坐標系與地理坐標系的夾角，因為該點的水平磁場不一定指向磁北，但在一維地球和源的條件下，取等是成立的），Hx和Hy為臺站觀測的地理北向和地理東向分量，Hx-sv、Hy-sv是減去長期變化場（由國際地磁參考場IGRF 計算得到）的地理北向和地理東向分量，H-sv是減去長期場之后的水平磁場強度，DIGRF是由IGRF 計算得到的磁偏角，可以認為是地磁北和地理北方向的夾角。

圖3 地理坐標系旋轉到地磁坐標系示意Fig.3 Schematic diagram of rotation from geographical coordinate system to geomagnetic coordinate system

2.4 時頻轉換

在進行地磁測深C-響應求取時，使用一個短的自回歸濾波器對每一個時間序列進行預白化，可以有效減小頻譜偏差。將數據分段，截取每個頻率對應的時間段。由于常規錐形截取時間序列會造成頻譜泄漏和影響分辨率，故采用Slepian 序列窗來進行數據截取（Chave et al，2004），可以通過調節時間帶寬來滿足不同的周期要求。

由于需要在頻率域分析地磁響應，故使用離散傅里葉變換將截取的時間序列轉到頻率域，并采用重疊平均法來增加計算的可靠性（Chave et al，2004）。

2.5 基于Bounded Influence 方法的C-響應處理

采用Bounded Influence（BI）方法計算地磁測深的C-響應函數。Semenov 等（2012）通過對比，驗證了在地磁測深數據處理中，采用臺站自參考方法與遠參考具有良好的一致性。因此，采用自參考方法進行C-響應的求取。

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

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

為了判斷不同質量數據的臺站C-響應的可靠性，通過求取平方相關系數和數據誤差來進行C-響應質量評價。在實際數據處理中，可以去除相關系數較低或數據誤差較大的響應數據。

在計算得到C-響應之后，本文按照Semenov 等（2012）的方法計算了C-響應不同周期的平方相關系數。

其中，H*、Z*分別表示H和Z的復共軛。C-響應的數據誤差計算式（Schmucker，1999）為

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

3 臺站數據處理結果

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

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

為了進一步檢驗本方法求取的地磁臺站C-響應的正確性，利用從中國國家地磁臺網中心獲取的北京和蘭州地磁臺站記錄數據求取C-響應，所得結果見圖5。由圖5 可見，2個臺站數據的響應曲線穩定，且在短周期段具有較低的數據誤差和較高的一致性；而在長周期段，受限于數據記錄長度，無論是平方一致性曲線還是數據誤差以及C-響應曲線的連續性均相對較差。整體而言，利用該方法可以為中國地磁臺站C-響應求取提供良好的技術支持。

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

4 結論

利用自參考Bounded Influence 方法求取地磁測深C-響應值可有效減少實測地磁場數據的環境干擾和噪聲影響。C-響應曲線的“ 飛點”明顯減少，同時，響應曲線的平方相關系數更高、誤差更小。本方法為高分辨率三維地磁測深反演研究奠定了基礎。