基于監督下降法的大地電磁二維反演及應用研究

付興,譚捍東,董巖,汪茂

(中國地質大學(北京) 地球物理與信息技術學院,北京 100083)

0 引言

大地電磁測深法(MT)是一種通過觀測天然交變電磁場來對地電構造進行研究的物探方法,大地電磁法的頻域很寬,其探測深度距地表有幾百米到數百公里不等,具有分辨率高,穿透能力強,勘探深度大以及效率高等實際優點,被廣泛應用于金屬礦產普查、勘探地熱和油氣田等領域[1]。

大地電磁二維反演常采用高斯牛頓法[2]、OCCAM法[3]、非線性共軛梯度法[4]等,這些方法在最小化目標函數時各有優劣,但都無法在迭代過程中避免雅可比矩陣或海森矩陣的計算,且存在陷入局部極小值的風險,導致反演效果變差。

近些年隨著計算機行業的發展,機器學習成功在多種領域取得不錯的效果[5]。機器學習能夠挖掘變量與標簽之間的關聯,并學習其關聯進行預測,大大減少反演時間,提升計算效率[6]。

監督下降法(SDM)在人臉識別技術中被廣泛使用[7-8],后被帶入到多領域進行研究[9-11]。采取監督下降法解決地球物理勘探領域的問題也取得了不錯的成效[12],模型與正演響應合并構成的訓練集在進行一維瞬變電磁(TEM)反演中發揮了重要的作用[13]。監督下降法分為離線訓練和在線預測兩個階段,在離線訓練階段,通過學習大量模型及其響應得出平均下降步長。在預測階段,采用學習的步長進行迭代預測。平均步長具有較好的指導作用,避免了迭代陷入局部極小的風險,提升了反演的效果。此外,監督下降法并不需要計算雅可比矩陣和海森矩陣,避免了內存浪費的同時,提升了反演的效率。

本文采用監督下降法理論改進大地電磁二維數據反演算法,編寫相關代碼并設計理論模型驗證算法的正確性,對其抗噪能力及實用性進行檢驗。相較于傳統的非線性共軛梯度反演,基于監督下降法的反演具有收斂速度快、反演效果好、抗噪能力強等特點。

1 大地電磁二維正反演簡介

1.1 大地電磁二維有限單元法正演

忽略位移電流和磁導率的影響,時間因子為e-iωt,大地電磁場滿足麥克斯韋方程組:

(1)

式中:σ為介質的電導率;μ0為介質的磁導率;E為電場強度;H為磁感應強度。由式(1)可推導出赫姆霍茲方程組:

(2)

(3)

式中:AB為上邊界,CD為下邊界。

(4)

對圖1所有單元進行離散,式(3)變為:

(5)

圖1 部分單元網格Fig.1 Sketch of partial grid

式中:K為總體系數矩陣,求取式(5)的極值,可得

Ku=0,

(6)

可得出各節點的u,相應的視電阻率及阻抗相位為:

(7)

1.2 大地電磁二維非線性共軛梯度反演

常規的大地電磁二維反演方法有高斯牛頓法、OCCAM法、非線性共軛梯度法。本文將采用非線性共軛梯度反演與基于監督下降法的反演進行對比,下面概要介紹非線性共軛梯度反演方法。

定義大地電磁二維反演目標函數如下[15]:

Φ(m)=(dobs-F(m))TV-1(dobs-F(m))+
λ(m0-m)TLTL(m0-m),

(8)

式中:m為模型參數;F(m)為求取阻抗張量的正演函數;dobs為模型正演響應或實測數據;V為協方差矩陣;λ為正則化參數;L為二次差分拉普拉斯算子;m0為先驗信息,目標函數具有數據擬合差最小和模型最光滑的雙重約束。

目標函數的梯度為:

g=-2ATV-1e+2λLTL(m0-m),

(9)

其中:A表示雅可比矩陣,數據誤差向量:e=dobs-F(m)。

非線性共軛梯度反演的具體步驟如下[16-17]:

1)i=0時,設定初始模型mi,計算梯度ri=gi;

3)求解ri=0,得出最優搜索步長αi;

4)mi+1=mi+αiui,ri+1=gi+1;

6)取i=i+1,回到步驟3),直至收斂,反演結束。

2 基于監督下降法的大地電磁二維反演

監督下降法運用機器學習理論,采取半監督學習方式,分為訓練和預測兩個階段。訓練階段通過對訓練集的學習,得到一組平均下降方向。預測階段運用平均下降方向迭代計算數據殘差,進而進行迭代。監督下降法的流程如圖2。

2.1 訓練階段

在進行大地電磁反演計算時,定義目標函數如下:

(10)

對式(10)的正演函數F(m)進行一階Taylor展開,求取對m的極小值,可得:

Δm=(JTJ)-1JTΔd,

(11)

其中:J為雅可比矩陣;Δm=m-m1;Δd=dobs-F(m);α=(JTJ)-1JT,α為迭代步長。由于正演數據存在非線性,常規反演通常采取不斷迭代的方式來求取Φ(m)關于m的極小值,計算雅可比矩陣及海森矩陣是在進行常規反演時較為消耗計算時間和占用內存空間的一步。況且,上述目標函數在此時容易陷入局部極小,所以,需要大量的模型迭代來獲得迭代最優方向。

此時,相較于每次迭代均需計算雅可比矩陣或海森矩陣且每步均需搜索下降步長的常規大地電磁反演,引入機器學習的思想,通過學習的平均下降方向,使目標函數快速趨于極小值的方法可以節約大量計算時間以及大量內存空間,同時避免目標函數陷入局部極小值。

合成一組訓練集(標簽集MP和數據集DP),為了計算下降方向,求下式極小值:

(12)

(13)

(14)

αi=(ΔDiTΔDi+λI)-1ΔDiTΔMi,

(15)

式中:i為迭代次數;λ為阻尼系數(可與相加的內容呈現比例關系),Mi+1表達式如下:

Mi+1=Mi+ΔDiαi,

(16)

在此階段,數據擬合差RMSD與模型擬合差RMSM為:

(17)

當模型擬合差趨于穩定、足夠小或者達到了預設的最大迭代次數時,訓練階段結束。

2.2 預測階段

將上階段的平均步長用于該階段的計算,具體過程如下:

mi+1=mi+[dobs-F(mi)]αi,

(18)

從i=1開始計算,m1與訓練過程的初始模型相同。當F(mi)與dobs間的數據殘差為0時,停止迭代,輸出反演結果。

為解決反演的多解性問題,考慮在預測階段加入吉洪諾夫正則化來解決這一問題。

(19)

式(19)為進行正則化的目標函數,其中υv、υh為垂直方向以及水平方向的正則化系數。

(20)

式中:

(21)

求取式(19)的極小值進行偏導處理,可得下式,化簡后,可得

Gimi=di,

(22)

(23)

求解式(22),解出的mi,利用mi進行下次迭代,直至數據擬合差停止變小或小于一定數值,數據擬合差公式如下:

(24)

最大訓練次數的步長預測完后,得出的數據擬合差未達到標準或者預測結果未達到期望,采取重啟步長算法,即從第一個步長開始,進行新一輪的預測,直至數據擬合差足夠小或不再減少。

3 理論模型合成數據二維反演算例

大地電磁的二維正演采用有限單元法,將地下介質剖分成33×26的網格單元,橫向網格中心區域均勻剖分,兩側不斷變大?？v向為26個不斷增大的非均勻網格,設計40個頻點。

將部分不同電阻率的規則形狀的異常體在背景電阻率為100 Ω·m的均勻半空間內移動,形成了4 355個不同模型,如圖3所示。

圖3 部分訓練集Fig.3 Some training samples

將上述模型分別進行正演計算,使用正演結果中TM模式的數據信息,將上述模型與正演響應構成一個訓練集,進行離線訓練,本階段采用matlab及C++混合編程的方法,模型試算工作采用的CPU型號為i5-12500H,開啟10線程并行,訓練共耗時2 h。

訓練過程共迭代8次,初始迭代誤差為0.423,迭代8次后誤差縮小為0.087。迭代過程如圖4所示。從訓練集中隨機取出80組模型進行驗證工作,其中數據誤差低于15%的有80組(如圖5所示),可見離線訓練所求平均下降方向滿足反演要求。

圖5 誤差分布Fig.5 The histogram of the normalized model misfit

3.1 理論模型1

模型1(圖6)包含4個異常體,背景電阻率為100 Ω·m,異常體1是大小為800 m×800 m、埋深為400 m、電阻率值為1 000 Ω·m的高阻異常體,異常體2是大小為800 m×700 m、埋深為300 m、電阻率值為10 Ω·m的低阻異常體,異常體3是大小為600 m×700 m、埋深為300 m、電阻率值為1 000 Ω·m的高阻異常體,異常體4是大小為1 000 m×1 000 m、埋深為400 m、電阻率值為10 Ω·m的低阻異常體。

圖6 理論模型1示意(從左到右分別為異常體1、2、3、4)Fig.6 Schematic diagram of model 1 (abnormal bodies 1, 2, 3 and 4 are arranged from left to right)

圖7a為模型1正演響應的反演結果,圖7b為10%高斯誤差合成數據的反演結果,圖8為兩者的數據擬合曲線?？梢?本算法具有一定的抗噪能力。

a—高斯誤差為0%的預測結果 ;b—高斯誤差為10%的預測結果

圖8 迭代過程的數據擬合差Fig.8 The normalized data misfit in iteration

圖9a為非線性共軛梯度法的二維反演結果,用時5.695 s。高阻異常體1、3均未能恢復理論模型電阻率,能還原異常體大致位置;低阻異常體2、4未能恢復理論模型電阻率,大致還原異常體位置。

a—非線性共軛梯度反演結果;b—監督下降法反演結果

圖9b為監督下降算法反演的結果,用時3.145 s。高阻異常體1、3理論模型電阻率的恢復效果優于圖9a,且能夠還原異常體位置;低阻異常體2、4的反演結果基本可以恢復理論模型電阻率,也能夠還原異常體位置。

3.2 理論模型2

模型2(圖10)包含4個異常體,背景電阻率為100 Ω·m,異常體1是大小為1 400 m×1 000 m、埋深為1 200 m、電阻率值為1 000 Ω·m的高阻異常體,異常體2是埋深300 m、電阻率值為10 Ω·m的階梯狀低阻異常體,異常體3是大小為1 200 m×260 m、埋深為40 m、電阻率值為10 Ω·m的低阻異常體,異常體4是大小為600 m×1 300 m、埋深為850 m、電阻率值為1 000 Ω·m的高阻異常體。

圖10 理論模型2示意(從左到右分別為異常體1、2、3、4)Fig.10 Schematic diagram of model 2 (abnormal bodies 1, 2, 3 and 4 are arranged from left to right)

圖11a為非線性共軛梯度法的二維反演的結果,用時5.6 s。高阻異常體1未能恢復理論模型的電阻率,也未能還原異常體位置;非線性共軛梯度反演方法僅恢復了低阻異常體2地下1 000 m以上的信息,但基本無法得出1 000 m以下的有效信息,可能受低阻異常體3的影響較大;低阻異常體3大致恢復理論模型的電阻率,并還原異常體位置;高阻異

a—非線性共軛梯度反演結果;b—監督下降法反演結果

常體4未恢復理論模型的電阻率,也未還原異常體大致位置。

圖11b為監督下降算法反演的結果,初始迭代誤差為0.546,迭代10次后迭代誤差為0.063 1,用時2.8 s。高阻異常體1理論模型電阻率的恢復效果優于11a,且能還原異常體位置;監督下降算法可以反演出低阻階梯狀異常體2的地下2 000 m以上的異常體信息,并恢復異常體理論模型電阻率,2 000～3 000 m的部分沒有很好的體現,可能是低阻異常體3屏蔽的影響;低阻異常體3能恢復理論模型的電阻率,也能還原異常體位置;高阻異常體4理論模型電阻率恢復效果優于11a,可能是低阻異常體3的影響,也未還原異常體位置。

3.3 理論模型3

模型3(圖12)為復雜地層模型,背景電阻率為100 Ω·m,存在地層平均層厚500 m,頂部埋深100 m,電阻率值10 Ω·m的中部突起地層異常體。

圖13a為非線性共軛梯度法的二維反演的結果,用時6.2 s。圖13b為監督下降算法反演的結果,初始迭代誤差為0.842,迭代10次后迭代誤差為0.085,用時3.4 s。監督下降算法的形態更接近真實模型,對異常體的電阻率恢復效果更好。

a—非線性共軛梯度反演結果;b—監督下降法反演結果

總體上,監督下降法的反演效果優于非線性共軛梯度法的反演效果。訓練集中并沒有訓練較復雜的模型,但我們仍可以反演理論模型的高阻異常體、復雜異常體以及復雜地層異常模型,說明監督下降法具有處理較復雜模型的能力。

4 實測數據應用

為驗證算法的實用性,將本文開發的反演算法應用于西藏高原南部雅魯藏布江縫合帶地區實測數據的反演。測線南起錯那縣,終止于墨竹工卡縣,全長約240 km,測點29個,觀測頻率范圍為0.000 1～250 Hz。測線橫跨雅魯藏布江縫合帶,縫合帶附近存在大規模的岡底斯花崗巖[18]。

圖14為原始數據擬斷面圖,本文主要對TM模式的數據進行反演,對比監督下降算法及非線性共軛梯度算法反演結果,以此驗證監督下降法的實用性。

a—TM模式視電阻率擬斷面;b—TM模式相位擬斷面;c—TE模式視電阻率擬斷面;d—TE模式相位擬斷面

圖15a為非線性共軛梯度法的反演結果,圖15b為監督下降法的反演結果,兩者反演結果基本一致。圖16為反演結果與實測數據對比,可見大構造擬合程度較高。在時間方面,非線性共軛梯度大致需要30 s,監督下降法的訓練時間為2 h左右,預測時間大約25 s。

a—非線性共軛梯度法反演結果;b—監督下降法的反演結果

5 結論

本文實現了基于監督下降法的大地電磁二維反演,學習了一組訓練集,對三組理論模型和實測數據進行了反演,檢驗了反演算法的可行性和實用性。反演結果表明基于監督下降法的大地電磁二維反演收斂速度、反演效果優于非線性共軛梯度反演,具有較好的抗噪能力及實用性。

1)訓練得到的平均步長對異常體的預測起了指導性作用,避免了雅可比矩陣及海森矩陣的計算,大大降低了陷入局部極小值的風險,減少反演多解性。

2)就預測階段來說,需要的時間小于常規反演所需的時間,對于批量數據的處理具有一定優勢。

3)監督下降法具有較高泛化能力,可以用于更多地球物理勘探方法。

監督下降法總體耗時(訓練+預測)較常規反演略長,同時平均步長的儲存對硬盤的需求較大,并行降低耗時與減少儲存是其未來的研究方向。