王琪琪，馬 濤，王紅兵，譚 鑫，高 鑫，汪強強

(中鐵第一勘察設計院集團有限公司 陜西鐵道工程勘察有限公司，陜西 西安 710000)

1 引 言

隨著大地電磁數(shù)據(jù)處理解釋技術的迅速發(fā)展，MT(Magnetotelluric)的應用領域已經(jīng)從簡單的地質構造向復雜的地質構造發(fā)展[1-3]。但在這些區(qū)域構造離地表很近的地方往往存在著三維局部電性不均勻體，當進行大地電磁法測深，局部不均勻體的邊界上會形成積累電荷，使周圍的電場突然變大或者減小，從而引起靜態(tài)效應[4-6]。靜態(tài)效應往往會使視電阻率擬斷面圖出現(xiàn)陡立密集的等值線異常，很容易將這種異常解釋為地下存在直立的大斷裂或縱向延伸的巖脈，同時靜態(tài)效應也會嚴重影響反演后的定量解釋，如果直接用發(fā)生靜位移后的視電阻率曲線進行反演，會得到錯誤的巖層電阻率值和界面深度信息[7]。所以正確認識分析靜態(tài)效應的響應特征曲線，選擇合適的方法識別并進行校正，對大地電磁測深數(shù)據(jù)的處理與反演解釋至關重要[8,9]。

近年來，大地電磁法阻抗張量分解研究工作不斷深入，得到了快速發(fā)展。趙國澤等[10]將Bahr分解應用到實測電磁數(shù)據(jù)的處理中，對區(qū)域異常分布做了一定的說明；王書明[11]首次將廣義逆方法引入到GB分解，并用局部畸變模型進行驗證。尹兵祥等[12]提出一種在基于三維構造的正則分解方法，對其參數(shù)的物理意義和地質含義進行了解釋；楊長福等[13]用視電阻率和相位代替了阻抗張量去實現(xiàn)GB分解，使方程和未知參數(shù)的個數(shù)減少，提高了計算效率。晉光文等[14]利用Swift和Bahr分解方法對川西-藏東地區(qū)的大地電磁阻抗數(shù)據(jù)進行了校正，分離了區(qū)域阻抗和畸變因子，并分析了該地區(qū)的構造分布；蔡軍濤等[15]通過復數(shù)定義提出一種“共軛阻抗變換”分解方法，分解后的阻抗具有和相位張量相似的性質，并得到了區(qū)域主軸方位角和畸變因子，同時根據(jù)轉換后的阻抗定義了新的維性判別參數(shù)。謝成良等[16]用相位張量得到電性主軸和維性信息，建立合理的GB分解初始模型，同時將各向異性參數(shù)作為未知參數(shù)進行計算。

隨著大地電磁法勘探規(guī)模和計算機處理能力的增大，多測點同步觀測將會成為一種趨勢，而已有的多頻多測點阻抗分解程序需要計算的參數(shù)會隨測點和頻點數(shù)呈線性增長。針對傳統(tǒng)局部優(yōu)化算法，如果未知參數(shù)沒有初值約束，視電阻率的計算結果可能會嚴重偏離真實值，陷入局部最優(yōu)，而利用多個測點間的聯(lián)系來消除靜態(tài)效應成為了可能[17,18]。

本文利用靜態(tài)因子在區(qū)域分布的統(tǒng)計學性質實現(xiàn)了1D靜態(tài)校正,并編寫了基于粒子群算法和高斯牛頓法的畸變校正程序，通過正演模型驗證了該算法的有效性，發(fā)現(xiàn)該方法適用于存在多個淺部異常體的區(qū)域1D電磁數(shù)據(jù)。

2 基本原理

2.1 G B分解

1989年，Groom和Bailey[19]對阻抗張量進行了如下分解：

Zm=RCZRT

(1)

其中，R是旋轉矩陣;Z是區(qū)域未受畸變的一維阻抗張量;C是一個與頻率無關的實數(shù)畸變矩陣，C=gTSA；g是增益因子，為無量綱參數(shù)；T是扭曲矩陣；S是剪切矩陣；A是各向異性因子矩陣，它們的矩陣分別表示如下：

式(2)～式(6)中，θ是主軸方向角，單位為度；Z1為一維區(qū)域阻抗張量元素，單位為Ω；t是扭曲因子；e是剪切因子；s是各向異性因子；t,e,s均為無量綱參數(shù)。其中t、e組成了畸變矩陣C的可確定部分，g、s組成了不可確定部分，所以實測阻抗又可寫為下式：

(7)

2.2 畸變校正

李洋[20]提到對于局部優(yōu)化算法，初值的確定有很大的影響，本文先通過粒子群全局優(yōu)化算法得到一個包含畸變參數(shù)和區(qū)域阻抗的全局最優(yōu)解空間，然后將全局最優(yōu)解空間作為高斯牛頓法的初值進行計算。該方法能夠更好地避免陷入局部最優(yōu)，得到更真實的畸變參數(shù)，使得校正后區(qū)域阻抗更可靠。算法采用的目標函數(shù)如下所示：

利用式(7)得到如下參數(shù):

其中,σ=A+B,δ=A-B。建立如下目標函數(shù):

(12)

2.2.1 粒子群算法

粒子群算法是一種全局隨機搜索算法[21]。該算法模擬鳥群覓食行為，將鳥看做粒子，通過粒子個體之間共享信息來找尋最優(yōu)解。每個粒子有速度和位置兩個屬性，粒子通過自身當前在搜索路徑中的最優(yōu)值和整個粒子群中的最優(yōu)值來調整速度和位置，其更新方式如下所示：

其中Vi、xi為第i個粒子的速度和位置，ω為慣性權重，c1、c2為學習因子，控制粒子的認知模式，ξ和η為0～1間的隨機數(shù)，pbest、gbest分別為粒子個體最優(yōu)值和全局最優(yōu)值，r為約束因子。當目標函數(shù)式滿足規(guī)定的閾值或者達到最大迭代次數(shù)時迭代結束。

2.2.2 高斯牛頓法

目標函數(shù)公式(12)二階泰勒展開如下所示：

(15)

其中,a′是每次迭代得到的畸變參數(shù)和區(qū)域阻抗;c為a′對應的目標函數(shù)值;H為海森矩陣。則γ2(a)梯度為：

(16)

(17)

因為在H(a(k))矩陣中二階項很小，則H(a(k))≈JTJ，代入上式得：

(18)

當步長足夠小或者迭代次數(shù)達到閾值時停止迭代，此時a(k+1)為最優(yōu)解。

2.3 增益因子

上述步驟對GB分解中確定參數(shù)進行了擬合，可以得到僅受靜態(tài)效應影響的阻抗數(shù)據(jù)，本節(jié)主要對不確定參數(shù)進行校正，具體步驟如下：

當某一測區(qū)的測點數(shù)足夠多時，地表不均勻體產(chǎn)生的電流畸變在統(tǒng)計計算上將會被光滑[22]，利用這一假設，得到多測點多頻點下的SSQ旋轉不變量：

(19)

其中,ri為第i個測點;ω為角頻率，單位為rad/s。區(qū)域構造為1D時，整個測區(qū)的SSQ旋轉不變量為：

(20)

對應的每個測點視增益表示為：

(18)

gi為第i個測點對應的視增益因子。

3 正演模型分析

建立一個如圖1所示的二層模型,第一層厚度為2 km，電阻率為100 Ω·m，第二層電阻率為1 000 Ω·m，淺部存在三個低阻異常體，電阻率為1 Ω·m，x和y方向延伸400 m，z方向延伸500 m。

圖1 正演模型

由圖2可以看出，當?shù)?00次之前，畸變參數(shù)的擬合誤差發(fā)生快速下降，之后隨著迭代次數(shù)的增多誤差不再發(fā)生變化，約為3.8×10-4，表明該算法收斂速度較快，擬合結果比較精確。

圖2 高斯牛頓法迭代過程中均方根誤差曲線

圖3和圖4顯示了頻率分別為350 Hz和0.07 Hz時不同測點位置上畸變參數(shù)和靜態(tài)因子的橫向分布，圖中冒號后面的數(shù)字表示不同畸變參數(shù)隨測點變化的標準差?？梢园l(fā)現(xiàn)整體上增益因子g和各向異性因子s隨測點的變化比較大，g存在兩個極大值和三個極小值，極大值對應異常體間的測點，極小值對應三個異常體正上方，s的分布與之相反。扭曲因子t和剪切因子e隨測點變化較小，標準差都在0.005以下，其中頻率為320 Hz時靜態(tài)因子橫向分布的標準差小于0.07 Hz，表明當頻率較高時,地表異常體對視電阻率的影響不再是靜態(tài)畸變，而是作為一個局部異常體。

圖3 頻率為320 Hz時畸變參數(shù)和靜態(tài)因子分布

圖4 頻率為0.07 Hz時畸變參數(shù)和靜態(tài)因子分布

圖5為27號測點校正前后的視電阻率和相位曲線，測點位于異常體一側，其中xy_correct、yx_correct分別表示校正后的電阻率和相位，xy_distort、yx_distort分別表示存在淺部異常體時的電阻率和相位，xy_nodistort、yx_nodistort分別表示不存在淺部異常體時的電阻率和相位。由圖5(a)和圖5(b)可以發(fā)現(xiàn),地表異常體的存在使得xy模式和yx模式視電阻率發(fā)生了明顯分離，xy模式下移，yx模式上移，相位曲線在高頻段出現(xiàn)了分離，受到明顯靜態(tài)效應影響，校正后的視電阻率與未含異常體的視電阻率幾乎重合，原來分離的視電阻率發(fā)生閉合；圖5(c)和圖5(d)中10 Hz處靜態(tài)因子存在一個拐點p，與圖5(a)和圖5(b)對照發(fā)現(xiàn)，頻率高于10 Hz時，視電阻率和相位曲線包含的是淺部異常體的信息，低于10 Hz時，區(qū)域構造受異常體影響曲線發(fā)生偏移，此時圖5(c)中拐點p所指向的位置，可以作為判斷靜態(tài)效應影響頻率范圍的一個依據(jù)。

圖5 27號測點和畸變參數(shù)隨頻率變化曲線

圖6為31號測點校正前后的視電阻率和相位曲線，測點位于異常體正上方。xy模式和yx模式受靜態(tài)效應的影響發(fā)生向上偏移，和27號測點一樣，校正后的xy模式和yx模式視電阻率曲線相互重合，表現(xiàn)出明顯的1D性質，說明該方法很好地消除了畸變的影響。

圖6 31號測點校正前后曲線和畸變參數(shù)隨頻率變化曲線

4 結 論

本文針對部分1D區(qū)域構造，引入SSQ旋轉不變量來計算得到增益因子，結合GB分解實現(xiàn)了區(qū)域阻抗的恢復。同時以正演模擬數(shù)據(jù)為例，獲得了能夠反映真實地下介質分布的視電阻率曲線和相位以及相關畸變參數(shù)，方便后續(xù)反演工作的進行，驗證了該方法的有效性。最后針對論文內容總結了以下幾點：

1)在多測點多頻點GB分解的基礎上，將粒子群算法和高斯牛頓法結合，避免計算結果陷入局部最優(yōu)。

2)通過建立的畸變模型發(fā)現(xiàn)，相比視電阻率相位受畸變影響更小，并且影響大小隨頻率發(fā)生變化，所以相位可以作為判斷區(qū)域維性的重要參數(shù)。同時針對均勻淺部異常體，不確定畸變參數(shù)(增益因子和各向異性因子)呈“正態(tài)分布”，在異常體中心最大，逐漸向兩側減小，同時發(fā)現(xiàn)不確定參數(shù)存在一個拐點，可以用來確定靜態(tài)效應的影響頻段。

3)通過正演模擬數(shù)據(jù)表明，本文提出的校正方法針對于1D區(qū)域構造能夠很好地壓制多個地表異常體的畸變影響，校正后的視電阻率和相位曲線幾乎重合，能夠更真實地反映地下介質分布情況。