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

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

1 引 言

隨著大地電磁數(shù)據(jù)處理解釋技術的迅速發(fā)展，MT(Magnetotelluric)的應用領域已經(jīng)從簡單的地質(zhì)構(gòu)造向復雜的地質(zhì)構(gòu)造發(fā)展[1-3]。但在這些區(qū)域構(gòu)造離地表很近的地方往往存在著三維局部電性不均勻體，當進行大地電磁法測深，局部不均勻體的邊界上會形成積累電荷，使周圍的電場突然變大或者減小，從而引起靜態(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]提出一種在基于三維構(gòu)造的正則分解方法，對其參數(shù)的物理意義和地質(zhì)含義進行了解釋；楊長福等[13]用視電阻率和相位代替了阻抗張量去實現(xiàn)GB分解，使方程和未知參數(shù)的個數(shù)減少，提高了計算效率。晉光文等[14]利用Swift和Bahr分解方法對川西-藏東地區(qū)的大地電磁阻抗數(shù)據(jù)進行了校正，分離了區(qū)域阻抗和畸變因子，并分析了該地區(qū)的構(gòu)造分布；蔡軍濤等[15]通過復數(shù)定義提出一種“共軛阻抗變換”分解方法，分解后的阻抗具有和相位張量相似的性質(zhì)，并得到了區(qū)域主軸方位角和畸變因子，同時根據(jù)轉(zhuǎn)換后的阻抗定義了新的維性判別參數(shù)。謝成良等[16]用相位張量得到電性主軸和維性信息，建立合理的GB分解初始模型，同時將各向異性參數(shù)作為未知參數(shù)進行計算。

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

本文利用靜態(tài)因子在區(qū)域分布的統(tǒng)計學性質(zhì)實現(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是旋轉(zhuǎn)矩陣;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)值來調(diào)整速度和位置，其更新方式如下所示：

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

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旋轉(zhuǎn)不變量：

(19)

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

(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，表明該算法收斂速度較快，擬合結(jié)果比較精確。

圖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號測點校正前后的視電阻率和相位曲線，測點位于異常體一側(cè)，其中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ū)域構(gòu)造受異常體影響曲線發(fā)生偏移，此時圖5(c)中拐點p所指向的位置，可以作為判斷靜態(tài)效應影響頻率范圍的一個依據(jù)。

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

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

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

4 結(jié) 論

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

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

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

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