無人機頻率域半航空電磁法三維反演

王堃鵬， 羅威 ， 曹輝， 王緒本， 藍星， 段長生

1 成都理工大學地球物理學院， 成都 610059 2 中國地質調查局成都地質調查中心, 成都 610081 3 四川省冶勘設計集團有限公司， 成都 610051 4 贛中南地質礦產勘查研究院， 南昌 330029

0 引言

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

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

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

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

為了實現上述目標，本文正演采用交錯網格有限差分，反演采用有限內存擬牛頓法(LBFGS)，對于結構約束我們使用的方式不需要進行泰勒展開與近似，可以快速實現基于交叉梯度的結構約束方案．最后，本文通過建立四個理論模型驗證三維反演的效果．

1 接收系統與正演模擬

1.1 無人機半航空電磁法裝置

圖1 發射與接收系統示意圖Fig.1 Diagram of the transmission and receiving system

圖1展示了本文的發射與接收系統，兩個正交的接地雙極源被放置在不同的坐標軸上．發射源1放置在Y軸上，并平行于X軸，發射源2放置在X軸上，且平行于Y軸.當發射源1工作時，關閉發射源2，無人機懸停在空中采集垂直磁感應強度Bz1.當發射源2工作時，關閉發射源1，無人機懸停在空中采集垂直磁感應強度Bz2.圖1所示的系統與Liu 等(2020)有一定的不同，此外我們在圖1中還顯示了地面的大地電磁觀測臺站，目的就是要進一步拓展整個系統的勘探深度.

1.2 正演模擬

對于圖1的無人機系統，這里采用一次場與二次場分離的辦法實現正演，二次場方程如下：

(1)

其中，ω為圓頻率，μ0為真空磁導率，σ為介質電導率，σp為背景電導率，Es為二次場，Ep為背景場．針對大地電磁法，使用總場法計算：

(2)

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

2 反演實現

2.1 無人機反演目標函數及梯度式

本文采用LBFGS法進行反演迭代計算，設目標函數為

φ=φd_UAV+λφm

(3)

本文采用仿射線性參數變換(Egbert and Kelbert，2012；Kelbert et al.，2014)對真實模型參數進行如下轉換：

(4)

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

(5)

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

(6)

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

2.2 與大地電磁聯合反演的目標函數及梯度式

為了提高無人機頻率域半航空電磁法的反演深度，我們在反演中融入了大地電磁法，設聯合反演目標函數為

φ=φd_UAV+φd_MT+λφm

(7)

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

同樣在參數變換的策略下，目標函數(7)的梯度表達式為

(8)

上式中的梯度計算，這里依然采用“擬正演”方案實現．

2.3 基于交叉梯度的結構耦合

在2.2節聯合反演目標函數的基礎上，我們希望可以進一步提高對先驗結構信息的利用．基于此，本文提出在目標函數(7)的基礎上進一步添加交叉梯度項：

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

(9)

其中，φcg為交叉梯度項，λcg為交叉梯度項權重．

本文的交叉梯度項φcg定義如下：

(10)

上式中，mv為外部輸入的已知結構信息，如速度結構．將上式按三個方向進一步展開，有

(11)

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

(12)

其中，m與mv分別為某個電阻率參數和約束模型參數．

對于某個具體的tx(i,j,z)計算，這里以圖2來說明，將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中的內部單元循環，可以完成tx向量所有元素的計算．最終的tx向量是由式(14)的矩陣累加形成，拓展到整個三維模型單元為

tx=Wx_cgmx.

(16)

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

ty=Wy_cgm,

tz=Wz_cgm.

(17)

在式(16)與(17)中，矩陣Wx_cg,Wy_cg,Wz_cg包含了網格信息以及約束模型，因此交叉梯度以類似于模型協方差矩陣的形式對電阻率參數進行了結構約束．在式(16)與(17)的基礎上，交叉梯度項φcg最終的矩陣形式為

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

+(Wz_cgm)T(Wz_cgm).

(18)

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

(19)

最終在目標函數(3)(7)(9)及相應的梯度表達式下，我們使用LBFGS(Nocedal and Wright，2006)完成上述反演．

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

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

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

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

計算剖分網格為80×80×39，這里給出50 Hz和150 Hz的計算結果(圖5)，圖5中的黑色方框為異常體在地面的投影.從圖5的結果可以看到，Bz1對Y方向的邊界靈敏，Bz2對X方向的邊界靈敏.圖5的結果顯示，應當使用兩個正交發射源分別采集Bz1和Bz2可以獲得最佳的分辨率，我們會在下一節反演中進一步的證實這個結論.

圖5 無人機測量的Bz1與Bz2振幅Fig.5 The magnitude of Bz1 and Bz2 recorded by UAV

4 理論模型反演分析

4.1 無人機半航空頻率域電磁法反演

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

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

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

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

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

圖8 水平切片圖第一列為真實模型，第二、三、四列分別為Bz1、Bz2、Bz1+Bz2反演結果.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 垂直切片圖第一列為真實模型，第二、三、四列分別為Bz1、Bz2、Bz1+Bz2反演結果.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 反演抗噪測試

在本節我們將利用一個低阻棱柱體模型進行反演抗噪測試，在正演數據中添加不同程度的噪音，進一步驗證無人機數據反演的穩定性和可靠性.如圖10所示，本節的低阻模型尺寸為1280 m×1280 m×200 m，埋深120 m.背景電阻率為100 Ωm，低阻異常體電阻率為10 Ωm.反演頻率、發射源位置及測點分布與前一節一致，本節分別對正演響應添加3%、5%、10%、20%、35%、50%的隨機噪音，所有反演的數據權重都選擇為場值的3%.

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

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

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

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

4.3 聯合大地電磁法反演

在4.1節我們證明了同時使用Bz1和Bz2反演效果最佳，然而對于無人機來說接收更低的頻率仍然十分困難，因此僅靠無人機還無法獲得更深的電阻率結構.獲得更低頻率的常用做法，是與大地電磁進行聯合反演，這樣可以使用頻率更低的觀測數據.

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

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

無人機的觀測系統與前一節一致，測點分布如圖4所示.大地電磁的觀測點遠遠小于無人機觀測點數目，其分布如圖14，共400個大地電磁臺站.因為大地電磁在這里關注的是更深的結構，因此臺站整體分布相對稀疏，圖14中大地電磁臺站在X和Y方向每隔400 m一個.

圖14 400個大地電磁測點分布圖Fig.14 The distribution of 400 MT measuring points

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

為了體現單獨反演與聯合反演的差異，我們首先使用無人機數據進行了反演測試，反演結果如圖15所示.從圖中結果可以看到，單獨的無人機數據反演可以較好地恢復淺部結構，但深部的兩個高低阻異常體都未能很好地反應出來.隨后，我們在圖16中給出了僅僅依靠大地電磁反演出來的結果，圖16的深部相較于圖15而言有了明顯的改觀，但此時圖16的淺部異常體效果不如圖15.最后，我們將無人機數據與大地電磁數據進行了聯合反演，結果如圖17所示，從圖中的結果可以看到，聯合反演具有相對最佳的效果，既能反應淺部構造，同時也能獲得深部信息.圖18展示了三個反演的擬合差曲線，從中可以看出三個反演都非常穩定.

圖15 無人機數據單獨反演Fig.15 The single inversion of UAV data

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

本節的反演實驗為實際勘探給出了一個可能的方案，即對于淺部結構，借助無人機的快捷，可以進行較為密集的觀測.而深部信息的獲得，則可在地面補充布置較為稀疏的大地電磁測點，這樣可以相對快速地獲得淺部和深部信息，提高野外工作效率.

圖17 聯合反演Fig.17 The joint inversion

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

4.4 基于交叉梯度結構約束的聯合反演

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

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

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

在本節測試中，反演網格剖分為80×80×50.無人機的計算頻率為300，150，80，50 Hz，反演數據為垂直磁感應強度的實部與虛部.大地電磁的計算頻率為50，10，1，0.5，0.1 Hz.正演數據加入了2%的高斯隨機噪音，反演初始模型為100 Ωm的均勻半空間.

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

圖20 反演結果(第一列無交叉梯度項，第二列使用交叉梯度項)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

本節的反演初步嘗試了使用結構約束的聯合反演，對于實際勘探而言，下一步還有很多工作需要驗證，比如無人機數據的準確性，聯合反演的實用性，結構約束的可靠性等等.但總體而言，我們認為無人機勘探潛力巨大，值得更深入的研究.

5 結論

無人機探測效率高，在研究淺層電阻率結構方面有很大的潛力.在本文中，我們為無人機頻率域半航空電磁法開發了一種有效的三維反演算法.在正演中，我們使用了交錯網格有限差分法，該方法在三維電磁法中很容易實現.在反演部分，我們采用了LBFGS法，它具有較高的穩定性和效率，不需要計算和存儲巨大的靈敏度矩陣.

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