大地電磁多參考站陣列數據處理方法

周 聰 湯井田 原 源 鄧居智 石福升 李 勇

(①核資源與環境國家重點實驗室(東華理工大學)，江西南昌 330013； ②自然資源部地球物理電磁法探測技術重點實驗室，河北廊坊 065000； ③有色金屬成礦預測與地質環境監測教育部重點實驗室(中南大學)，湖南長沙 410083；④自然資源部覆蓋區深部資源勘查工程技術創新中心，安徽合肥 230001)

0 引言

大地電磁法(Magnetotelluric，MT)是一類重要的地球物理方法[1-2]。隨著經濟的發展和應用環境的多樣化，MT方法所面臨的噪聲問題越來越復雜。如何有效壓制各種電磁噪聲，獲得高質量的阻抗估計，一直是MT方法最重要的任務之一[1,3-5]。

遠參考法是大地電磁法壓制噪聲最有效的策略之一。Goubau等[6]首先提出互功率譜估計法，以壓制常規單站處理阻抗估計中自功率譜計算所引起的噪聲。在此基礎上，Gamble等[7]提出了遠參考估計法，該方法要求在距離測點一定范圍內觀測磁場信號的變化，利用遠參考站與測站間信號相干、噪聲不相干的性質，以互功率譜代替自功率譜計算，以數據相干性壓制噪聲的影響。隨著人文噪聲的影響日益嚴重，遠參考法得到了極大的發展并逐步普及[8-12]。

遠參考技術的發展要求參考數據規模不斷地增大；數據的利用從雙磁道向電磁多測道、多站疊加及多站處理不斷發展。常規參考道數據一般選擇兩道水平磁場，這是因為磁場傳播過程相對更穩定，不同測站間磁場的相關性更好。但當參考站的磁場數據質量不佳時，電場數據也可被用于作為參考道進行互功率譜的計算，即所謂“電參考”[13]。Epishkin[14]提出一種利用所有測道的遠參考技術，可同時利用電場和磁場分別計算一套權重系數，提高了參考道數據的利用率，進而提高了MT數據處理質量。但該算法僅限于單參考站的處理。Varentsov等[15]提出利用遠參考站的磁場與測站磁場間的相干性分析，計算一系列閾值參數，據此對測站的含噪功率譜數據進行挑選，從而提高阻抗估計的穩健性，即所謂“遠參考磁場控制法(RRMC)”。張剛等[16]將這一思想進一步發展到多點遠參考數據，但其實質仍是不同參考站處理結果的疊加，因而對于持續性噪聲，仍難以通過功率譜數據的挑選提高數據質量。MT多站疊加技術、“偽遠參”及站間依賴關系等方法[17-19]拓展了遠參考法的應用，但這些方法未明確多參考站的應用方式。Larsen等[20]提出基于遠參考站的信噪分量方法，Oettinger等[21]在此基礎上進一步提出一種可利用雙遠參考數據壓制相干噪聲的技術，并可同時獲得相干噪聲的阻抗信息。該方法為雙遠參考站數據的合理利用提供了可行策略，但當遠參考站超過兩個時，該方法無法進一步提高數據利用率。

進一步提高參考數據的利用率是遠參考技術的發展方向之一。隨著陣列觀測技術的不斷發展[22]，多參考站數據的獲取已具備實際可操作性。多參考站數據不僅可以從測區外直接進行選點觀測，還可以從測區內的低噪數據中選擇。例如，從測區內低噪數據中選擇參考站數據作互參考處理[23]就是常用的一種選擇。一些學者研究了多個遠參考站條件下不同參考距離對處理結果的影響[24-26]，然而算法的實質仍是各參考站數據相互獨立地進行處理，未能充分利用多點遠參考數據間的相干信息。顯而易見，多點遠參考站間均可能滿足信號相干、噪聲不相干的條件。如何更加有效地利用多參考站數據，統一提取多站、多道間的相干信息，增加數據利用率，提高處理質量，仍是重要的研究方向。

陣列電磁數據處理方法的提出與發展為多參考站數據的利用提供了基礎[27-31]。常規陣列處理僅構建單一數據陣列[28]，天然場源極化參數的提取精度易受影響。本文在常規陣列電磁數據處理方法的基礎上，發展了多參考站條件下的陣列電磁數據處理技術，提出分別構建天然場占優的數據陣列和含噪的目標數據陣列，以陣列數據分析方法提取高精度天然電磁場極化參數，繼而提高目標數據陣列的天然場響應估計質量。實測數據試驗驗證了該方法的有效性。

1 陣列電磁數據處理模型

對于頻域電磁場，假設在I個激勵時窗內，系統輸入端共有L個電磁場源同時進行激勵，輸出端共有K個測道同步測量系統對所有輸入激勵的總響應。第l(l=1,2,…,L)個場源在第i(i=1,2,…,I)個時窗內的極化參數記為Ali，L×I個極化參數構成輸入陣列A；第k(k=1,2,…,K)個測道在第i個時窗的觀測數據記為Xki，K×I個觀測數據構成輸出陣列X。線性時不變系統條件下，系統頻率特性響應Ukl不隨時間變化，輸出與輸入滿足線性關系[27,30]

(1)

矩陣形式為

X=UA+R

(2)

式中：R為頻域輸出端噪聲矩陣，可存在于各個測道與時窗中；A為場源的極化信息，或稱為極化參數矩陣；K×L維矩陣U反映了地球電磁系統本身的特性，以及激勵源與道的空間信息，與激勵延時(觀測時窗)無關，稱為系統響應矩陣。式(1)和式(2)即為陣列電磁勘探模型的頻率域輸出—輸入關系表達式。

式(1)和式(2)亦是大地電磁法陣列數據處理的基本模型。通過合理的實施方案，可獲得陣列觀測數據。求解上述陣列方程組，并從中提取相應的地電解釋參數，可以建立陣列數據處理的基本理論。

2 多參考站陣列數據處理技術

由于天然場源及噪聲的隨機性和不可預測性，很難準確地知道場源個數，因此采用數據降維類算法(如主成分分析等)求解式(2)是合理的。一般而言，在大地電磁法和音頻大地電磁法涉及的觀測頻段范圍內，天然場源信號可被認為是均勻平面電磁波，因其發收距遠大于觀測尺度，不同方向場源的主要能量可以通過分解、合并等效為兩個正交水平場源的作用，或者說天然場源的主要極化可等效至兩個正交水平方向上。

根據上述特點，Egbert[27]設計了基于穩健主成分分析的常規陣列電磁數據處理方案。該方案將遠參考與目標測站置于同一陣列，在低噪環境下，可獲得較好的處理結果。然而，當環境中存在相干噪聲時，含噪數據會嚴重影響阻抗數據的估計質量。

多點遠參考數據為提高天然場源極化參數的估計質量提供了可能。利用多點遠參考數據，將遠參考測站與目標測站物理分離，分別構建兩類陣列數據集，即天然場信號占主導的低噪參考陣列數據集和含噪的目標陣列數據集。先提取高精度的天然場源極化參數，再利用該參數提高目標陣列數據集的處理質量。在此思路下，本文提出一種分步求解方案，即先估計場源極化參數，后求系統響應參數。

2.1 優選參考站數據，分別構建參考數據陣列和目標測站陣列

一般可通過在測區外布設兩個甚至多個遠參考站獲得多點遠參考數據。當測區范圍較大、噪聲分布不均勻時，還可優選測區內的同步低噪測站數據，或在低噪時段中加入參考數據陣列。參考站數據的選擇與常規單站遠參考、互參考數據的挑選方法相同，以阻抗估計誤差、信號相干度和極化方向等多種參數進行信、噪識別，挑選低噪測站和時段，需保證參考站本身的數據品質較高，并且盡可能與待處理的測站噪聲不相干。

實際中，根據參考測站的測道組合方式，主要有以下幾種參考陣列組合： ①經過挑選的磁道組合；②經過挑選的電道組合； ③經過挑選的混合測道組合，包括電道、磁道甚至其他可能的梯度測道等。記參考數據矩陣為Xr，對應上述三種組合方式可分別寫為

(3)

(4)

(5)

式中：電磁場分量的上角標1、2、… 、J表示經過挑選的參考測站序號，J為參考測站總數；E表示電場；H表示磁場。

同樣地，可以將待處理目標測站的張量電道、磁道數據寫入同一個數據矩陣中，記為Xn。

(2)利用參考數據陣列，建立陣列方程組

Xr=UrA+Rr

(6)

式中Rr為參考陣列的信號殘差。

(3)根據多個目標測站的數據可建立陣列方程組

Xn=UnA+Rn

(7)

式中Rn為目標陣列的信號殘差，高噪條件下，該殘差的方差較大。

2.2 估計天然場源極化參數

通過合適的數據選擇，可確保參考數據陣列中天然場占主導。如前文所述，天然場源的主要極化可等效至兩個正交水平方向上，即L≈2。此時遵循Egbert[27]和Smirnov等[29]提出的穩健主成分分析(Robust Principal Component Analysis,RPCA)方法，可提取出數據中的主成分信息

(8)

穩健估計式(7)的實現可直接采用現有的RPCA算法，如公開的Libra統計學代碼[31]。但經測試，該算法在數據量較大時，對計算資源要求較高，計算速度慢。因此，本文采用常規陣列處理中所用的RPCA算法[27,29]，其簡要流程如下。

(1)計算K個測道的不相干噪聲方差，構建K×K維對角矩陣CK。

(4)設置最大循環次數jmax、迭代終止條件閾值ε，對j(j=0,1,…,jmax)作循環計算

作SVD，Y(j)W(j)=u(j+1)s(j+1)[v(j+1)]T

如果‖u(j+1)s(j+1)s(j+1)[u(j+1)]*-I‖<ε,停止迭代； 其中：Y、r、w、W分別為預測矩陣、誤差、權系數及權系數矩陣； diag(·)表示求對角矩陣。

2.3 估計目標測站的阻抗張量

利用最小二乘法，對目標數據矩陣Xn及式(7)進行帶參考信號的穩健回歸分析，可求得系統響應

(9)

式中W′為對式(5)進行分析時的穩健估計權值矩陣。

獲得Un的估計后，再利用陣列響應參數的提取方法計算阻抗及其視電阻率和相位[30]。阻抗的計算公式為

ZMT=[UE(UH)T][UH(UH)T]-1

(10)

式中UE和UH分別為電場和磁場張量矩陣。上式省略了測站的角標編號。

分析以上過程可知，與常規方法不同，本文將多參考站數據集中于統一的陣列數據中進行分析處理，而不是簡單的多點參考數據疊加。由于以上陣列數據分析的引入，提高了多參考站數據的利用率和天然場信息提取精度。基于噪聲的統計特征，理論上參與計算的數據越多，噪聲壓制效果越好。

3 應用實例

3.1 實測多參考站觀測試驗

2016年10月，中南大學在中國東部地區開展了多點遠參考觀測及處理試驗，測點分布如圖1所示。共包含目標測站1個(S1)，位于安徽省霍山地區，該處背景噪聲較高；同步在背景噪聲相對較低的地區布設參3個參考站，分別位于安徽省滁州市定遠縣(R1)、安徽省黃山市(R2)、江西省景德鎮市(R3)，與S1的距離分別約為160、200、280km；各測站的張量布設方向均為正南或正北，以北向為x方向。R1、R2參考站觀測水平4分量電磁場(Ex、Ey、Hx、Hy)，R3僅觀測水平2分量磁場(Hx、Hy)。數據采集儀器采用加拿大鳳凰公司生產的V5-2000型大地電磁儀，采樣率為15Hz，同步采集時長約18h。

圖1 試驗測點位置分布示意圖

3.2 觀測數據分析

陣列數據處理的理論依據之一是認為不同空間測站的觀測數據中包含大量相干信息，并可從中提取出相應的場源極化參數。為此，分析了實例中同步陣列觀測的時間域、時頻譜數據，如圖2所示。不難看出，S1、R1、R2和R3盡管相距數百千米，所處的噪聲環境各不相同，但各測站的同步電磁場間時間序列的總體輪廓及細節變化均具有相似性，功率譜隨時間的變化趨勢也一致，只是強度存在差異。這說明電磁數據處理可考慮空間多測站間的同步信息，進行多站聯合處理。

圖2 同步陣列觀測的Hy時間序列(左)和時頻譜(右)片段

為進一步分析各測站的含噪情況，進行了不同頻率的極化方向對比分析，如圖3所示。可以看出，S1測站1Hz和0.1Hz的極化方向呈集中趨勢，并且在整個觀測時段趨勢一致，說明該測站處包含穩定的人文噪聲場，且影響貫穿整個觀測時段。由于噪聲超過了穩健估計的“崩潰點”[32]，常規穩健估計方法很難得到可靠的處理結果。并且，因噪聲包含在所有觀測時段中，數據挑選類方法往往也難以湊效。各參考測站的極化方向分布相對均勻，僅R2測站略呈集中趨勢。考慮到0.1～1Hz處于MT的“死頻帶”范圍[33]，結合后文的視電阻率、相位數據分析可知，各參考測站在0.1～1Hz范圍內未包含明顯的人文噪聲，但天然場信號強度較低，分散的極化方向可能為天然場和隨機噪聲的混合表現。

圖3 同步陣列觀測在1 Hz和0.1Hz時的極化方向對比

3.3 處理結果對比

經過時頻轉換，先對張量觀測的S1、R1和R2測站進行了常規Robust處理和常規陣列處理，獲得了視電阻率和相位結果(圖4)。常規陣列處理中R3的兩道磁場也參與了計算，只是因該點未觀測電場，因此未獲得視電阻率、相位數據。可以看出，S1測站在0.1～1Hz的MT“死頻帶”范圍內出現了明顯的畸變，視電阻率曲線在低于0.1Hz頻段明顯脫節，yx模式的數據表現尤為明顯。而參考站R1和R2在0.1Hz附近頻段也出現了一定程度的畸變，這正是MT“死頻帶”畸變的顯著表現。由于“死頻帶”內天然場信號強度極低，甚至可能低于儀器的本底噪聲，因此即使在人文噪聲較低的測區，也難以獲得高質量的天然場響應。這也是目前MT數據處理中最棘手的問題之一。經過常規陣列處理后，利用多測道數據壓制了隨機噪聲的干擾，S1測站在0.1～1Hz內的形態得到明顯改善，但yx模式視電阻率曲線在0.1Hz前后頻段仍表現出一定的脫節，顯示出相干噪聲的影響仍然存在。這說明在陣列規模不大的條件下，常規陣列處理可較好地壓制隨機噪聲的影響，但對于相干噪聲，其處理效果不容樂觀。

圖4 單站Robust處理與常規陣列處理的視電阻率(上)和相位(下)曲線

基于本文所述的多參考站陣列數據處理方法，對測站S1進行了處理。結合前述算法可知，多參考站陣列處理與常規陣列處理的區別主要在于天然場源極化參數的估計策略不同。圖5是兩種算法估計得到的極化參數。由于極化參數為復數，圖中顯示的是其絕對值。可以看出，盡管兩種算法得到的天然場源極化參數變化規律基本一致，但細節仍存在諸多不同，說明了本文方法有別于常規陣列處理，取得的處理結果也不同。進一步地，圖6給出了多參考站陣列處理的阻抗視電阻率計算結果。對比可見，常規陣列處理的S1點(圖4d上)yx模式視電阻率曲線的相對高(0.2～2Hz)、低(0.01～0.6Hz)兩個頻段出現脫節。而多參考站陣列處理(圖6g)曲線在全頻段連續，顯然更合理。由大地電磁響應曲線連續的特征[1]可判斷，圖4d上中0.2～2Hz頻段的曲線為畸變形態，包含相了干噪聲的影響。這說明多參考站陣列處理可取得優于常規陣列處理的結果。

圖5 不同處理估計得到的天然場源極化參數x方向(a)和y方向(b)絕對值對比

圖6 測站S1采用不同參考站信息處理的視電阻率曲線對比

圖6中還給出了利用常規遠參考方法處理所得的結果，這里采用了不同的參考道組合。對比利用3個遠參考站磁道分別對S1進行常規遠參考處理的結果(圖6g～圖6i)與圖4所示的S1常規處理結果，可以發現3個遠參考站數據對S1在0.1～1Hz的畸變均有一定程度的改善，人文噪聲的畸變影響得到了壓制。相對而言，R1的參考效果相對最差，視電阻率在多個頻點仍有較強的畸變；R2的參考效果略好，但yx模式的誤差棒最大。R3的參考效果最優，僅在0.8Hz、 0.2Hz左右幾個頻點處出現了跳變。這也說明，實際中，通過優選遠參考可以得到更好的處理效果。對比單參考處理與兩兩組合的多參考處理，可以發現，R1、R2組合多參考處理結果(圖6d)優于他們各自的單參考處理(圖6a、圖6b))，前者不僅減少了畸變的頻點，也壓制了誤差棒； R1、R3的多參考處理結果(圖6e)和R2、R3的多參考處理結果(圖6f)分別優于R1的單參考處理(圖6a)和R2的單參考處理(圖6b)，前者均有效壓制了0.8Hz處R3單參考處理結果里的跳變點(圖6c)。整體上，3個參考站的磁道組合多參考結果(圖6g)僅在0.2Hz頻點處出現了跳變，取得了最佳處理效果。

圖6h和圖6i分別是3個參考站的電道組合及混合測道多參考處理結果，分別對應式(4)、式(5)兩種參考數據組合情況。可以看出，遠參考電道數據的參考效果不如磁道數據，混合測道的處理結果優于電道組合，但不如磁道組合。這與常規遠參考處理的測道選擇認識一致[6,14,34]。實際中，可通過噪聲分析，選擇含噪較低的測道組合進行處理，并經過對比評價，優選最終的處理方案。特別需注意的是，電道參考可能會較大程度地改變曲線形態，實際應用時需進一步論證其可靠性。

4 結論

(1)多參考站陣列數據處理技術有別于常規遠參考、互參考及多點遠參考疊加處理方法，該方法將多參考站、多參考道數據整合于統一的參考陣列數據集，利用陣列處理方法提取天然場源極化參數，可有效提高多參考站條件下的數據利用率。

(2)實測含“死頻帶”畸變的MT數據處理結果表明，多參考站陣列數據處理技術可獲得優于常規單站處理、常規陣列處理以及常規單參考站處理的結果。當存在多參考站條件時，可采用多種測站及測道組合進行處理MT數據，并對比優選處理結果。

需指出的是，強干擾條件下，“死頻帶”畸變MT數據的處理仍是十分棘手的問題。本文提供了一種改善質量的可行策略，但提高遠參考站本身的質量并盡可能提高采集數據的信噪比仍是改善目標測站數據處理質量的最關鍵因素。