帶約束的多輻射場(chǎng)源半航空瞬變電磁一維自適應(yīng)正則化反演方法

張 瑩 瑩

(新疆大學(xué) 地質(zhì)與礦業(yè)工程學(xué)院，新疆 烏魯木齊 830047)

0 引言

半航空瞬變電磁法(semi-airborne transient electromagnetic method, semi-airborne TEM)，也稱為地空瞬變電磁法，采用地面發(fā)射、空中接收的工作方式，兼具地面發(fā)射功率強(qiáng)、勘探深度大和空中接收面積測(cè)量、工作效率高的雙重優(yōu)點(diǎn)，目前已在地?zé)嵴{(diào)查、火山結(jié)構(gòu)調(diào)查、地下巷道調(diào)查、地下水鹽漬化及地下水監(jiān)測(cè)、采空區(qū)探測(cè)、隧道勘察、古河道結(jié)構(gòu)探測(cè)等領(lǐng)域得到成功應(yīng)用[1-11]。半航空TEM的快速發(fā)展集中在近十年，目前對(duì)半航空TEM的研究仍處在初步階段，數(shù)據(jù)處理和解釋大都仍基于一維解釋方法[12]，研究?jī)?nèi)容主要集中在視電阻率求解、快速成像方面，如：陽貴紅分析了接收高度對(duì)電磁響應(yīng)的影響并設(shè)計(jì)了求解全區(qū)視電阻率的算法[13]；嵇艷鞠等基于三層BP神經(jīng)網(wǎng)絡(luò)和Levenberg-Marquardt算法實(shí)現(xiàn)了電性源半航空TEM神經(jīng)網(wǎng)絡(luò)視電阻率反演[14]；張瑩瑩等對(duì)多分量磁場(chǎng)響應(yīng)和磁場(chǎng)的時(shí)間導(dǎo)數(shù)響應(yīng)進(jìn)行分析并基于反函數(shù)原理實(shí)現(xiàn)了多分量的全域視電阻率定義[15]；吳啟龍研究了半航空TEM視電阻率成像，并將該算法應(yīng)用于復(fù)雜地形地區(qū)隧道勘察[10]；李貅、張瑩瑩等基于瞬變電磁擬地震成像技術(shù)提出了地空瞬變電磁法逆合成孔徑成像，并基于浮動(dòng)薄板原理研究了地空TEM快速成像算法[16-17]；呂仁斌利用查詢表方式求解視電阻率，并結(jié)合視深度公式研究了電導(dǎo)率深度成像[18]；王仕興等基于分段二分搜索算法建立起電導(dǎo)—電磁響應(yīng)數(shù)據(jù)庫(kù)，實(shí)現(xiàn)半航空TEM電導(dǎo)率深度快速成像[19]。這類研究的特點(diǎn)是計(jì)算速度快，不依賴初始模型，但其解釋結(jié)果大都很難提供與層厚相關(guān)的信息，且精度普遍較低，因此，諸如視電阻率計(jì)算和快速成像方法多用作定性分析。在一維反演方面，已有部分國(guó)內(nèi)外學(xué)者做了相關(guān)工作，如：張澎等引入并行技術(shù)實(shí)現(xiàn)了最平緩模型約束條件下的半航空時(shí)間域電磁數(shù)據(jù)自適應(yīng)正則化反演算法[20]；Abdallah等將頻率域電磁反演方法用于GREATEM的時(shí)間域數(shù)據(jù)，該法適用于對(duì)電阻率橫向變化大的數(shù)據(jù)進(jìn)行反演[21]；趙涵等以O(shè)ccam反演理論為基礎(chǔ)研究了多輻射場(chǎng)源地空瞬變電磁一維反演算法[22]；楊聰?shù)冉Y(jié)合自適應(yīng)正則化反演方法和阻尼最小二乘反演方法提出自適應(yīng)正則化—阻尼最小二乘反演算法，該法可在一定程度上提高半航空TEM對(duì)高阻的分辨率，提高反演精度[23]。目前已經(jīng)發(fā)表的一維反演方面的研究成果大都基于單個(gè)輻射源，實(shí)際上半航空TEM的輻射源可以是單一的，也可以是多個(gè)，磁性源也可視作多個(gè)輻射場(chǎng)源的疊加。多輻射場(chǎng)源能夠減少單個(gè)輻射源體積效應(yīng)的影響，在地下多方位激發(fā)地質(zhì)體，獲得多方位耦合信息,可通過調(diào)整電性源的位置或電流方向，有針對(duì)性地加強(qiáng)某一分量的采集信號(hào)強(qiáng)度，削弱隨機(jī)噪聲，提高半航空TEM對(duì)異常體的分辨能力[15,24]。在地形復(fù)雜地區(qū)，半航空TEM具有重要優(yōu)勢(shì)，若再結(jié)合多輻射場(chǎng)源技術(shù)，對(duì)于實(shí)現(xiàn)高分辨率深部探測(cè)具有重要意義。與航空電磁法類似，半航空TEM同樣面臨數(shù)據(jù)量大、數(shù)據(jù)信噪比低等問題，可以借鑒航空電磁資料反演中綜合各種先驗(yàn)約束信息并輔以適當(dāng)?shù)淖顑?yōu)化的正則化反演方法，降低反演多解性，實(shí)現(xiàn)半航空TEM數(shù)據(jù)穩(wěn)定快速的一維反演[25]。

本文結(jié)合半航空TEM的研究現(xiàn)狀和具體特點(diǎn)，以多輻射場(chǎng)源半航空TEM為例，實(shí)現(xiàn)了垂直分量帶約束的自適應(yīng)Occam反演算法，在Occam反演的基礎(chǔ)上[26]，采用陳小斌等提出的CMD自適應(yīng)調(diào)節(jié)方案改進(jìn)了拉格朗日乘子的處理方式[27]，省去傳統(tǒng)Occam反演中搜索拉格朗日乘子的大量計(jì)算過程，并利用自然邊界條件和模型修正量可行下降方向法對(duì)反演過程進(jìn)行約束，減少反演過程的計(jì)算量，保證反演過程的穩(wěn)定性和計(jì)算效率。典型層狀地電模型反演結(jié)果表明，該法可以通過較少的迭代次數(shù)獲得可靠的、抗干擾能力強(qiáng)的反演結(jié)果，為實(shí)現(xiàn)多輻射場(chǎng)源半航空TEM數(shù)據(jù)穩(wěn)定快速一維反演打下理論基礎(chǔ)。

1 多輻射場(chǎng)源TEM一維正演方法

圖1 地表電偶極子層狀模型示意Fig.1 Schematic diagram of an electric dipole on the surface of layered earth

(1)

根據(jù)電矢量位與磁感應(yīng)強(qiáng)度的對(duì)應(yīng)關(guān)系：

(2)

可得層狀大地地表電偶極子在空中產(chǎn)生的頻率域垂直分量磁感應(yīng)強(qiáng)度表達(dá)式：

(3)

(4)

考慮到實(shí)際工作中使用的發(fā)射源可能長(zhǎng)達(dá)數(shù)千米，在近區(qū)和過渡區(qū)，場(chǎng)源不能被看作電偶極子，本文采用剖分疊加、有限求和的方式，將發(fā)射源剖分成多段電偶極子，然后將每一段電偶極子產(chǎn)生的頻率域磁場(chǎng)疊加在一起近似計(jì)算總場(chǎng)源的頻率域響應(yīng)。對(duì)于每一段剖分的電偶極子，z軸始終指向下，選定一個(gè)電性源所在的坐標(biāo)系作為標(biāo)準(zhǔn)坐標(biāo)系(如圖2中的發(fā)射源1)，可得層狀大地多輻射場(chǎng)源在空中產(chǎn)生的頻率域垂直分量磁感應(yīng)強(qiáng)度表達(dá)式[28]：

圖2 多輻射場(chǎng)源剖分示意Fig.2 Coordinate sketch of multi-source subdivision

(5)

式中：PEij=Iidsij，Ii表示第i個(gè)電性源的供電電流，dsij表示第i個(gè)電性源剖分的第j個(gè)電偶極子的長(zhǎng)度；φij和rij分別表示測(cè)點(diǎn)M在地表的投影P與第i個(gè)電性源剖分的第j個(gè)電偶極子的夾角和對(duì)應(yīng)的偏移距；n為電性源總個(gè)；mi為第i個(gè)電性源的剖分段數(shù)。

當(dāng)采用階躍電流作為發(fā)射波形時(shí)，時(shí)間域響應(yīng)可以由頻率域響應(yīng)通過正弦變換得到：

(6)

2 多輻射場(chǎng)源瞬變電磁一維反演方法

2.1 Occam反演方法

Occam反演是一種改進(jìn)的帶平滑約束的最小二乘反演方法，在反演時(shí)要求實(shí)測(cè)數(shù)據(jù)與模型數(shù)據(jù)的擬合誤差在一定條件下達(dá)到極小，同時(shí)也要求反演模型的粗糙度達(dá)到極小，反演的目標(biāo)函數(shù)為[26]：

(7)

R=‖?m‖2；

(8)

其中?為粗糙度矩陣，形式如下:

(9)

在初始模型m1附近作線性化，將非線性問題轉(zhuǎn)為線性，則：

F[m1+Δ]=F[m1]+J1Δ。

(10)

式中：Δ是模型修正量；m1是反演一次迭代的解；J1是初始模型的偏導(dǎo)數(shù)矩陣，即雅克比矩陣，其元素

Jij=?Fi[m]/?mj，

(11)

mk+1={μ?T?+(WJ[mk])TWJ[mk]}-1·

{WJ[mk]}TWdk，

(12)

其中：

dk=d-F[mk]+J[mk]mk。

(13)

反演擬合誤差為實(shí)測(cè)數(shù)據(jù)和響應(yīng)數(shù)據(jù)的全域視電阻率均方根相對(duì)誤差：

(14)

反演終止標(biāo)準(zhǔn)為上述均方根相對(duì)誤差小于給定值。

2.2 拉格朗日乘子調(diào)節(jié)方案

根據(jù)式(7)，拉格朗日乘子的大小決定了反演的主要擬合對(duì)象。在反演中，反演效率、反演矩陣的病態(tài)程度和反演模型收斂性都與拉格朗日乘子有關(guān)。經(jīng)典的求取最佳拉格朗日乘子的方法是一維線性搜索，但該法每次迭代需要多次正演和求解反演方程，導(dǎo)致計(jì)算量很大；且每次迭代可能存在多個(gè)解，選取不當(dāng)可能導(dǎo)致反演不收斂。本文采用陳小斌等提出的CMD自適應(yīng)調(diào)節(jié)方案[27]，該方案在滿足總目標(biāo)函數(shù)極小的同時(shí)，還能近似滿足數(shù)據(jù)目標(biāo)函數(shù)極小和模型約束目標(biāo)函數(shù)極小的條件：

(15)

式中Δdk是第k次實(shí)測(cè)數(shù)據(jù)與模型數(shù)據(jù)的殘差。該方案不需要額外的正演計(jì)算，調(diào)解方案簡(jiǎn)單，可以有效減少計(jì)算量。

2.3 可行下降方向法約束

反演方程式(12)求解的是一個(gè)數(shù)學(xué)問題，為了盡量減少多解性的影響，引入可行下降方向法[29]。但反演過程中每次迭代的結(jié)果并不總是滿足自然規(guī)律，可能是極大值、極小值或負(fù)數(shù)，這與實(shí)際地質(zhì)情況顯然不符。為此，設(shè)定m下和m上分別為mk的下界和上界，即可行域，對(duì)每一次迭代的結(jié)果mk進(jìn)行約束，判斷每一層反演的電阻率是否超出給定的下界m下和上界m上，文中取m下=0.1 Ω·m，m上=1 000 Ω·m。若超出可行域，即不滿足自然邊界條件，給出如下可行下降方向法：

1) 對(duì)于可行域內(nèi)的mk，利用式(13)求出mk+1，計(jì)算出模型修正量Δmk=mk+1-mk。

2) 令

(16)

其中ε為給定的小正數(shù)。

3) 令

(17)

這樣求得的mk+1必在可行域內(nèi)。反演要達(dá)到更高的精度，要求相應(yīng)足夠小的步長(zhǎng)才有可能收斂，否則，在接近極小點(diǎn)附近時(shí)，模型修正量大，可能使反演結(jié)果快速陷入局部極小，導(dǎo)致迭代在不能滿足精度要求后進(jìn)入發(fā)散軌跡。在此建議模型修正量的最大值在50 Ω·m以內(nèi)。

這首詩暗用林和靖梅妻鶴子的典故，但不直白，頗含蓄蘊(yùn)藉，耐人尋味。他將長(zhǎng)安二月未見梅的原因詼諧地歸結(jié)于梅花還在等她的主人到來。董玘的七律也有個(gè)別寫得形象而富有哲理，如《題畫》：

2.4 雅可比矩陣優(yōu)化

由式(12)可知反演迭代的主要工作量在于雅克比矩陣的計(jì)算，若使用傳統(tǒng)的差分方法計(jì)算雅克比矩陣的矩陣元，若有N個(gè)反演參數(shù)就必須進(jìn)行N次正演，相應(yīng)的反演時(shí)間基本上就是正演耗時(shí)的N+1倍。瞬變電磁測(cè)深Occam反演的地層參數(shù)包括幾何參數(shù)和電參數(shù)，一般同時(shí)待反演的參數(shù)能達(dá)到幾十個(gè)，這就導(dǎo)致雅克比矩陣的計(jì)算量很大，因此研究雅克比矩陣的快速算法對(duì)于瞬變電磁資料反演具有重要意義。在此，借鑒文獻(xiàn)[30]中的方法對(duì)靈敏度矩陣進(jìn)行更新，減少雅可比矩陣的計(jì)算量：

(18)

3 層狀模型算例

為檢驗(yàn)上述反演算法的可靠性，對(duì)層狀模型的一維正演數(shù)據(jù)進(jìn)行反演。多輻射場(chǎng)源采用與電流方向相反的平行源，裝置布置示意見圖3，電性源長(zhǎng)度均為200 m，間隔200 m，發(fā)射電流為階躍脈沖，電流大小10 A，電流方向如圖中箭頭所示，在off-time取了20個(gè)時(shí)間道電磁響應(yīng)作為反演數(shù)據(jù)，接收線圈接收面積1 m2，匝數(shù)1，飛行高度100 m，以坐標(biāo)原點(diǎn)O(0,0)為例進(jìn)行分析。反演過程中，設(shè)置最大迭代次數(shù)20次，反演地層30層，首層厚度為5 m，層厚按對(duì)數(shù)等間隔進(jìn)行離散，反演的初始模型均為電阻率100 Ω·m的均勻半空間。

圖3 平行源裝置布置Fig.3 Arrangement diagram of parallel source device

3.1 三層模型

選擇三層模型中典型的H和K型地電模型進(jìn)行反演試算，H型模型各層電阻率分別為100、10、100 Ω·m，層厚為200、100 m；K型模型各層電阻率分別為100、1 000、100 Ω·m，層厚為200、100 m。圖4所示是H型模型反演結(jié)果，圖中顯示隨著反演迭代次數(shù)增加，拉格朗日乘子和擬合誤差快速下降，反演算法快速穩(wěn)定收斂，反演結(jié)果逐漸逼近真實(shí)模型，最終反演得出的各地層的深度和電阻率均與真實(shí)模型吻合較好。通過10次反演迭代，擬合均方根相對(duì)誤差降至約1%，已能獲得較為理想的反演結(jié)果。

圖4 H型模型反演結(jié)果Fig.4 Inversion results of H model

圖5所示是K型模型反演結(jié)果。圖中顯示隨著反演迭代次數(shù)的增加，拉格朗日乘子、擬合誤差、反演結(jié)果均表現(xiàn)出與H型模型相似的變化趨勢(shì)，反演結(jié)果基本能反映出真實(shí)模型的特征。所不同的是，由于淺部低阻層的屏蔽作用，最終反演結(jié)果中間高阻層的電阻率值無法反演出真實(shí)模型的電阻率值，這也與瞬變電磁法對(duì)高阻層不夠靈敏的特點(diǎn)相吻合。

圖5 K型模型反演結(jié)果Fig.5 Inversion results of K model

3.2 四層模型

設(shè)計(jì)了更為復(fù)雜的四層HK模型對(duì)反演算法進(jìn)行驗(yàn)證。HK型模型各層電阻率分別為200、50、300、100 Ω·m，層厚為100、100、100 m。隨著反演迭代次數(shù)的增加，拉格朗日乘子、擬合誤差、反演結(jié)果表現(xiàn)出與H型和K型模型相似的變化趨勢(shì)，反演過程快速穩(wěn)定，最終能夠得到與真實(shí)模型特征吻合的反演結(jié)果。圖6顯示，本方法對(duì)第二層(低阻層)的反演效果很好，與真實(shí)模型很接近，但對(duì)高阻層沒有準(zhǔn)確反演出真電阻率值；隨著迭代次數(shù)增加，拉格朗日乘子逐漸減小，數(shù)據(jù)擬合誤差逐漸減小，最后趨近于0.4%，經(jīng)過10次迭代，本文算法已經(jīng)可以獲得較好的反演結(jié)果。

圖6 HK型模型反演結(jié)果Fig.6 Inversion results of HK model

4 含噪數(shù)據(jù)算例

在野外實(shí)際工作中，由于各種噪聲源摻雜，半航空瞬變電磁測(cè)量的數(shù)據(jù)包含各種噪聲，如隨機(jī)噪聲和相關(guān)噪聲。為了評(píng)價(jià)噪聲對(duì)上述反演算法的影響，按文獻(xiàn)[25]的方法對(duì)一維正演所得的瞬變電磁響應(yīng)添加高斯噪聲模擬含噪數(shù)據(jù)進(jìn)行一維反演，t時(shí)刻一維正演瞬變響應(yīng)數(shù)據(jù)中加入如下所示高斯隨機(jī)噪聲：

(19)

式中：cnoise為常數(shù)，表示噪聲大??；tNd和VNd分別表示最后一個(gè)(第Nd個(gè))延時(shí)及該時(shí)刻的理論瞬變響應(yīng)值。這種噪聲表現(xiàn)為晚期信號(hào)信噪比低的特征，與實(shí)際情況相符。

4.1 三層模型

圖7 cnoise=0.3高斯隨機(jī)噪聲H型模型反演結(jié)果 Fig.7 Inversion results of H model with cnoise=0.3 Gaussian noise

圖8所示為添加cnoise=0.3高斯噪聲的K型模型反演結(jié)果，由于電磁場(chǎng)對(duì)高阻的不敏感性，對(duì)中間高阻層的反演結(jié)果一般，不能準(zhǔn)確刻畫高阻層電阻率值。隨著迭代次數(shù)的增加，拉格朗日乘子逐漸降低，但受噪聲影響，擬合均方根相對(duì)誤差在迭代8次時(shí)就逐漸趨于穩(wěn)定。對(duì)比圖5a和圖8a可見，噪聲對(duì)該模型反演結(jié)果的影響同樣主要體現(xiàn)在深部，對(duì)應(yīng)模型深度大于500 m的反演結(jié)果。

圖8 cnoise=0.3高斯隨機(jī)噪聲K型模型反演結(jié)果Fig.8 Inversion results of K model with cnoise=0.3 Gaussian noise

4.2 四層模型

為了分析噪聲強(qiáng)弱對(duì)反演結(jié)果的影響，采用前文中的HK型四層模型，在一維正演瞬變響應(yīng)中分別加入cnoise=0.1、0.3的高斯隨機(jī)噪聲，利用本文提出的方法對(duì)地下電阻率進(jìn)行反演。圖9所示為含cnoise=0.1高斯噪聲的反演結(jié)果。對(duì)比圖6a和圖9a可以看出，由于噪聲信號(hào)強(qiáng)度較弱，含噪與不含噪數(shù)據(jù)的反演結(jié)果差別不大，能夠較好地刻畫出地下電阻率變化趨勢(shì)，對(duì)于淺部低阻層反演效果較好，層界面位置也較為清楚，但是對(duì)于較深的高阻層反演效果一般，沒有準(zhǔn)確反演出真電阻率值，僅反映出在相應(yīng)層位存在高阻層，深部地層電阻率反演結(jié)果與模型真電阻率吻合很好。

圖10所示為含cnoise=0.3高斯噪聲的反演結(jié)果，由圖可見，擬合誤差和拉格朗日乘子隨著迭代次數(shù)逐漸下降，在8次迭代后，擬合均方根相對(duì)誤差就逐漸趨于穩(wěn)定。與圖9a不同的是，由于噪聲信號(hào)強(qiáng)度較大，雖然能夠刻畫出地下電阻率變化趨勢(shì)，但深部地層電阻率反演結(jié)果與模型電真阻率差別較大。對(duì)比圖9b和圖10b可以看出，噪聲信號(hào)強(qiáng)度會(huì)影響擬合均方根相對(duì)誤差，噪聲信號(hào)越強(qiáng)，最終的擬合誤差越大。

圖9 cnoise=0.1高斯隨機(jī)噪聲HK型模型反演結(jié)果 Fig.9 Inversion results of HK model with cnoise=0.1 Gaussian noise

圖10 cnoise=0.3高斯隨機(jī)噪聲HK型模型反演結(jié)果Fig.10 Inversion results of HK model with cnoise=0.3 Gaussian noise

5 結(jié)論

本文基于垂直分量研究了帶約束的多輻射場(chǎng)源地空瞬變電磁一維自適應(yīng)正則化反演方法，并將其分別應(yīng)用于含噪和不含噪一維層狀模型數(shù)據(jù)的反演，得出如下結(jié)論。

1)該算法基于Occam反演，采用CMD自適應(yīng)拉格朗日乘子調(diào)節(jié)方案改進(jìn)拉格朗日乘子的處理方式，并利用自然邊界條件和模型修正量可行下降方向法對(duì)反演過程進(jìn)行約束，減少反演過程的計(jì)算量，提高反演效率，計(jì)算速度快，收斂穩(wěn)定，經(jīng)過10次左右的迭代，已能獲得較為理想的反演結(jié)果。

2)層狀模型反演結(jié)果表明該算法可以獲得光滑穩(wěn)定的反演結(jié)果，對(duì)于典型的H、K型三層模型和HK型四層模型均能獲得較好的反演結(jié)果，該算法對(duì)于復(fù)雜地電構(gòu)造仍具有較好的適應(yīng)性和有效性；受低阻層屏蔽作用和對(duì)高阻體不靈敏的影響，該算法對(duì)高阻層的反演僅能獲得下部地層的平均電阻率值。

3)噪聲信號(hào)大小會(huì)影響反演結(jié)果，噪聲越大，對(duì)反演結(jié)果的影響越大，由于本文添加的高斯隨機(jī)噪聲表現(xiàn)為晚期數(shù)據(jù)信噪比低，噪聲對(duì)反演結(jié)果的影響主要體現(xiàn)在深部地層電阻率反演結(jié)果上；噪聲同樣會(huì)影響擬合均方根相對(duì)誤差，噪聲越大，最終的擬合誤差越大。