時間域航空電磁橫向分段組合約束反演

林興元 陸從德* 鐘炳輝 范崇霄

(①成都理工大學地球物理學院，四川成都 610059； ②成都理工大學信息科學與技術學院，四川成都 610059)

0 引言

時間域航空電磁法(Airborne Time-domain Electromagnetic Method,ATEM)利用飛機上搭載的發射線圈向地下發送一次脈沖磁場(一次場)，地下的地質體由一次場激發出感應渦流，產生隨時間變化的感應場(二次場)，同時使用飛機上搭載的接收線圈采集從地下返回的電磁信號，通過分析采集到的電磁信號，探測地下地質體[1]。ATEM具有良好的靈活性和機動性，勘查速度快、成本相對較低，受地形、地貌條件的影響較小等優勢，能夠應用于油氣[2]和礦產資源勘探[3]、水文地質調查[4]以及環境監測[5-6]等領域。由于ATEM數據量巨大，且二維[7]和三維反演[8-10]需要大量的計算時間，因此很少應用于實際工程。一維反演因其計算量相對較小，在實際的ATEM數據處理和解釋中有著廣泛的應用[4-5]。

ATEM數據是全時連續采集，數據密度大，測點間距小，在大部分地質環境中，尤其是沉積環境下，相鄰測點間的反演結果應具有較強的橫向連續性。傳統的一維反演方法，如Occam反演[5,11]、自適應正則化反演[12-13]、阻尼特征參數反演[14]等只考慮單個測點的情況，既使通過選取合理的初始值對反演進行優化[15]，其反演結果的橫向連續性也較差。為此，Auken等[16]在反演直流電法數據時利用相鄰測點間的橫向粗糙度對模型進行約束(Laterally constrained inversion,LCI)，獲得了較好的反演結果；此外，Auken等[17]還應用LCI對古河道模型的瞬變電磁數據進行反演，并與一維單點反演結果進行對比，結果表明LCI可以在一定程度上提高反演分辨率，同時還能抑制數據噪聲的影響。Vallee等[18]將橫向約束應用于ATEM數據，并詳細地討論了橫向約束中一階導數約束和二階導數約束對反演結果的影響。為了突出層參數橫向連續的重要性，蔡晶等[19]引入加權橫向約束的概念。然而，由于LCI同時對多個測點進行計算,隨著測點數的增多，模型參數的數量也隨之增加，導致了計算時間急劇增加。Siemon等[20]對頻率域航空電磁數據進行反演，并討論了LCI模型參數的數量與計算時間之間的關系：當模型參數數量超過1000時，反演時間會急劇增加。Yu等[21]用前一個測點的反演結果作為先驗信息約束當前測點，以此保證反演結果在橫向上的連續性。但是這種橫向先驗信息約束很大程度上依賴于前一個測點的反演結果，因此反演結果具有明顯的方向性。

ATEM勘探中，單條測線的測點數量從幾千到幾萬不等，數據量巨大，因此對整條測線上所有測點同時進行計算，會耗費大量的時間。為了減少計算時間和提高反演精度，本文對整條測線上的測點進行分段計算，對每一段同時使用橫向和垂向粗糙度約束以確保反演結果的連續性，段與段之間使用Yu等[21]提出的橫向先驗信息進行約束，確保反演結果在段間的連續性。

本文首先詳細闡述這種橫向分段組合約束反演的實現方案，然后使用模擬數據和實測數據對本文反演方案進行測試，最后分析上述三種約束對反演結果的影響，驗證了本文反演方案的實用性和有效性。

1 反演方法

1.1 數據擬合

ATEM數據反演通過不斷調整迭代模型或迭代模型參數m，使正演響應g(m)與觀測數據dobs實現最佳擬合。將觀測數據通過一階泰勒級數展開，建立模型參數與觀測數據之間的線性逼近關系[16,22]

dobs?g(m)+G(mtrue-m)

(1)

式中：mtrue表示真實模型數據；G是雅可比矩陣，其第u行、第v列元素為

(2)

式中：du表示第u個正演響應數據；mv表示第v個模型參數。令

δm=mtrue-m

(3)

δdobs=dobs-g(m)

(4)

則式(1)變為

Gδm-δdobs=eobs

(5)

式中eobs為截斷誤差。本文采用阻尼最小二乘法對式(5)進行求解。由于只考慮了數據擬合誤差函數，因此反演結果具有多解性。為此，除了擬合誤差函數，還需對模型進行約束，以降低反演的多解性。

1.2 橫向粗糙度約束

假設地質模型在橫向上的電性特征變化較平緩，為了保證反演結果在橫向上具有較好的連續性，需要施加橫向粗糙度約束[16]

Rlsmtrue=els

(6)

式中：els表示模型參數在橫向上的誤差，其值越小，代表模型在橫向方向上越光滑；Rls是模型的橫向約束矩陣。將式(6)的左、右兩邊同時減去Rlsm，得到

Rlsmtrue-Rlsm=-Rlsm+els

(7)

則有

Rlsδm+Rlsm=els

(8)

橫向約束矩陣Rls采用模型在橫向上的一階差分矩陣，即在橫向相鄰的兩模型參數間使用一階差分，其他位置均為0

Rls=

(9)

式中：S1=(M-1)N;T1=MN。其中M為測點個數，N為模型層數。

1.3 橫向先驗信息約束

橫向先驗信息約束是將前一個點的反演結果作為先驗信息約束當前點的反演[21]。先驗信息約束假設模型與先驗信息的差距盡可能小，即

mtrue=mprior+elp

(10)

式中：mprior為模型的先驗信息，本文指上一段最后一個測點的反演結果；elp是模型與先驗信息的誤差，elp越小，代表模型與先驗信息越接近。將式(10)左、右兩邊同時減去m，有

mtrue-m=mprior-m+elp

(11)

得到

δm=mprior-m+elp

(12)

1.4 垂向粗糙度約束

垂向粗糙度約束[23]能使模型在垂直方向上具有較好的連續性。其推導過程與橫向約束一樣，因此與式(8)類似，即

Rvsδm+Rvsm=evs

(13)

式中下標“VS”代表垂直方向，各變量的含義與式(8)相同。與式(9)類似，Rvs的取值為

(14)

式中S2=MN,T2=M(N-1)。

1.5 橫向分段組合約束反演

圖1 橫向分段組合約束反演示意圖

根據圖1中的橫向分段組合約束反演的計算策略，將數據擬合方程(式(5))、橫向粗糙度約束方程(式(8))、橫向先驗信息約束方程(式(12))以及垂向粗糙度約束方程(式(13))聯合，可以得到總的反演方程

(15)

令

(16a)

(16b)

則式(15)可以寫為

G′δm-δd′=e′

(17)

根據文獻[24],式(17)的解為

δmest=[(G′)TWG′]-1(G′)TWδd′

(18)

式中W為模型參數和數據項的權值矩陣，可以寫成

(19)

其中Wobs、Wls、Wvs和Wlp均為對角加權矩陣，Wobs為數據項的加權矩陣，Wls為模型項的橫向約束加權矩陣，Wvs為模型項的垂向約束加權矩陣，Wlp為先驗模型的加權矩陣。為了使式(18)解更穩定，引入阻尼因子，將其改為迭代形式

(20)

式中：λn為第n次迭代的阻尼因子，求解方法見文獻[25]；mn為第n次迭代的模型參數。將式(16)和式(19)代入式(20)，有

(21)

2 實驗驗證

2.1 理論數據反演

為了驗證該方法的有效性，本文使用2.5維正演程序[26]模擬ATEM響應數據。ATEM數據采集系統見圖2，發射和接收系統主要參數為：發射線圈面積為123.9m2，匝數為5；發射基頻為25Hz，發射波形為三角波，上升沿時長為1.5ms，下降沿時長為1.50ms，峰值電流為300A；接收和發射線圈離地高度均為30m。理論模型包括3層：第1層的厚度為40m，電阻率為100Ω·m；第2層電阻率為5Ω·m，兩邊厚度為40m，中間厚度為100m，在中間位置有一寬度大約為200m、厚度逐漸變化的區域；第3層是基底，電阻率為100Ω·m。測線長度為1600m，測點間距為25m，測點個數為65。在測量發射電流關斷后0.04～4.00ms時間段內按照對數等間距抽取21道的感應電動勢(圖3)。為了盡量接近實測數據，對模擬的電磁響應數據添加5%的高斯白噪聲。

圖2 航空電磁采集系統示意圖

圖3 理論模型的模擬電磁響應曲線 圖中曲線從上到下分別對應衰減早期到晚期的感應電動勢

反演初始模型是電阻率為100Ω·m的均勻半空間，共26層，每層厚度均為10m。圖4b～圖4g為不同方法的反演結果，其中圖4b～圖4d和圖4g采取分段反演方法，即測線上的所有測點分為7段，前6段的測點數均為10，最后一段的測點個數為5，而圖4e采用單測點依次反演，圖4f則是所有測點同時反演。表1為不同反演方法的計算時間及反演結果與理論模型的均方根誤差。本文所用的計算平臺CPU型號為E3-1505M(8核)，內存為16GB。

從圖4b～圖4d、圖4g可以看出：無約束的分段反演結果(圖4b)整體效果較差，邊界模糊，中間部位電阻率值波動較大。進行了垂向粗糙度約束(圖4c)的反演結果得到了一定的改善，中間部位基本連成一塊，但是橫向連續性依然不理想。圖4d在圖4c的基礎上增加了橫向粗糙度約束，可見反演結果在橫向上的連續性得到了改善，但是由于段與段之間沒有進行橫向約束，反演結果出現了分塊現象。圖4g是在圖4d的基礎上在段與段之間增加了橫向先驗信息約束，其結果與圖4a中的理論模型比較接近。對比圖4d與圖4g可以發現，橫向先驗信息約束可以消除由分段引起的反演結果的分塊現象，提高反演的精度。

圖4 模型數據反演結果

(a)理論模型; (b)無約束分段反演; (c)垂向粗糙度約束分段反演; (d)垂向和橫向粗糙度約束分段反演; (e)垂向粗糙度約束和橫向先驗信息約束單點反演; (f)垂向和橫向粗糙度約束反演; (g)橫向分段組合約束反演

圖中黑線為理論模型界面

結合表1可以發現，在ATEM數據反演中，增加適當的約束可在提高反演結果準確性的同時減少計算時間。

對比圖4e與圖4g可以看出： 圖4e整體連續性較好，但由于采用前一個測點的反演結果作為當前測點的先驗信息來保持模型在橫向上的連續性，這只是單向約束，導致反演結果具有明顯的方向性，因而與圖4a中理論模型相差較大； 圖4g的橫向分段組合約束反演是在段與段之間采用單向的橫向先驗信息約束，而在段內采用雙向約束的橫向粗糙度約束，因此，反演結果受單向橫向先驗信息約束的影響較小，沒有出現類似于圖4e中明顯的方向性。因此，圖4e雖計算時間較少(表1)，但是反演結果的異常體的邊界不準確，誤差較大； 圖4g的計算時間稍長，但是反演精度較高，能真實地反映地下介質的電性分布特征。

對比圖4f與圖4g可以看出，二者與理論模型(圖4a)較為接近，邊界輪廓比較清晰，整體的連續性也較好。圖4f的計算時間是圖4g的近80倍，且誤差也稍大。因此，橫向分段組合約束反演可以在確保反演精度的同時大幅度減少計算時間。

表1 理論模型不同反演方法計算時間及均方根誤差(RMSE)

2.2 實測數據反演

為了驗證該橫向分段組合約束反演方案的實用性，對美國地質調查局(USGS)網站上公開的野外ATEM勘探數據(圖5)[27]進行反演。

每10個測點分為一段，反演結果見圖6。對比圖6a～圖6d可以發現，無約束的反演結果(圖6a)整體連續性較差，僅能夠反應地質體的大致位置，邊界較模糊；增加了垂向約束的反演結果(圖6b)在垂直方向上的連續性增強，但是地質體的形狀和輪廓仍然較為模糊，且第三層的反演效果不佳；增加了橫向粗糙度約束的反演結果(圖6c)的整體連續性增加，但是分段導致相鄰段之間的反演結果差別較大，不能很好地突出地質體的邊界和輪廓；圖6d整體較光滑，但是地質體的位置和形狀不能很好地得到反映；橫向分段組合約束的反演結果(圖6e)所反映的地質體位置、形狀和邊界較清晰。

圖5 實測電磁響應曲線

圖6 實測數據不同方式的反演結果

(a)無約束的分段反演結果; (b)垂向粗糙度約束的分段反演結果; (c)垂向和橫向粗糙度約束的分段反演結果; (d)垂向粗糙度約束和橫向先驗信息約束的單點反演結果; (e)橫向分段組合約束反演的反演結果

對比文獻[26]提供的地質資料，對圖6分析如下：

(1)圖中①處存在著一條山脈，分布著大量的基性火成巖，其電阻率大于1000Ω·m，與反演結果基本一致；

(2)圖6e中①、②之間的A處存在斷層,而圖6e中的電阻率在①、②之間的斷層A處的左右對比明顯，能夠很好地展現出斷層的垂向變化；

(3)圖中的②、③、④位于盆地中，其中②處的表層是由松散砂礫石組成的沖積扇沉積物，電阻率為80～150Ω·m,反演結果能夠很好地體現，③位于露頭花崗巖附近，圖中這個位置也呈現較高電阻率，④處富含粘土，具有較低的電阻率；

(4)圖中的⑤處位于盆地邊緣，主要由長英質火山巖組成，其電性與花崗巖相近，即具有較高的電阻率，與反演剖面結果一致；

(5)圖中B處存在一個橫切的斷裂帶，這可能是導致⑤處內部電阻率等值線不連續的原因；

(6)圖中⑥處表層黃色位置為長英質火山巖，具有較高的電阻率。藍色位置由較新的中等沖積扇沉積物組成，主要包含分選不良、致密的膠結砂和礫石，具有較低的電阻率(<10Ω·m)。

由此可見，實測ATEM數據反演所得到的構造、斷層及地質體的位置、形狀和邊界較清晰，反演結果與實際的地質資料吻合較好，說明本文提出的反演策略可有效地應用于實測ATEM數據的反演。

3 結論

本文將垂向粗糙度約束、橫向粗糙度約束以及橫向先驗信息約束三者相結合，對ATEM數據進行分段反演。模擬數據和實際數據的反演結果說明：

(1)適當的約束條件可提高ATEM數據反演結果的準確性，同時也能減少計算時間；

(2)橫向先驗信息約束可消除由分段引起的段間不連續現象；

(3)橫向粗糙度約束能有效地降低由橫向先驗信息約束引起的方向性。