基于多測(cè)點(diǎn)多頻點(diǎn)的大地電磁一維區(qū)域阻抗恢復(fù)

王琪琪，馬 濤，王紅兵，譚 鑫，高 鑫，汪強(qiáng)強(qiáng)

(中鐵第一勘察設(shè)計(jì)院集團(tuán)有限公司 陜西鐵道工程勘察有限公司，陜西 西安 710000)

1 引 言

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

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

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

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

2 基本原理

2.1 G B分解

1989年，Groom和Bailey[19]對(duì)阻抗張量進(jìn)行了如下分解：

Zm=RCZRT

(1)

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

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

(7)

2.2 畸變校正

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

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

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

(12)

2.2.1 粒子群算法

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

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

2.2.2 高斯牛頓法

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

(15)

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

(16)

(17)

因?yàn)樵贖(a(k))矩陣中二階項(xiàng)很小，則H(a(k))≈JTJ，代入上式得：

(18)

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

2.3 增益因子

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

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

(19)

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

(20)

對(duì)應(yīng)的每個(gè)測(cè)點(diǎn)視增益表示為：

(18)

gi為第i個(gè)測(cè)點(diǎn)對(duì)應(yīng)的視增益因子。

3 正演模型分析

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

圖1 正演模型

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

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

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

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

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

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

圖5 27號(hào)測(cè)點(diǎn)和畸變參數(shù)隨頻率變化曲線

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

圖6 31號(hào)測(cè)點(diǎn)校正前后曲線和畸變參數(shù)隨頻率變化曲線

4 結(jié) 論

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

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

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

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