張文波，張瑩，李建慧

(中國地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，湖北 武漢 430074)

0 引言

瞬變電磁法在地下水調(diào)查、金屬礦產(chǎn)勘探、煤礦防治水等領(lǐng)域發(fā)揮著重要作用。瞬變電磁法三維反演是當(dāng)前的研究熱點(diǎn)，但實(shí)用化的資料處理解釋手段仍以一維反演為主[1-2]。

電磁法一維反演的常用方法有最小構(gòu)造反演、Occam反演、橫向約束反演(laterally constrained inversion，LCI)、空間約束反演(spatially constrained inversion，SCI)，以及模擬退火法、深度學(xué)習(xí)等多種人工智能優(yōu)化算法。一維最小構(gòu)造反演是求取多層介質(zhì)模型的最光滑解，即在一定擬合差條件下，使得模型粗糙度最小，因此，目標(biāo)函數(shù)不僅包含數(shù)據(jù)擬合項，也包含模型正則化項[3-4]。Occam反演與最小構(gòu)造反演類似，不同之處在于正則化因子的選取：Occam反演通過一元函數(shù)優(yōu)化來選取正則化因子，而最小構(gòu)造反演中正則化因子通常固定不變或者通過一些非常簡單的算法選取[5]，因此，Occam反演更加穩(wěn)定，對正則化初始值依賴程度低，但其計算量偏大。LCI方法同時反演電阻率和層厚度，并在目標(biāo)函數(shù)中加入相鄰測點(diǎn)之間電阻率、層厚度和界面深度約束項[6]，因此單條測線所有測點(diǎn)數(shù)據(jù)需整體反演。SCI與LCI類似，不僅考慮同一條測線相鄰測點(diǎn)之間的約束，還需考慮其他測線相鄰測點(diǎn)之間的約束[7]，因此測區(qū)內(nèi)所有測點(diǎn)數(shù)據(jù)需整體反演。

目前，國內(nèi)外多個課題組均已開發(fā)了瞬變電磁法一維反演軟件。如：加拿大英屬哥倫比亞大學(xué)Oldenburg課題組在國際上率先實(shí)現(xiàn)了用于回線源瞬變電磁法數(shù)據(jù)處理的一維最小構(gòu)造反演[8]，開發(fā)了實(shí)用化的軟件包EM1DTM，并成功應(yīng)用于多個野外實(shí)例[9-10]；德國科隆大學(xué)Tezkan課題組采用阻尼最小二乘法和Occam算法開展了回線源瞬變電磁法、長偏移距瞬變電磁法、射頻大地電磁法等方法的單一反演或2種以上方法聯(lián)合反演[11-14]，并開發(fā)了EMUPLUS軟件包；丹麥奧胡斯大學(xué)Auken課題組基于阻尼最小二乘法開發(fā)了瞬變電磁法一維LCI程序，用于解決反演電阻率和層界面橫向連續(xù)性差等問題[15]，還開發(fā)了SCI程序，用于解決反演電阻率和層界面水平方向連續(xù)性差的問題[7]；Kirkegaard和Auken在并行計算、迭代法求解方程組、靈敏度矩陣近似等方面優(yōu)化了算法，進(jìn)而提升了反演效率[16]。上述反演算法和策略均集成于AarhusInv軟件之中。此外，IX1D、Maxwell等商業(yè)軟件也都包含瞬變電磁法一維反演模塊。

國內(nèi)多個課題組也開展了相關(guān)研究。吉林大學(xué)殷長春課題組先后開發(fā)了適用于航空瞬變電磁法資料處理解釋的一維LCI和SCI程序[17-18]，采用Occam算法對m序列發(fā)射波形多道瞬變電磁法數(shù)據(jù)開展了反演研究[19]，還采用深度學(xué)習(xí)方法實(shí)現(xiàn)了航空瞬變電磁法的一維反演快速成像[20]；北京大學(xué)黃清華課題組開展了地面回線源瞬變電磁法的一維最小構(gòu)造反演研究，該算法適用于水平方向發(fā)生形變的復(fù)雜發(fā)射回線[21]；中國礦業(yè)大學(xué)(北京)程久龍課題組開展了粒子群優(yōu)化——阻尼最小二乘混合算法研究，并應(yīng)用于礦井瞬變電磁法探測[22]；山東大學(xué)孫懷鳳等采用模擬退火法開展了回線源瞬變電磁法一維反演，并取得了良好效果[23]。

綜上所述，瞬變電磁法一維反演技術(shù)已趨于成熟，并在資料處理中發(fā)揮著重要作用。近年來，筆者團(tuán)隊初步開發(fā)了一套地面回線源瞬變電磁法一維反演系統(tǒng)。該一維反演系統(tǒng)基于Fortran語言編寫，并融合了以下2個開源程序：Occam反演程序[24]與電磁場從頻率域向時間域轉(zhuǎn)換的正弦和余弦變換程序[25]。為了提高計算效率，我們對漢克爾變換、矩陣與矢量乘積、一般矩陣乘積等環(huán)節(jié)開展了基于OpenMP的并行計算。

1 反演算法原理

上述地面回線源瞬變電磁法一維反演系統(tǒng)包括基于高斯牛頓法的最小構(gòu)造反演和Occam反演，基于阻尼最小二乘法的橫向約束反演和空間約束反演(圖1)。其中，一維正演基于自研算法開發(fā)[26-27]，靈敏度矩陣采用擾動法計算。

圖1 一維反演框架Fig.1 The framework of the 1D inversion system

1.1 最小構(gòu)造反演與Occam反演

基于Tikhonov正則化方式構(gòu)建反演目標(biāo)函數(shù)φ，有：

φ=‖Wd(dobs-dpre)‖2+λ‖Wm(m-m*)‖2。

(1)

式中:右邊第一項為實(shí)測數(shù)據(jù)與正演模型響應(yīng)的擬合差函數(shù)，第二項為模型正則化函數(shù)，‖·‖表示l2范數(shù)；dobs為實(shí)測數(shù)據(jù)向量；dpre為模型參數(shù)向量m相應(yīng)的正演模型響應(yīng)向量，由瞬變電磁法一維正演計算；Wd為數(shù)據(jù)加權(quán)矩陣，此處為一個對角矩陣，其元素為實(shí)測數(shù)據(jù)標(biāo)準(zhǔn)差的倒數(shù)；Wm為模型粗糙度算子；λ為正則化因子；m*為參考模型，先驗(yàn)信息包含其中。這兩種反演算法中，m為各層介質(zhì)電導(dǎo)率。

對于目標(biāo)函數(shù)式(1)，采用高斯牛頓法[28]最優(yōu)化時，第k+1次迭代的正規(guī)方程為：

(2)

式中：Jk為正演函數(shù)在模型參數(shù)mk處的一階導(dǎo)數(shù)，即雅克比矩陣；δmk+1為模型更新量。最小構(gòu)造反演中，通常采用“冷卻法”更新正則化因子[29]。算法流程見圖2。

圖2 一維最小構(gòu)造反演流程Fig.2 The flow chart for the 1D minimum-structure inversion

Occam反演中，上述目標(biāo)函數(shù)和正規(guī)方程可保持不變，其與最小構(gòu)造反演的主要區(qū)別在于如何更新正則化因子，具體策略見文獻(xiàn)[5]和[28]。

1.2 LCI和SCI

LCI和SCI中，構(gòu)建如下目標(biāo)函數(shù)φ：

φ=‖Wd(dobs-dpre)‖2+α‖Rp(m-m*)‖2+

β‖Rh(m-m*)‖2。

(3)

式中:右邊第一項為實(shí)測數(shù)據(jù)與正演模型響應(yīng)的擬合差函數(shù)；第二項為模型電導(dǎo)率和層厚約束項，Rp和α分別為其約束矩陣和權(quán)重；第三項為模型層界面深度約束項，Rh和β分別為其約束矩陣和權(quán)重。其中，Rp和Rh具體形式可參考文獻(xiàn)[6],其他參數(shù)與上文一致。

對于目標(biāo)函數(shù)式(3)，采用阻尼最小二乘法[28]最優(yōu)化時，第k+1次迭代的正規(guī)方程為：

(4)

式中：I為單位矩陣，λ為阻尼因子，仍采用“冷卻法”更新。圖3為LCI和SCI流程圖，合成正規(guī)方程后，其流程與最小構(gòu)造反演一致。

圖3 一維LCI和SCI流程Fig.3 The flow chart for the 1D LCI and SCI

LCI和SCI的區(qū)別在于參與約束的測點(diǎn)不同。如圖4所示,LCI中對于某一測點(diǎn)(紅色)，參與約束的為該測點(diǎn)所在測線中的2個相鄰點(diǎn)(綠色)，而SCI中對于某一測點(diǎn)(紅色)，參與約束的為該測點(diǎn)的所有鄰近點(diǎn)(綠色)。本研究中，以某一測點(diǎn)為圓心，按特定搜索半徑自動尋找各個方向與該測點(diǎn)距離最近的測點(diǎn)，并將這些測點(diǎn)作為約束測點(diǎn)。

圖4 LCI和SCI參與約束測點(diǎn)示意Fig.4 The diagram showing the observation points used as constrained points in the LCI and SCI

2 應(yīng)用實(shí)例

2.1 實(shí)例概述

內(nèi)蒙古自治區(qū)那仁寶力格煤田勘探區(qū)內(nèi)地勢平坦，新近系火山噴出巖比較發(fā)育，火山巖漿噴溢時的通道穿越煤層時不但對煤層產(chǎn)生破壞，而且對煤的變質(zhì)也會有不同程度的影響。為了保證今后煤田的順利開采，瞬變電磁法被用于確定巖漿上涌通道，以及通道周邊玄武巖與圍巖的界面。

從火山口和火山錐的存在和分布形態(tài)看，巖漿以中心式噴發(fā)為主，裂隙式噴發(fā)次之；巖性組合主要有橄欖拉斑玄武巖和玄武巖?？碧絽^(qū)內(nèi)，玄武巖體電阻率可達(dá)上千歐姆米，而砂巖和泥巖地層電阻率僅為幾十歐姆米，二者電性差異明顯[30]。這也是采用瞬變電磁法圈定玄武巖巖體在沉積巖中分布范圍的前提。

如圖5所示，野外工作采用了大回線源裝置，發(fā)射回線邊長為600m×300m，其中600m邊長沿測線方向布置。本次野外工作布設(shè)了11條測線，線距為100m，點(diǎn)距為40m，共計486個測點(diǎn)。160測線為控制測線，穿越了整個玄武巖露頭，其他測線均位于玄武巖露頭中心區(qū)域。對于160測線，鋪設(shè)1次回線源只采集6個測點(diǎn)數(shù)據(jù)；對于其他測線，鋪設(shè)1次回線源可采集2條相鄰測線的12個測點(diǎn)數(shù)據(jù)。采用加拿大Phoenix公司生產(chǎn)的V8多功能電法儀，觀測時間序列為分布于0.5～20ms之間的30個時刻。

圖5 那仁寶力格煤田測點(diǎn)分布示意Fig.5 The distribution diagram of measuring points in the Narenbaolige coalfield

2.2 單點(diǎn)反演

為了驗(yàn)證算法的正確性，首先將本研究開發(fā)的最小構(gòu)造反演和Occam反演程序與商業(yè)軟件IX1D對比。對于最小構(gòu)造反演和Occam反演，初始模型為36層層狀介質(zhì)，每層電阻率均為101.5Ω·m，目標(biāo)擬合差為3。IX1D軟件使用上述初始模型時，在多個測點(diǎn)無法擬合實(shí)測數(shù)據(jù)，因此使用了軟件自動估算的初始模型。由于IX1D軟件無法設(shè)置目標(biāo)擬合差，其在多數(shù)測點(diǎn)的最終均方根擬合差(RMS)大于3(圖6所示)；而本研究開發(fā)的最小構(gòu)造反演收斂性良好，所有測點(diǎn)均方根擬合差均在3以內(nèi)。

圖6 IX1D軟件和本研究最小構(gòu)造反演結(jié)果的RMS對比Fig.6 The comparison of RMS for IX1D software and the minimum-structure inversion method presented here

圖7為160線由單點(diǎn)反演地電模型拼接而成的電阻率二維斷面。圖中，高阻異常體(玄武巖巖體)形態(tài)對稱，呈“鍋底狀”分布，并且3種方法的反演電阻率分布范圍基本一致。其中，IX1D反演結(jié)果在相鄰測點(diǎn)之間跳躍性較強(qiáng)，玄武巖巖體與沉積巖界面連續(xù)性較差。這一現(xiàn)象與部分測點(diǎn)最終均方根擬合差偏大有關(guān)。此外，3種反演結(jié)果與鉆孔揭露的玄武巖厚度吻合良好，僅在鉆孔ZK3處差異明顯。該鉆孔終孔深度為716.6m，未揭露下伏沉積巖，結(jié)合該鉆孔位于高阻異常體中心，推斷其極有可能位于火山通道分布范圍之內(nèi)。因此，一維反演能夠獲得玄武巖巖體除巖漿上涌通道區(qū)域外的分布形態(tài)，卻無法直接反映巖漿上涌通道信息。這一現(xiàn)象與電磁場的傳播規(guī)律有關(guān)，即隨著探測深度的增大，電磁成像分辨率隨之下降，也與玄武巖巖體“倒錐形”的真實(shí)形態(tài)有關(guān)。

圖8為測區(qū)最小構(gòu)造反演結(jié)果俯視圖，色標(biāo)與圖7保持一致。圖中，測線中心區(qū)域電阻率最大，色標(biāo)為桔紅色至紅色，測線兩端電阻率偏低，顯示為綠色甚至藍(lán)色。顏色的這種分布規(guī)律一定程度上反映了玄武巖厚度。另外，盡管多數(shù)測點(diǎn)與其四周相鄰測點(diǎn)的電阻率在同一個色標(biāo)范圍內(nèi)，但它們之間電阻率差異仍然較大。

圖7 160線單點(diǎn)反演電阻率斷面Fig.7 The section view stitched from the single-point 1D inversion models for survey line 160

圖8 測區(qū)最小構(gòu)造反演電阻率平面(z=0 m)Fig.8 The plane view of the 1D minimum-structure inversion results (z=0 m)

2.3 LCI反演

根據(jù)最小構(gòu)造反演和Occam反演結(jié)果，研究區(qū)域電性結(jié)構(gòu)相對簡單，界面連續(xù)性良好，適合開展LCI反演。反演中，初始模型設(shè)置為三層層狀介質(zhì)，每層電阻率均為101.5Ω·m，第一、第二層的厚度均為100m；約束權(quán)重α和β依次設(shè)置為10、100和1000；目標(biāo)擬合差仍為3。在移動工作站(處理器為4核8線程、2.7GHz主頻)、Ubuntu操作系統(tǒng)條件下，單條測線一次反演耗時約為10min，最終擬合差均小于3。

圖9為160線LCI反演電阻率斷面， 其中玄武巖巖體(紅、黃和綠色區(qū)域)形態(tài)與上一節(jié)中單點(diǎn)反演結(jié)果基本一致，依然無法直接反映巖漿上涌通道信息。圖9中，隨著約束權(quán)重的增大，相鄰測點(diǎn)之間的玄武巖巖體電阻率差異縮小，以及玄武巖與沉積巖界面深度連續(xù)性不斷增強(qiáng)，并且反演獲取的玄武巖巖體厚度也略有增大，但在鉆孔ZK3處仍未達(dá)到400m。這種橫向連續(xù)性的增強(qiáng)是LCI反演的特點(diǎn)，也是與單點(diǎn)反演結(jié)果的主要區(qū)別。

圖9 160線LCI反演電阻率斷面Fig.9 The section view stitched from LCI inversion models for survey line 160

圖10為測區(qū)LCI反演結(jié)果俯視圖(α=β=1000),圖中同一測線相鄰測點(diǎn)之間的電阻率和層厚連續(xù)性良好，但測線之間相鄰測點(diǎn)的電阻率差異依然較大。

圖10 測區(qū)LCI反演結(jié)果俯視圖(z=0 m，α=β=1 000)Fig.10 The plane view of the LCI results (z=0 m, α=β=1 000)

2.4 SCI反演

SCI反演同樣適用于該研究區(qū)域。反演中，初始模型設(shè)置與上一節(jié)LCI反演相同；約束權(quán)重α和β依次設(shè)置為10、100和1000；目標(biāo)擬合差仍為3。整個測區(qū)一次反演耗時約為60min。圖11為160線SCI反演電阻率斷面，其與LCI反演結(jié)果(圖9)非常相似，同一條測線相鄰測點(diǎn)之間電阻率和層厚連續(xù)性隨約束權(quán)重α和β變化的規(guī)律也一致。

圖11 160線SCI反演電阻率斷面Fig.11 The section view stitched from SCI inversion models for survey line 160

圖12～圖14為不同約束權(quán)重下測區(qū)SCI反演結(jié)果俯視圖。如圖所示,隨著約束權(quán)重的增大，任一測點(diǎn)與四周相鄰測點(diǎn)的電阻率連續(xù)性不斷增強(qiáng)，目標(biāo)函數(shù)中空間約束項作用逐漸凸顯。

圖12 測區(qū)SCI反演結(jié)果俯視圖(z=0 m，α=β=10)Fig.12 The plane view of the SCI results (z=0 m, α=β=10)

圖13 測區(qū)SCI反演結(jié)果俯視圖(z=0 m，α=β=100)Fig.13 The plane view of the SCI results (z=0 m,α=β=100)

圖14 測區(qū)SCI反演結(jié)果俯視圖(z=0 m，α=β=1 000)Fig.14 The plane view of the SCI results (z=0 m, α=β=1 000)

3 結(jié)論

本文詳細(xì)敘述了筆者團(tuán)隊研發(fā)的地面回線源瞬變電磁法一維反演系統(tǒng)，并將其應(yīng)用于內(nèi)蒙古那仁寶力格煤田玄武巖巖體形態(tài)探測實(shí)例。該實(shí)例中，單點(diǎn)反演、LCI和SCI刻畫了玄武巖巖體除巖漿上涌通道區(qū)域之外的分布形態(tài)。下一步，將以上述一維反演結(jié)果為初始模型，開展該實(shí)例的三維反演研究，力爭獲得更加精確的玄武巖巖體形態(tài)。