無(wú)人機(jī)頻率域半航空電磁法三維反演

王堃鵬， 羅威 ， 曹輝， 王緒本， 藍(lán)星， 段長(zhǎng)生

1 成都理工大學(xué)地球物理學(xué)院， 成都 610059 2 中國(guó)地質(zhì)調(diào)查局成都地質(zhì)調(diào)查中心, 成都 610081 3 四川省冶勘設(shè)計(jì)集團(tuán)有限公司， 成都 610051 4 贛中南地質(zhì)礦產(chǎn)勘查研究院， 南昌 330029

0 引言

航空Z(yǔ)軸傾子法(ZTEM)正在逐漸成為一種快速普查淺層電阻率結(jié)構(gòu)的重要方法(Hübert et al.，2016；Wang et al.，2016；Lee et al.，2018)，該方法利用直升機(jī)在一定飛行高度測(cè)量天然場(chǎng)源的垂直磁場(chǎng)，具有測(cè)量效率高的優(yōu)點(diǎn)．ZTEM最大的優(yōu)勢(shì)是不需要直升機(jī)攜帶發(fā)射源，然而ZTEM觀測(cè)的天然場(chǎng)源數(shù)據(jù)易受到環(huán)境噪音的污染．因此，我們希望可以使用可控源替代天然場(chǎng)源．最常見的利用人工源替代天然場(chǎng)源的勘探法，就是可控源音頻大地電磁法(CSAMT)(林昌洪等，2012；Lin et al., 2018; He et al., 2019; Wang et al., 2019)．CSAMT采用有限長(zhǎng)線源產(chǎn)生一個(gè)人工電磁場(chǎng)，在強(qiáng)噪音環(huán)境下，通常被用來(lái)替代音頻大地電磁法(AMT)．

為了改善ZTEM，一種新的航空電磁系統(tǒng)在2015SEG會(huì)議被首次提出(Kang et al., 2015)，實(shí)際可以認(rèn)為是人工源頻率域半航空電磁法．設(shè)計(jì)者在早期計(jì)劃利用直升機(jī)測(cè)量垂直磁場(chǎng)(Nittinger et al., 2017)，地面放置發(fā)射源，而現(xiàn)在新的實(shí)驗(yàn)和分析正在使用無(wú)人機(jī)進(jìn)行(Zhou et al., 2016；Gao et al., 2019；Lin et al., 2019；Liu et al., 2020)．利用無(wú)人機(jī)測(cè)量空中數(shù)據(jù)已在多個(gè)領(lǐng)域進(jìn)行了嘗試(Koparan et al.，2018；Zhai et al.，2019；Kachroo et al.，2019)，而將無(wú)人機(jī)應(yīng)用于半航空電磁法同樣具有明顯的優(yōu)勢(shì)：①利用無(wú)人機(jī)可以測(cè)量來(lái)自人工源更強(qiáng)的垂直磁場(chǎng)；②無(wú)人機(jī)易于操作，可以懸停在空中更好的記錄數(shù)據(jù)，且發(fā)射頻率可以更低；③無(wú)人機(jī)不需要攜帶發(fā)射源．

目前無(wú)人機(jī)頻率域半航空電磁法主要還在硬件研發(fā)階段．Zhou等(2016)與Lin等(2019)介紹了他們的硬件系統(tǒng)，并初步對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了分析．Gao等(2019)描述了硬件開發(fā)中一種新的發(fā)射波形．Liu等(2020)利用 COMSOL Multiphysics分析軟件進(jìn)一步研究了數(shù)據(jù)分布特征．可以看出，無(wú)人機(jī)頻率域半航空電磁法還有待進(jìn)一步發(fā)展，但相關(guān)成果已經(jīng)顯示出該領(lǐng)域廣闊的前景．

本文對(duì)Liu等(2020)的系統(tǒng)進(jìn)行了調(diào)整，兩個(gè)水平的接地雙極源放置在不同的坐標(biāo)軸上．為了更好的提高無(wú)人機(jī)半航空電磁法的反演效果，我們將它與大地電磁法進(jìn)行聯(lián)合反演，并引入交叉梯度實(shí)現(xiàn)結(jié)構(gòu)約束．首先，大地電磁具有更低的頻率，適當(dāng)?shù)脑黾哟蟮仉姶艤y(cè)深點(diǎn)，可以極大的提升勘探深度．利用大地電磁增加勘探深度已在其它領(lǐng)域有了諸多研究，比如Lee等(2018)實(shí)現(xiàn)了大地電磁與ZTEM的聯(lián)合反演，Amatyakul et al.(2017)研究了大地電磁與直流電法的聯(lián)合反演．其次，交叉梯度結(jié)構(gòu)約束可以進(jìn)一步提高電磁法的分辨能力，交叉梯度理論是Gallardo和Meju(2003)首次提出，并且他們很快就實(shí)現(xiàn)了直流電法與地震的二維聯(lián)合反演(Gallardo and Meju，2004)，取得了不錯(cuò)的效果．從此交叉梯度逐漸成為提高反演分辨率的重要方法，彭淼等(2013)實(shí)現(xiàn)了大地電磁與地震走時(shí)資料三維聯(lián)合反演，Wang等(2017)實(shí)現(xiàn)了CSAMT與磁法的二維聯(lián)合反演，閆政文等(2020)與張镕哲等(2019)實(shí)現(xiàn)了多個(gè)方法的聯(lián)合反演，吳萍萍等(2020)實(shí)現(xiàn)了電阻率法和背景噪聲法三維聯(lián)合反演．從這些研究可以看出，交叉梯度可以有效提高反演分辨率．

為了實(shí)現(xiàn)上述目標(biāo)，本文正演采用交錯(cuò)網(wǎng)格有限差分，反演采用有限內(nèi)存擬牛頓法(LBFGS)，對(duì)于結(jié)構(gòu)約束我們使用的方式不需要進(jìn)行泰勒展開與近似，可以快速實(shí)現(xiàn)基于交叉梯度的結(jié)構(gòu)約束方案．最后，本文通過(guò)建立四個(gè)理論模型驗(yàn)證三維反演的效果．

1 接收系統(tǒng)與正演模擬

1.1 無(wú)人機(jī)半航空電磁法裝置

圖1 發(fā)射與接收系統(tǒng)示意圖Fig.1 Diagram of the transmission and receiving system

圖1展示了本文的發(fā)射與接收系統(tǒng)，兩個(gè)正交的接地雙極源被放置在不同的坐標(biāo)軸上．發(fā)射源1放置在Y軸上，并平行于X軸，發(fā)射源2放置在X軸上，且平行于Y軸.當(dāng)發(fā)射源1工作時(shí)，關(guān)閉發(fā)射源2，無(wú)人機(jī)懸停在空中采集垂直磁感應(yīng)強(qiáng)度Bz1.當(dāng)發(fā)射源2工作時(shí)，關(guān)閉發(fā)射源1，無(wú)人機(jī)懸停在空中采集垂直磁感應(yīng)強(qiáng)度Bz2.圖1所示的系統(tǒng)與Liu 等(2020)有一定的不同，此外我們?cè)趫D1中還顯示了地面的大地電磁觀測(cè)臺(tái)站，目的就是要進(jìn)一步拓展整個(gè)系統(tǒng)的勘探深度.

1.2 正演模擬

對(duì)于圖1的無(wú)人機(jī)系統(tǒng)，這里采用一次場(chǎng)與二次場(chǎng)分離的辦法實(shí)現(xiàn)正演，二次場(chǎng)方程如下：

(1)

其中，ω為圓頻率，μ0為真空磁導(dǎo)率，σ為介質(zhì)電導(dǎo)率，σp為背景電導(dǎo)率，Es為二次場(chǎng)，Ep為背景場(chǎng)．針對(duì)大地電磁法，使用總場(chǎng)法計(jì)算：

(2)

對(duì)于方程(1)與方程(2)，我們采用交錯(cuò)網(wǎng)格有限差分進(jìn)行模擬計(jì)算．方程(1)的背景場(chǎng)采用解析解獲得，邊界條件設(shè)為電場(chǎng)的切向分量為零．方程(2)的邊界條件來(lái)自于二維大地電磁正演，計(jì)算速度快，具體的離散過(guò)程可以參考作者之前的論文(Wang et al.，2019；羅威等，2019)．

2 反演實(shí)現(xiàn)

2.1 無(wú)人機(jī)反演目標(biāo)函數(shù)及梯度式

本文采用LBFGS法進(jìn)行反演迭代計(jì)算，設(shè)目標(biāo)函數(shù)為

φ=φd_UAV+λφm

(3)

本文采用仿射線性參數(shù)變換(Egbert and Kelbert，2012；Kelbert et al.，2014)對(duì)真實(shí)模型參數(shù)進(jìn)行如下轉(zhuǎn)換：

(4)

在上式變換下，最終的真實(shí)模型參數(shù)可以用以下式子獲得：

(5)

在參數(shù)變換策略下，目標(biāo)函數(shù)(3)的梯度式為

(6)

針對(duì)式(6)所示的梯度計(jì)算，本文采用“擬正演”方案實(shí)現(xiàn)(Newman and Alumbaugh，2000；Commer and Newman，2008)．

2.2 與大地電磁聯(lián)合反演的目標(biāo)函數(shù)及梯度式

為了提高無(wú)人機(jī)頻率域半航空電磁法的反演深度，我們?cè)诜囱葜腥谌肓舜蟮仉姶欧?，設(shè)聯(lián)合反演目標(biāo)函數(shù)為

φ=φd_UAV+φd_MT+λφm

(7)

其中，MT代表大地電磁數(shù)據(jù)，其余參數(shù)與(3)式一致．

同樣在參數(shù)變換的策略下，目標(biāo)函數(shù)(7)的梯度表達(dá)式為

(8)

上式中的梯度計(jì)算，這里依然采用“擬正演”方案實(shí)現(xiàn)．

2.3 基于交叉梯度的結(jié)構(gòu)耦合

在2.2節(jié)聯(lián)合反演目標(biāo)函數(shù)的基礎(chǔ)上，我們希望可以進(jìn)一步提高對(duì)先驗(yàn)結(jié)構(gòu)信息的利用．基于此，本文提出在目標(biāo)函數(shù)(7)的基礎(chǔ)上進(jìn)一步添加交叉梯度項(xiàng)：

φ=φd_UAV+φd_MT+λφm+λcgφcg，

(9)

其中，φcg為交叉梯度項(xiàng)，λcg為交叉梯度項(xiàng)權(quán)重．

本文的交叉梯度項(xiàng)φcg定義如下：

(10)

上式中，mv為外部輸入的已知結(jié)構(gòu)信息，如速度結(jié)構(gòu)．將上式按三個(gè)方向進(jìn)一步展開，有

(11)

根據(jù)閆政文等(2020)，式(11)中tx,ty,tz三個(gè)向量的任意元素可以表述為

(12)

其中，m與mv分別為某個(gè)電阻率參數(shù)和約束模型參數(shù)．

對(duì)于某個(gè)具體的tx(i,j,z)計(jì)算，這里以圖2來(lái)說(shuō)明，將tx(i,j,z)以中心差分離散：

(13)

圖2 交叉梯度離散示意圖Fig.2 The discrete diagram of cross-gradient

將式(13)以矩陣形式描述，則有

(14)

其中，a11、a22、a33、a44、a55形式如下：

a33=-a11-a22-a44-a55.

(15)

在圖2中的內(nèi)部單元循環(huán)，可以完成tx向量所有元素的計(jì)算．最終的tx向量是由式(14)的矩陣?yán)奂有纬?，拓展到整個(gè)三維模型單元為

tx=Wx_cgmx.

(16)

同樣地，ty,tz以矩陣形式表述為

ty=Wy_cgm,

tz=Wz_cgm.

(17)

在式(16)與(17)中，矩陣Wx_cg,Wy_cg,Wz_cg包含了網(wǎng)格信息以及約束模型，因此交叉梯度以類似于模型協(xié)方差矩陣的形式對(duì)電阻率參數(shù)進(jìn)行了結(jié)構(gòu)約束．在式(16)與(17)的基礎(chǔ)上，交叉梯度項(xiàng)φcg最終的矩陣形式為

φcg=(Wx_cgm)T(Wx_cgm)+(Wy_cgm)T(Wy_cgm)

+(Wz_cgm)T(Wz_cgm).

(18)

為了配合變換后的目標(biāo)函數(shù)，最終含交叉梯度的總目標(biāo)函數(shù)梯度式，有如下描述：

(19)

最終在目標(biāo)函數(shù)(3)(7)(9)及相應(yīng)的梯度表達(dá)式下，我們使用LBFGS(Nocedal and Wright，2006)完成上述反演．

3 半航空頻率域電磁法正演分析

本節(jié)建立了如圖3所示的正演模型，背景電阻率100 Ωm，低阻模型電阻率為10 Ωm.模型的頂面埋深為120 m，整體尺寸為1600 m×1600 m×200 m.根據(jù)圖1所示的系統(tǒng)，這里將發(fā)射源1放置在y軸負(fù)半軸(x=0 m,y=-10 km,z=0 m)，而發(fā)射源2放置在x軸負(fù)半軸(x=-10 km,y=0 m,z=0 m).整個(gè)測(cè)量范圍在-3920 m到3920 m，在X和Y方向每隔160 m一個(gè)點(diǎn)(圖4).垂直磁感應(yīng)強(qiáng)度Bz1來(lái)自發(fā)射源1，垂直磁感應(yīng)強(qiáng)度Bz2來(lái)自發(fā)射源2，采集高度均在100 m空中，測(cè)點(diǎn)分布2500個(gè).

圖3 三維模型俯視圖Fig.3 The top view of 3-D model

圖4 測(cè)點(diǎn)分布圖，所有測(cè)點(diǎn)均在空中100 m高度Fig.4 The distribution of measuring points．All the points are at a height of 100 m in the air

計(jì)算剖分網(wǎng)格為80×80×39，這里給出50 Hz和150 Hz的計(jì)算結(jié)果(圖5)，圖5中的黑色方框?yàn)楫惓ｓw在地面的投影.從圖5的結(jié)果可以看到，Bz1對(duì)Y方向的邊界靈敏，Bz2對(duì)X方向的邊界靈敏.圖5的結(jié)果顯示，應(yīng)當(dāng)使用兩個(gè)正交發(fā)射源分別采集Bz1和Bz2可以獲得最佳的分辨率，我們會(huì)在下一節(jié)反演中進(jìn)一步的證實(shí)這個(gè)結(jié)論.

圖5 無(wú)人機(jī)測(cè)量的Bz1與Bz2振幅Fig.5 The magnitude of Bz1 and Bz2 recorded by UAV

4 理論模型反演分析

4.1 無(wú)人機(jī)半航空頻率域電磁法反演

本節(jié)建立了如圖6所示的理論模型，該模型是為了進(jìn)一步證明第3節(jié)正演中得出的結(jié)論.圖6模型的背景電阻率為100 Ωm，四個(gè)高阻異常體模型為1000 Ωm.異常體的頂面埋深為120 m，尺寸為960 m×960 m×200 m，四個(gè)頻率(300，150，80，50 Hz)被用于反演測(cè)試.

圖6 反演模型1俯視圖Fig.6 The top view of the first inversion model

測(cè)點(diǎn)位置及范圍仍然沿用圖4所示的分布，兩個(gè)水平正交的雙極源被分別放置在X軸負(fù)半軸以及Y軸負(fù)半軸，源1的具體坐標(biāo)為(x=0 m,y=-10 km,z=0 m)，源2為(x=-10 km,y=0 m,z=0 m).正演及反演網(wǎng)格為80×80×39，正演數(shù)據(jù)添加2%的高斯隨機(jī)噪音.Bz1和Bz2數(shù)據(jù)的反演權(quán)重為各自振幅的2%，反演數(shù)據(jù)為Bz1和Bz2的實(shí)部與虛部，初始模型為100 Ωm均勻半空間.為了更好地展示Bz1和Bz2在反演中的作用，我們選擇三組反演數(shù)據(jù)進(jìn)行測(cè)試：①僅使用Bz1反演；②僅使用Bz2反演；③Bz1和Bz2同時(shí)參與反演.

最終的反演結(jié)果如圖8和9所示，圖7給出了三種數(shù)據(jù)體的反演擬合差曲線，Bz1、Bz2、Bz1+Bz2擬合差分別為0.993，0.979，0.997，整體反演呈現(xiàn)出穩(wěn)定收斂，證明LBFGS適用于無(wú)人機(jī)頻率域半航空電磁法反演.為了更好地展示細(xì)節(jié)，這里縮小了成圖范圍.在圖5的正演響應(yīng)特征中，我們可以看到Bz1與Bz2在X與Y方向的顯著差異，Bz1對(duì)Y方向邊界靈敏，而Bz2對(duì)X方向邊界靈敏，這個(gè)現(xiàn)象在圖7和圖8的反演結(jié)果中也再一次體現(xiàn)，而Bz1和Bz2的組合反演得到了最佳的反演效果.對(duì)于無(wú)人機(jī)而言，使用兩個(gè)發(fā)射源的工作量?jī)H僅是多設(shè)置了一個(gè)發(fā)射源.盡管顯著增加了無(wú)人機(jī)測(cè)量的工作量，但相比于地面測(cè)量而言，仍然有極高的效率.因此，我們認(rèn)為實(shí)際情況下，應(yīng)盡量使用兩個(gè)發(fā)射源，以增加對(duì)復(fù)雜異常體的識(shí)別能力.

圖7 擬合差曲線Fig.7 The rms curves

圖8 水平切片圖第一列為真實(shí)模型，第二、三、四列分別為Bz1、Bz2、Bz1+Bz2反演結(jié)果.Fig.8 The horizontal sliceThe first column is the true model. The second、third、fourth column is the inversion results of Bz1、Bz2、Bz1+Bz2, respectively.

圖9 垂直切片圖第一列為真實(shí)模型，第二、三、四列分別為Bz1、Bz2、Bz1+Bz2反演結(jié)果.Fig.9 The vertical sliceThe first column is the true model. The second、third、fourth column is the inversion results of Bz1、Bz2、Bz1+Bz2, respectively.

4.2 反演抗噪測(cè)試

在本節(jié)我們將利用一個(gè)低阻棱柱體模型進(jìn)行反演抗噪測(cè)試，在正演數(shù)據(jù)中添加不同程度的噪音，進(jìn)一步驗(yàn)證無(wú)人機(jī)數(shù)據(jù)反演的穩(wěn)定性和可靠性.如圖10所示，本節(jié)的低阻模型尺寸為1280 m×1280 m×200 m，埋深120 m.背景電阻率為100 Ωm，低阻異常體電阻率為10 Ωm.反演頻率、發(fā)射源位置及測(cè)點(diǎn)分布與前一節(jié)一致，本節(jié)分別對(duì)正演響應(yīng)添加3%、5%、10%、20%、35%、50%的隨機(jī)噪音，所有反演的數(shù)據(jù)權(quán)重都選擇為場(chǎng)值的3%.

圖10 反演模型2俯視圖Fig.10 The top view of the second inversion model

反演擬合差曲線如圖11所示，對(duì)于3%的噪音而言，由于數(shù)據(jù)權(quán)重也與之匹配，僅迭代25次就收斂，擬合差為0.984，為了防止出現(xiàn)嚴(yán)重過(guò)擬合的現(xiàn)象，余下反演都限制了最大迭代次數(shù)為25次.最終的反演結(jié)果如圖12所示，從圖中的結(jié)果可以看到無(wú)人機(jī)數(shù)據(jù)隨著噪音的逐漸增加，不僅收斂越來(lái)越難，反演結(jié)果也隨之出現(xiàn)愈來(lái)愈多的假異常.所有的反演對(duì)于異常體都具有非常明顯的反應(yīng)，并且隨機(jī)噪音增加至20%時(shí)，整體的反演仍然較為理想，這說(shuō)明本文的無(wú)人機(jī)數(shù)據(jù)反演具有較強(qiáng)的抗噪性.需要指出的時(shí)，實(shí)際情況的地質(zhì)結(jié)構(gòu)復(fù)雜，噪音水平仍然需要做到最大程度的降低.

圖11 擬合差曲線Fig.11 The rms curves

圖12 不同噪音的反演結(jié)果Fig.12 The inversion results of different noise

4.3 聯(lián)合大地電磁法反演

在4.1節(jié)我們證明了同時(shí)使用Bz1和Bz2反演效果最佳，然而對(duì)于無(wú)人機(jī)來(lái)說(shuō)接收更低的頻率仍然十分困難，因此僅靠無(wú)人機(jī)還無(wú)法獲得更深的電阻率結(jié)構(gòu).獲得更低頻率的常用做法，是與大地電磁進(jìn)行聯(lián)合反演，這樣可以使用頻率更低的觀測(cè)數(shù)據(jù).

為了研究無(wú)人機(jī)數(shù)據(jù)與大地電磁法的聯(lián)合反演效果，我們首先建立了如圖13所示的模型.背景電阻率為100 Ωm，淺部和深部的高阻異常體均為1000 Ωm，低阻異常體為10 Ωm.淺部的高阻異常體頂面埋深為120 m，深部的高、低阻異常體頂面埋深均為920 m.淺部低阻異常體的尺寸為1600 m×1600 m×200 m，兩個(gè)深部異常體的尺寸均為3200 m×2400 m×1200 m.

圖13 反演模型3Fig.13 The third inversion model

無(wú)人機(jī)的觀測(cè)系統(tǒng)與前一節(jié)一致，測(cè)點(diǎn)分布如圖4所示.大地電磁的觀測(cè)點(diǎn)遠(yuǎn)遠(yuǎn)小于無(wú)人機(jī)觀測(cè)點(diǎn)數(shù)目，其分布如圖14，共400個(gè)大地電磁臺(tái)站.因?yàn)榇蟮仉姶旁谶@里關(guān)注的是更深的結(jié)構(gòu)，因此臺(tái)站整體分布相對(duì)稀疏，圖14中大地電磁臺(tái)站在X和Y方向每隔400 m一個(gè).

圖14 400個(gè)大地電磁測(cè)點(diǎn)分布圖Fig.14 The distribution of 400 MT measuring points

在本節(jié)測(cè)試中，我們使用了一個(gè)較細(xì)的網(wǎng)格進(jìn)行正演計(jì)算，網(wǎng)格剖分為136×136×60.無(wú)人機(jī)的計(jì)算頻率為300，150，80，50 Hz，反演數(shù)據(jù)依然是實(shí)部和虛部.大地電磁的計(jì)算頻率為50，10，1，0.5，0.1 Hz，共5個(gè)頻率.為了驗(yàn)證反演的穩(wěn)健性，反演采用稀疏網(wǎng)格，整體剖分為80×80×50，正演數(shù)據(jù)加入了3%的高斯隨機(jī)噪音，反演初始模型為100 Ωm的均勻半空間.

為了體現(xiàn)單獨(dú)反演與聯(lián)合反演的差異，我們首先使用無(wú)人機(jī)數(shù)據(jù)進(jìn)行了反演測(cè)試，反演結(jié)果如圖15所示.從圖中結(jié)果可以看到，單獨(dú)的無(wú)人機(jī)數(shù)據(jù)反演可以較好地恢復(fù)淺部結(jié)構(gòu)，但深部的兩個(gè)高低阻異常體都未能很好地反應(yīng)出來(lái).隨后，我們?cè)趫D16中給出了僅僅依靠大地電磁反演出來(lái)的結(jié)果，圖16的深部相較于圖15而言有了明顯的改觀，但此時(shí)圖16的淺部異常體效果不如圖15.最后，我們將無(wú)人機(jī)數(shù)據(jù)與大地電磁數(shù)據(jù)進(jìn)行了聯(lián)合反演，結(jié)果如圖17所示，從圖中的結(jié)果可以看到，聯(lián)合反演具有相對(duì)最佳的效果，既能反應(yīng)淺部構(gòu)造，同時(shí)也能獲得深部信息.圖18展示了三個(gè)反演的擬合差曲線，從中可以看出三個(gè)反演都非常穩(wěn)定.

圖15 無(wú)人機(jī)數(shù)據(jù)單獨(dú)反演Fig.15 The single inversion of UAV data

圖16 大地電磁單獨(dú)反演Fig.16 The single inversion of MT

本節(jié)的反演實(shí)驗(yàn)為實(shí)際勘探給出了一個(gè)可能的方案，即對(duì)于淺部結(jié)構(gòu)，借助無(wú)人機(jī)的快捷，可以進(jìn)行較為密集的觀測(cè).而深部信息的獲得，則可在地面補(bǔ)充布置較為稀疏的大地電磁測(cè)點(diǎn)，這樣可以相對(duì)快速地獲得淺部和深部信息，提高野外工作效率.

圖17 聯(lián)合反演Fig.17 The joint inversion

圖18 擬合差曲線Fig.18 The rms curves

4.4 基于交叉梯度結(jié)構(gòu)約束的聯(lián)合反演

在4.3節(jié)，我們進(jìn)一步證明了無(wú)人機(jī)數(shù)據(jù)與大地電磁聯(lián)合反演的有效性，可同時(shí)獲得淺部和深部結(jié)構(gòu).然而，對(duì)于實(shí)際情況而言，面對(duì)更復(fù)雜的異常體，有可能這種數(shù)據(jù)體仍顯不足，對(duì)真實(shí)電阻率模型的分辨能力仍然有限.因此，利用某些已知的結(jié)構(gòu)(如速度結(jié)構(gòu))，約束電磁法反演逐漸成為一種流行的方式.在這種背景下，交叉梯度理論被提出和廣泛使用，并且已經(jīng)被證明可以有效實(shí)施結(jié)構(gòu)約束反演.

在4.3節(jié)的支撐下，本節(jié)嘗試進(jìn)行基于交叉梯度結(jié)構(gòu)約束的無(wú)人機(jī)頻率域半航空電磁法與大地電磁法的聯(lián)合反演.首先建立圖19所示的電阻率模型與速度模型，這里的速度模型與電阻率模型具有相同的結(jié)構(gòu).淺部異常體的電阻率和速度分別為500 Ωm及8000 m·s-1，深部異常體的電阻率和速度分別為10 Ωm及2000 m·s-1，背景電阻率和背景速度分別為100 Ωm及6000 m·s-1.無(wú)人機(jī)的觀測(cè)系統(tǒng)與前一節(jié)一致，測(cè)點(diǎn)分布如圖4所示，大地電磁的觀測(cè)點(diǎn)與圖14一致.

圖19 反演模型4與速度模型Fig.19 The fourth inversion model and the velocity model

在本節(jié)測(cè)試中，反演網(wǎng)格剖分為80×80×50.無(wú)人機(jī)的計(jì)算頻率為300，150，80，50 Hz，反演數(shù)據(jù)為垂直磁感應(yīng)強(qiáng)度的實(shí)部與虛部.大地電磁的計(jì)算頻率為50，10，1，0.5，0.1 Hz.正演數(shù)據(jù)加入了2%的高斯隨機(jī)噪音，反演初始模型為100 Ωm的均勻半空間.

最終的反演結(jié)果如圖20所示，首先對(duì)于淺部異常體而言，由于無(wú)人機(jī)數(shù)據(jù)的密集性且異常體形態(tài)相對(duì)簡(jiǎn)單，交叉梯度并沒有對(duì)淺部構(gòu)造取得更明顯的約束效果.然而，面對(duì)深部的低阻異常體，無(wú)交叉梯度約束的反演，對(duì)構(gòu)造形態(tài)的準(zhǔn)確識(shí)別開始變差，尤其是對(duì)低阻中間的凹槽恢復(fù)能力較弱.使用交叉梯度項(xiàng)進(jìn)行結(jié)構(gòu)約束后，深部結(jié)構(gòu)開始有了明顯改善，這說(shuō)明交叉梯度項(xiàng)在聯(lián)合反演中有效，同時(shí)也證明了本文程序的正確性.圖21展示了無(wú)交叉梯度與使用交叉梯度的擬合差曲線，從圖中可以看到二者在局部迭代的差異，這也同時(shí)說(shuō)明交叉梯度項(xiàng)確實(shí)是在影響整體反演的走勢(shì).

圖20 反演結(jié)果(第一列無(wú)交叉梯度項(xiàng)，第二列使用交叉梯度項(xiàng))Fig.20 The inversion results (the first column has no cross-gradient term, the second column has cross-gradient term)

圖21 擬合差曲線Fig.21 The rms curves

本節(jié)的反演初步嘗試了使用結(jié)構(gòu)約束的聯(lián)合反演，對(duì)于實(shí)際勘探而言，下一步還有很多工作需要驗(yàn)證，比如無(wú)人機(jī)數(shù)據(jù)的準(zhǔn)確性，聯(lián)合反演的實(shí)用性，結(jié)構(gòu)約束的可靠性等等.但總體而言，我們認(rèn)為無(wú)人機(jī)勘探潛力巨大，值得更深入的研究.

5 結(jié)論

無(wú)人機(jī)探測(cè)效率高，在研究淺層電阻率結(jié)構(gòu)方面有很大的潛力.在本文中，我們?yōu)闊o(wú)人機(jī)頻率域半航空電磁法開發(fā)了一種有效的三維反演算法.在正演中，我們使用了交錯(cuò)網(wǎng)格有限差分法，該方法在三維電磁法中很容易實(shí)現(xiàn).在反演部分，我們采用了LBFGS法，它具有較高的穩(wěn)定性和效率，不需要計(jì)算和存儲(chǔ)巨大的靈敏度矩陣.

為了進(jìn)一步增強(qiáng)無(wú)人機(jī)勘探的反演能力，我們?nèi)谌肓伺c大地電磁的聯(lián)合反演，以及基于交叉梯度方案的結(jié)構(gòu)約束.在正演模擬和第一個(gè)反演測(cè)試中，我們比較了垂直磁感應(yīng)強(qiáng)度Bz1和Bz2之間的分辨率差異，并證明同時(shí)使用Bz1和Bz2是最理想的，因?yàn)锽z1對(duì)Y方向邊界敏感，而Bz2對(duì)X方向邊界敏感.在第二個(gè)反演測(cè)試中，我們對(duì)正演響應(yīng)添加了不同程度的噪音，結(jié)果顯示無(wú)人機(jī)數(shù)據(jù)的反演具有較強(qiáng)的抗干擾能力，這在一定程度上增加了未來(lái)實(shí)測(cè)數(shù)據(jù)反演的可靠性.在第三個(gè)反演測(cè)試中，我們分別進(jìn)行了三個(gè)反演試驗(yàn)，研究結(jié)果表明：①無(wú)人機(jī)的單獨(dú)反演可以很好地恢復(fù)淺層結(jié)構(gòu)；②大地電磁的單獨(dú)反演可以很好地恢復(fù)深層結(jié)構(gòu)；③無(wú)人機(jī)與大地電磁的聯(lián)合反演可以同時(shí)獲得淺部和深部的結(jié)構(gòu).最后，在第四個(gè)模型中，我們依靠交叉梯度實(shí)現(xiàn)了針對(duì)聯(lián)合反演的結(jié)構(gòu)約束.結(jié)果證明，在實(shí)際情況下若能夠從其它地球物理方法獲得可靠的結(jié)構(gòu)信息，交叉梯度可以顯著改善電磁法反演對(duì)邊界的識(shí)別.