大地電磁法的物理原理
大地電磁法利用太陽風與磁層交互作用以及赤道雷暴活動產生的天然電磁場作為信號源，通過同步測量地表電場 E 和磁場 H 的時變信號來推斷地下介質的電導率分佈。其核心物理在於電磁波的趨膚效應：不同頻率的電磁波穿透至不同深度。對於均勻半空間，趨膚深度 δ ≈ 503·√(ρ/f) 米，其中 ρ 為電阻率，f 為頻率。這意味著通過測量從 10⁴ Hz 到 10⁻⁴ Hz 的超寬頻帶響應，可以獲得從近地表數百米到上地函數百公里的電導率資訊。
MT 方法的基礎數學關係由阻抗張量 Z 描述：在頻域中，水平電場分量與磁場分量之間存在線性關係 E = Z·H。阻抗張量的四個複數元素 Z_xx、Z_xy、Z_yx、Z_yy 包含了地下電性結構的全部資訊。在二維構造中，對角元素為零，而非對角元素直接對應 TE 和 TM 模式的視電阻率和相位。
阻抗張量分析與維度判定
在進行反演之前，必須對 MT 數據進行嚴格的維度分析，以確定地下結構適合用一維、二維還是三維模型來描述。最常用的工具是相位張量分析：相位張量 Φ = X⁻¹·Y，其中 X 和 Y 分別為阻抗張量的實部和虛部。相位張量的橢圓率 β 表徵了地下結構的三維性程度：|β| < 3° 通常認為可近似為二維結構。
電性走向分析則用於確定二維構造的走向方向。基於阻抗張量分解，走向角可通過旋轉阻抗張量使得對角元素最小化來確定。在實際數據處理中，由於存在固有的 90° 模糊性，通常需要結合地質資訊或利用垂直磁場傳遞函數的感應向量來消除走向歧義。
三維反演與地函導電率成像
三維 MT 反演的目標是找到一個地下電導率模型 m，使得正演響應 F(m) 與觀測數據 d 之間的差異最小化。這是一個高度非線性且不適定的逆問題，通常通過正則化來求解。目標函數包含數據擬合項和模型約束項：Φ(m) = ||W_d·(F(m) - d)||² + λ·||W_m·(m - m_prior)||²。其中 W_d 為數據協方差加權矩陣，W_m 為模型平滑約束算子，λ 為正則化參數。
現代三維 MT 反演使用有限差分或有限元素法進行正演計算，結合非線性共軛梯度法或高斯-牛頓法進行迭代優化。對於包含數百個 MT 測點的區域性陣列，三維反演通常需要在高性能計算集群上運行數天。反演結果揭示了地殼和上地函中與部分熔融、揮發份含量及地函交代作用相關的導電異常體，為理解板塊構造和深部地質過程提供了關鍵約束。
MT 數據處理與阻抗估計
以下代碼實現了穩健的 MT 阻抗張量估計，使用迭代加權最小二乘法抑制非高斯噪聲。
import numpy as np from scipy.linalg import lstsq class RobustImpedanceEstimator: # Robust magnetotelluric impedance tensor estimation def __init__(self, max_iter=10, tolerance=1e-4): self.max_iter = max_iter self.tolerance = tolerance self.Z = None # estimated impedance tensor def estimate(self, Ex, Ey, Hx, Hy, Hx_remote, Hy_remote): # Robust remote-reference impedance estimation N = len(Ex) Ex, Ey = Ex.reshape(-1, 1), Ey.reshape(-1, 1) Hx_r = np.column_stack([Hx_remote, Hy_remote]) weights = np.ones(N) for iteration in range(self.max_iter): W = np.diag(np.sqrt(weights)) Hx_r_w = W @ Hx_r Ex_w = W @ Ex; Ey_w = W @ Ey Z_xx, Z_xy = lstsq(Hx_r_w, Ex_w)[0].flatten() Z_yx, Z_yy = lstsq(Hx_r_w, Ey_w)[0].flatten() self.Z = np.array([[Z_xx, Z_xy], [Z_yx, Z_yy]]) # Compute residuals and update Huber weights Ex_pred = Z_xx * Hx_remote + Z_xy * Hy_remote Ey_pred = Z_yx * Hx_remote + Z_yy * Hy_remote r = np.sqrt((Ex.flatten() - Ex_pred)**2 + (Ey.flatten() - Ey_pred)**2) sigma = np.median(r) * 1.4826 weights = np.where(r < 1.5 * sigma, 1, 1.5 * sigma / r) if sigma < self.tolerance: break return self.Z def apparent_resistivity(self, freq, mu0=4e-7*np.pi): # Calculate apparent resistivity from impedance omega = 2 * np.pi * freq rho_xy = np.abs(self.Z[0, 1])**2 / (omega * mu0) rho_yx = np.abs(self.Z[1, 0])**2 / (omega * mu0) return rho_xy, rho_yx # Example: synthetic 1000-sample dataset at 1 Hz np.random.seed(42) f = 1.0; N = 1000 Hx_r, Hy_r = np.random.randn(N), np.random.randn(N) Z_true = np.array([[0.01+0.02j, 0.15+0.08j], [-0.14-0.07j, -0.01-0.01j]]) Ex_s = Z_true[0,0]*Hx_r + Z_true[0,1]*Hy_r + 0.05*np.random.randn(N) Ey_s = Z_true[1,0]*Hx_r + Z_true[1,1]*Hy_r + 0.05*np.random.randn(N) estimator = RobustImpedanceEstimator() Z_est = estimator.estimate(Ex_s, Ey_s, np.random.randn(N), np.random.randn(N), Hx_r, Hy_r) rho_xy, rho_yx = estimator.apparent_resistivity(f) print(f"Apparent resistivity: ρ_xy={rho_xy:.1f} Ω·m, ρ_yx={rho_yx:.1f} Ω·m")
從岩石圈到地函轉換帶的電性約束
大地電磁法最獨特的貢獻在於其對地球深部揮發份和部分熔融的高靈敏度。實驗岩石物理表明，地函橄欖岩的電導率對微量氫（以 OH⁻ 形式存在於名義無水礦物中）極為敏感：僅數十 ppm 的水含量即可使電導率提升一個數量級。因此，MT 成像揭示的地函高導異常通常對應於俯衝板塊脫水、地函楔熔融或地函熱柱上升等關鍵地質過程。結合地震層析成像和重力數據，MT 為理解板塊俯衝帶的流體循環、大陸岩石圈的形成與破壞、以及深部碳循環提供了不可或缺的電性約束。
本文內容僅供學術研究參考。大地電磁數據的解釋存在固有的非唯一性，應結合多種地球物理和地質資料進行綜合約束。