一種時空陣列電磁數據處理模型

周聰，湯井田，原源，蘭學毅，郭冬

（1.東華理工大學核資源與環境國家重點實驗室，江西南昌，330013；2.中南大學有色金屬成礦預測與地質環境監測教育部重點實驗室，長沙，410083；3.安徽省勘查技術院，安徽合肥，230041）

電磁勘探系統包括輸入、輸出以及地球介質等，輸入端的場源條件、輸出端的觀測方式以及輸入-輸出間的相互關系等，是電磁勘探的基本研究內容，也是電磁勘探方法分類的重要依據。輸入端有多種場源形式，主要包括天然電磁場、不可控人文場以及可控人工場等，它們各自具有不同的特征。天然電磁場的場源主要包括太陽風、雷暴等，它們通常距離遙遠，場可近似視為以平面波形式垂直入射地球表面，具有很寬的頻譜，并且在各個方向上的強度基本相當[1-2]，但單頻點的信號弱，信噪比低。由于人類電磁活動加劇，人文電磁場的影響日益嚴重[3-6]，其以非平面波為主要成分，常導致基于平面波理論的常規數據處理結果出現畸變。特別地，在高頻“死頻帶”（1~5 kHz）及低頻“死頻帶”（1 Hz 附近）信號強度極低，阻抗更易受干擾影響[7]。為提高觀測數據信噪比，人工場源被引入進行信號發送，以增強場源信號強度[8-9]。由于人工場源的引入，場源的影響隨之出現，如產生非平面波效應、場源附加效應、場源陰影效應等[10-11]。為克服這些問題，研究者們提出了許多改進方案，主要包括增大發送功率以提高信噪比[12-13]、使用超大極距的固定發射臺以擴大波區的測量范圍[14-15]、采用編碼電磁信號及相關檢測技術提高數據的信噪比[16-18]以及采用張量發送和接收以擴大測量范圍提高信噪比[19-20]等。總體而言，這些改進方案可有效提高數據觀測質量，但難以同時利用天然場、人文場信息，也將導致采集成本增加。在輸出端，測量儀器與觀測手段的進步不斷推動著電磁勘探技術的發展，主要觀測方式逐步從一維、二維向三維過渡，對數據處理方法也不斷提出新的要求。近年來，分布式測量及三維陣列勘探受到了人們廣泛關注[21-24]，有效提高了觀測精度與空間分辨率，成為電磁觀測的新趨勢。然而，數據處理仍然只針對單個測站或部分場分量；天然場源電磁法數據處理一般針對單測站的張量電磁場進行，通過穩健估計提高處理質量，可控源電磁法數據處理常針對單測站的正交標量電磁場或單個電磁場分量[9-10]。多站數據的利用可以有效提高數據處理質量。遠參考技術[25]最早實現了不同測站間數據的利用。通過引入參考道的觀測數據，與測站的數據進行相關處理，以壓制噪聲，取得了很好的處理效果[2，26-30]。EGBERT等[31-33]提出了針對大地電磁法的頻率域陣列電磁數據處理技術，給出了基本數學模型與處理流程；GIUSEPPE 等[34]給出了該方法進行信噪分離的實例。但常規陣列電磁數據處理方法以頻域數據刪選的方式壓制相關噪聲的影響[31]，效果仍受制于高質量數據在整個觀測時段中的長度比例，難以處理持續型的含噪數據，且無法處理含人工源信號的混場源數據。為解決常規電磁數據處理方法輸入端相關噪聲難以壓制的瓶頸問題以及輸出端單站處理難以滿足陣列觀測數據處理需求的問題，本文作者提出一種陣列電磁數據處理新模型。假設多種電磁場源同時入射地球表面，多個測站同步觀測電磁場各分量，基于電磁勘探系統的線性時不變性，建立多輸入-多輸出的系統分析方法，以有效識別并分離天然電磁場、不可控人文電磁場以及可控人工電磁場，為頻域陣列電磁數據處理提供理論依據。

1 時空陣列電磁勘探模型

假設電磁勘探系統滿足線性時不變性[35]，并且其輸入端同時存在多個不同類型的場源（包括天然電磁場源、人工電磁場源和人文電磁場源），可以隨時間變化。即系統輸入端在空間上存在多個場源，每個場源對系統都有多個時刻的獨立激勵，各場源的激勵組合線性無關，而非簡單重復，構成空間陣列場源在不同時刻對系統的輸入時空陣列。輸出端存在多個測站，每個測站可包含多種類型的測道，記錄電場、磁場、電流密度以及電磁場梯度等分量。各測站同步觀測系統對輸入端所有場源在不同時刻激勵的響應，構成空間陣列測道對系統的輸出時空陣列。圖1所示為該系統的典型空間部署示意圖。

由于系統輸入、輸出兩端在時間、空間這2個維度均呈陣列分布，可稱該系統為時空陣列電磁勘探系統。圖2所示為該系統的簡化示意圖，包含多源輸入、大地及多道輸出3部分，對地電響應的研究可歸結為對多場源-多時刻激勵-多道接收的系統研究。

圖1 陣列電磁勘探系統的典型空間部署示意圖Fig.1 Typical spatial deployment of array electromagnetic exploration system

圖2 一種陣列輸入-陣列輸出的電磁勘探系統示意圖Fig.2 A schematic of an array input-output electromagnetic exploration system

2 時空陣列方程組

基于圖2所示的陣列電磁勘探模型，進行如下假設：

1）主要研究準靜態極限條件下的頻率域電磁場（如電場、磁場、磁場梯度、電流密度等），不考慮位移電流的影響。

2）系統中包含L個場源，共激勵有限長的觀測時間T，T 可分為I 個等長的激勵時窗；假設第l（l = 1，2，…，L）個場源在第i（i = 1，2，…，I）個時窗內的時域、頻域輸入極化電流強度分別為sli（t）和Sli（ω）。

3）系統中包含K個測道（包括不同分量的電場、磁場及其梯度分量等）和I個觀測時窗；對每個時窗進行傅里葉變換，得到相應的頻譜數據。第k（k =1，2，…，K）測道在第i（i = 1，2，…，I）個時窗的時域、頻域觀測值分別記為xki（t）和Xki（ω）。

4）系統滿足線性時不變條件，即大地系統的時域單位沖激響應和頻率特性響應不隨時窗的變化而變化，并且不同成分的輸出信號、輸出端噪聲滿足線性疊加關系。

對于l場源在i時窗內的輸入信號sli（t），記第k個測道對應的單位沖激響應為ψkl（t），輸出響應為則數字離散線性時不變系統輸入與輸出間的關系為

其中：*表示卷積；t 和τ 為采樣時刻；N 為i 時窗的采樣長度。注意ψkl（t）略去了時窗標號i。

由于k 測道處的觀測信號xki（t）為所有場源響應的疊加，考慮輸出端噪聲時，xki（t）可表述為

其中：rki（t）為輸出端噪聲。

式（2）即為圖2所示系統在線性時不變條件下的輸入-輸出關系。根據卷積定理，頻域觀測數據Xki（ω）滿足

其中：Ψkl（ω）和Rki（ω）分別為頻域系統響應、輸出端噪聲。

由于頻率域的各個頻點相互獨立，并進行獨立處理，為了表達方便，下文中均省去頻點符號，如將Sli（ω）記為Sli。對上述關系進行整體考慮，將方程寫成矩陣形式，有

各矩陣的上角標為其維數。

根據場源類型不同，可對上述關系進行拓展。假設天然場源的個數為L1，未知人文場源的個數為L2，可控人工場源的個數為L3，場源總數L =L1+ L2+ L3。分別以A，B和C標識天然場源、人文場源及可控人工場源的輸入矩陣，U，V和W分別為與天然場源、人文場源及可控人工場源相對應的頻域系統響應矩陣，則根據系統的線性疊加關系，時空陣列方程組可寫為

其中：

各矩陣的上角標為維數，上角標“*”表示矩陣轉置。

式（4）和（5）描述了線性時不變電磁勘探系統的輸入-輸出關系，稱為時空陣列方程組。不難理解，Ψ，U，V及W 反映了系統本身的特性以及各激勵源的空間信息，與激勵源隨時窗的變化無關，定義為空間模數矩陣，其中，元素相應稱為空間模數（spatial modes）。為方便論述，將時空陣列輸入矩陣S，A，B 及C 稱為極化參數矩陣，并將其中的元素相應稱為極化參數（polarization parameters），將時空陣列輸出矩陣X稱為時空陣列數據矩陣。在線性時不變電磁勘探系統的基礎上，對輸入端進行時空陣列激勵，對輸出端進行時空陣列觀測，通過時空陣列方程組描述輸入-輸出關系；通過直接求解時空陣列方程組，獲得場源的極化參數或各測道的空間模數，繼而提取場源信息或地電信息，即構成了時空陣列電磁數據處理模型。

需指出的是，與EGBERT[31]所提出的陣列電磁數據處理模型相比，本文模型由系統的線性時不變性推導出時空陣列方程組，并且考慮了存在可控人工源的輸入條件，更具一般性。

3 時空陣列方程組求解方案

基于時空陣列電磁數據處理模型，實現時空陣列方程組的求解，可分析場源信息或地電信息。對于電磁勘探，一般更關心地電結構的獲取，即從有限規模的觀測時空陣列數據矩陣X中提取出地電結構的目標參數（或簡稱目標參數），如卡尼亞視電阻率，其過程可描述如下：

其中：M為頻域的目標參數集合。

考慮采空區充填施工需要排壓孔，2個采井先用作排壓孔，待采空區充填穩定后，采用級配較好的中粗砂對2個采井回填治理，用振動棒將回填深度內的砂料振搗密實。考慮2個采井較深，先期充填的不規范，治理后還有繼續塌陷的可能，井口暫不封堵，采用鋼筋混凝土蓋板保護，根據塌陷情況，隨時回填處理。待塌陷完全穩定后，再用毛石混凝土封堵井口(圖4)。

為求取M，可將上述過程分為3 步：第一步，通過觀測合理的時空陣列，構建時空陣列數據矩陣X；第二步，利用數學方法，從X 中估計出Ψ，并保障結果的穩健性；第三步，從Ψ 中提取可供解釋的目標參數M，需剔除位置信息的影響（如采用各種視電阻率的定義方案以及目標參數的提取方案等）。

對上述第一步，可以利用所有的時空同步觀測值構建時空陣列數據矩陣，也可以將數據加以分類以構建不同種類的時空陣列數據矩陣[36]，進而得到不同類型的待求解方程組。一般地，根據數據集中主成分的不同，時空陣列方程組可以分為三大類，其解法也有所不同。

3.1 場源極化參數已知

在某些條件下，當觀測場中可控人工場占主導，而天然場和人文場能夠忽略時，或者利用其他預處理手段能夠獲得場源極化參數的估計值時，場源極化參數可能為已知信息。

當極化參數矩陣已知時，時空陣列方程組的求解較簡單，并可以求得空間模數矩陣的唯一值。假設通過適當的觀測方案以及前期處理已獲得極化參數矩陣S 的估計，考察時空陣列方程組（4）式可知，在數學上，該式為線性回歸模型。

最小二乘法是線性回歸模型的經典求解方法，其一般表達式為

其中：上角標“?”表示共軛轉置。

對于不相關噪聲R，可通過多時窗的疊加進行壓制。當觀測疊加次數足夠多，不相關噪聲的平方和疊加極小，即

其中：min表示極小。

3.2 場源極化參數未知

假設時空陣列數據矩陣中天然場占主導，L ≈L1，則時空陣列方程組為

式（9）中，觀測數據矩陣X 的維數為K × I，而待求解的空間模數矩陣U的維數為K × L，當I > L時，上述問題在數學上可歸為數據降維（dimension reduction）問題。基于特征值求解的方法是數據降維及矩陣分解的主要手段，奇異值分解是最簡單和最基本的方法。為了分析方便，以奇異值分解為例，對式（9）的求解進行說明。對矩陣X 進行奇異值分解：

在上述求解過程中，場源數量L 的確定是關鍵。由于天然場源及不可控人文場源的隨機性和不可預測性，很難準確知道其場源個數。一般地，在較寬的觀測頻段范圍內（如0.1~10 000 Hz），天然場源信號可被認為是均勻平面電磁波，其發收距因遠大于觀測深度及測點距離而可視作相等，不同方向場源的主要能量可以通過分解、合并等效為2個正交水平場源的作用，或者說天然場源的主要極化可等效視作2 個正交水平方向上的極化作用[31]。對于不可控的人文干擾場源，盡管難以確定場源的數量，但實際上無需求解準確的對應于每個人文場源的系統響應，只需求得人文場源的總體信號強度并予以分離即可保障天然場和可控人工場響應估計質量提高。因此，采用主成分分析法評估其等效能量繼而估計其等效場源數目及等效場源極化參數是合理的。

需指出的是，由于A未知，時空陣列方程組的分解并不唯一，因為

但張量空間模數的轉換函數具備唯一性，取

可見，當場源極化參數未知時，應求取張量空間模數的轉換函數T，并從T中提取目標參數。

3.3 場源極化參數部分已知

當輸入端包含可控人工場源，而觀測數據中的天然場或/和人文場也不可忽略時，極化參數部分已知。以包含天然場和可控人工場的時空陣列觀測數據集為例，對時空陣列方程組進行討論。

假設時空陣列數據矩陣中天然場和可控人工場占主導，人文場的影響可忽略，L ≈L1+ L3，則時空陣列方程組為

其中，X 和C 已知，而U，A，W 及R 均未知，并且滿足K × I ≥L（K + I）。求解式（14）可采用2 種方法。

一種方式是數學求解。式（14）可寫為非線性最小二乘的形式：

運用非線性最小二乘的分析方法或多變量含誤差的統計學方法可對式（15）進行分析求解。另外，由于高頻端常處于人工源的遠區，其響應一般受非平面波效應影響較小，可獲得地電參數的高質量估計值，繼而得到稍低頻點的W 初始估計值，對式（14）采用迭代計算可有望獲得U和W的更高質量估計。

另一種方法是結合物理方法進行求解[36]。由于場源極化參數反映的是場源隨時窗的變化，而天然場源的距離遠大于觀測深度及測點距離空間尺度，可假設一定區域內各測站處的天然場強度隨時窗具有相同的變化，即天然場源極化參數相當。在相對安靜的地區布設同步遠參考測站，保證遠參考測站中以天然場信號為主，進而選擇合適的數據集，構建主要包含天然場信號的參考觀測數據矩陣Xr，可得

其中：Ur和Rr分別為相應的空間模數和輸出端的不相關噪聲。注意Xr的構建可僅考慮單參考站的2道水平磁場測道，也可視情況加入水平電場和垂直磁場測道，當存在多遠參考站時，還可利用多站多道參考數據。

4 目標參數提取方法

針對不同的研究對象，需要獲取的解釋信息并不一致。如通訊電磁領域常關心場源位置及變化信息，對應地，需求解極化參數矩陣。電磁勘探領域關心的是地球系統本身的結構信息，因此，需求解空間模數矩陣，繼而提取目標參數。需指出的是，此處所指的目標參數是指視電阻率、傾子矢量、磁場梯度響應函數等地球物理反演及分析所用的參數，并非真實的地電結構物性參數。目標參數的提取方案主要分為以下兩大類。

1）當場源極化參數S及Ψ中包含的空間位置信息（即測點與場源間的距離及相對方位）已知時，可推導出Ψ 與地電參數（σ，ε，μ，η）間的顯性表達式，繼而通過對Ψ 求反函數的方式獲得目標參數的定義：

2）當S及Ψ中包含的空間位置信息未知時，無法從Ψ 與地電參數（σ，ε，μ，η）間的關系中獲得目標參數的定義，此時，可先求取轉換函數T：

其中：Ψ1和Ψ2分別為2 種類型的空間模數張量。根據式（13），T 具有唯一性，不再包含時變場源極化信息的影響。

通過合理定義，可從T中獲取目標參數[36]：

其中：函數f表示M與T間的映射關系，如以同一測站的電場和磁場空間模數張量求取阻抗張量、以不同測站的電磁場空間模數張量求取站間轉換函數等。

需說明的是，若系統中包含的可控人工場源極化方向單一，則由于方程組求解得到的人工源空間模數僅對應于單一方向，此時可提取人工源的標量參數，如場分量視電阻率、標量阻抗視電阻率等。

5 案例分析

基于本文提出的模型開展了數值模擬試驗[36]，論證了模型的可行性。2018年7月，在青海柴達木盆地大柴旦地區開展了實際觀測與處理試驗。試驗區地處戈壁灘，人文背景噪聲較低，是較理想的電磁觀測與處理試驗場地。為模擬多場源環境，在試驗區中心點約2 km 外采用加拿大鳳凰公司的V8 可控源發送機間斷發送不同頻率的方波信號。接收端同時投入5臺設備，觀測4個同步陣列測站與1個遠參考站，遠參考站距試驗區約100 km。每個測站均以音頻大地電磁法觀測方式進行水平張量觀測，響應頻率范圍為1~104Hz，觀測時間達200 min。

本例中，假設外加人工源信號場源未知，以模擬未知人文場源。依據本文所述的時空陣列電磁數據處理模型，可建立時空陣列方程組。基于前述方程組求解方案，采用差分求解策略[37]，以4個測站的水平磁場差分數據提取未知人文場源極化參數，以遠參考站的多測道陣列數據提取天然場源極化參數，繼而求解空間模數及目標參數（阻抗視電阻率、相位）。為進行對比，還對觀測數據進行常規數據處理。

觀測陣列中某測點處理結果如圖3所示，其中圖3（a）和圖3（c）所示為常規穩健估計處理結果，圖3（b）和圖3（d）所示為時空陣列電磁數據處理結果。從圖3 可以看出：試驗點常規處理結果視電阻率-相位曲線（ρxy，φxy）在1~30 Hz處產生了畸變，為典型的“近源型”畸變曲線[3]，顯然是受到了外加人工場源的影響；視電阻率-相位曲線（ρyx，φyx）沒有明顯“近源型”畸變，其原因是人工源極化方向單一，對該模式方向的電磁場影響較弱，但ρyx和φyx曲線仍有脫節及“飛點”現象，數據顯然存在偏差。而采用時空陣列電磁數據處理方法所得到的結果，不再表現出“近源型”畸變形態；ρxy曲線低頻段的虛假高阻異常得到了壓制，反映出該地深部地層的低阻特征，5~100 Hz頻段內相位曲線趨于0的畸變數據得到了校正，相位數據恢復正常；ρyx曲線在10 Hz 附近頻段的脫節及“飛點”得以校正，曲線整體更光滑，相位曲線低頻段趨于0的畸變數據也得到了校正，5 Hz以上的相位均恢復了正常，結果更合理。對比常規處理結果，本文方法對含相關噪聲明顯的ρxy，φxy低頻畸變數據進行了恢復，對含噪相對不明顯的ρyx，φyx數據質量進行了改善，取得了更可信的效果。需指出的是，在低頻端數個頻點處，即使按本文方法處理，結果也仍存在“飛點”，其原因是：低頻端處于大地電磁的“死頻帶”范圍內，信號強度極低；低頻端的功率譜數據樣本更少，陣列方程組的求解精度相對較低。這也說明在成本允許的前提下，在低信噪比地區，應盡可能延長采集時間，獲得更多的低頻功率譜，進而提高低頻響應估計質量。該實例說明本文所提出的時空陣列電磁數據處理模型具有較好的實用價值。

6 結論

1）提出了一種時空陣列電磁數據處理模型，適用于同時存在天然場源、可控人工源和不可控人文場源的陣列輸入條件以及多站、多道、多觀測時窗的陣列輸出條件。基于電磁勘探系統的線性時不變性推導出了時空陣列方程組，論證了時空陣列方程組的可解性與求解條件，闡述了目標參數的提取方法。本文模型為實現輸入端的信噪分離、輸出端的多站多道數據綜合分析提供了理論依據。

2）時空陣列電磁數據處理模型基于地球電磁系統的線性時不變性，具有廣泛的適用性。不同于常規數據處理方法，本文模型將輸入端所有能量不可忽略的信號均視為有效信號，繼而尋求相應的信噪分離方案，從理論上為相關噪聲的壓制提供了新思路；通過陣列數據的分類分析，可充分挖掘數據中各類信息，提高數據利用率。

3）本文模型旨在為多輸入-多輸出的電磁勘探提供一種新的數據分析方案，但并非在所有條件下均適用；模型的可行性和適用性有待進一步論證。對于天然場、人工場和人文場均不可忽略時的場源環境，其時空陣列方程組的高質量求解策略仍需進一步研究。