張文波，張瑩，李建慧

(中國(guó)地質(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ù)擬合項(xiàng)，也包含模型正則化項(xiàng)[3-4]。Occam反演與最小構(gòu)造反演類(lèi)似，不同之處在于正則化因子的選?。篛ccam反演通過(guò)一元函數(shù)優(yōu)化來(lái)選取正則化因子，而最小構(gòu)造反演中正則化因子通常固定不變或者通過(guò)一些非常簡(jiǎn)單的算法選取[5]，因此，Occam反演更加穩(wěn)定，對(duì)正則化初始值依賴程度低，但其計(jì)算量偏大。LCI方法同時(shí)反演電阻率和層厚度，并在目標(biāo)函數(shù)中加入相鄰測(cè)點(diǎn)之間電阻率、層厚度和界面深度約束項(xiàng)[6]，因此單條測(cè)線所有測(cè)點(diǎn)數(shù)據(jù)需整體反演。SCI與LCI類(lèi)似，不僅考慮同一條測(cè)線相鄰測(cè)點(diǎn)之間的約束，還需考慮其他測(cè)線相鄰測(cè)點(diǎn)之間的約束[7]，因此測(cè)區(qū)內(nèi)所有測(cè)點(diǎn)數(shù)據(jù)需整體反演。

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

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

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

1 反演算法原理

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

圖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)

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

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

(2)

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

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

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

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

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

(4)

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

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

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

圖4 LCI和SCI參與約束測(cè)點(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)地勢(shì)平坦，新近系火山噴出巖比較發(fā)育，火山巖漿噴溢時(shí)的通道穿越煤層時(shí)不但對(duì)煤層產(chǎn)生破壞，而且對(duì)煤的變質(zhì)也會(huì)有不同程度的影響。為了保證今后煤田的順利開(kāi)采，瞬變電磁法被用于確定巖漿上涌通道，以及通道周邊玄武巖與圍巖的界面。

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

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

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

2.2 單點(diǎn)反演

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

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

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

圖8為測(cè)區(qū)最小構(gòu)造反演結(jié)果俯視圖，色標(biāo)與圖7保持一致。圖中，測(cè)線中心區(qū)域電阻率最大，色標(biāo)為桔紅色至紅色，測(cè)線兩端電阻率偏低，顯示為綠色甚至藍(lán)色。顏色的這種分布規(guī)律一定程度上反映了玄武巖厚度。另外，盡管多數(shù)測(cè)點(diǎn)與其四周相鄰測(cè)點(diǎn)的電阻率在同一個(gè)色標(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 測(cè)區(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)相對(duì)簡(jiǎn)單，界面連續(xù)性良好，適合開(kāi)展LCI反演。反演中，初始模型設(shè)置為三層層狀介質(zhì)，每層電阻率均為101.5Ω·m，第一、第二層的厚度均為100m；約束權(quán)重α和β依次設(shè)置為10、100和1000；目標(biāo)擬合差仍為3。在移動(dòng)工作站(處理器為4核8線程、2.7GHz主頻)、Ubuntu操作系統(tǒng)條件下，單條測(cè)線一次反演耗時(shí)約為10min，最終擬合差均小于3。

圖9為160線LCI反演電阻率斷面， 其中玄武巖巖體(紅、黃和綠色區(qū)域)形態(tài)與上一節(jié)中單點(diǎn)反演結(jié)果基本一致，依然無(wú)法直接反映巖漿上涌通道信息。圖9中，隨著約束權(quán)重的增大，相鄰測(cè)點(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為測(cè)區(qū)LCI反演結(jié)果俯視圖(α=β=1000),圖中同一測(cè)線相鄰測(cè)點(diǎn)之間的電阻率和層厚連續(xù)性良好，但測(cè)線之間相鄰測(cè)點(diǎn)的電阻率差異依然較大。

圖10 測(cè)區(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。整個(gè)測(cè)區(qū)一次反演耗時(shí)約為60min。圖11為160線SCI反演電阻率斷面，其與LCI反演結(jié)果(圖9)非常相似，同一條測(cè)線相鄰測(cè)點(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)重下測(cè)區(qū)SCI反演結(jié)果俯視圖。如圖所示,隨著約束權(quán)重的增大，任一測(cè)點(diǎn)與四周相鄰測(cè)點(diǎn)的電阻率連續(xù)性不斷增強(qiáng)，目標(biāo)函數(shù)中空間約束項(xiàng)作用逐漸凸顯。

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

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

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