瞬變電磁法非線性優(yōu)化反演算法對比

徐正玉，付能翼，周 潔，付志紅

1.重慶大學(xué)電氣工程學(xué)院，重慶 400044 2.科羅拉多礦業(yè)學(xué)院地球物理系，美國 科羅拉多 80401 3.重慶璀陸探測技術(shù)有限公司，重慶 402660

0 引言

瞬變電磁法(transient electromagnetic method，TEM)是利用不接地回線或接地電極向地下發(fā)送脈沖電磁場，用接收線圈或電極觀測由該脈沖電磁場感應(yīng)產(chǎn)生的二次場隨時間和空間的變化規(guī)律來解決有關(guān)地質(zhì)問題的時間域電磁法。目前，該方法已成為淺地表工程地球物理調(diào)查的重要方法之一，被廣泛應(yīng)用于工程地質(zhì)調(diào)查、水文地質(zhì)調(diào)查以及環(huán)境調(diào)查等諸多領(lǐng)域[1-7]。

通過數(shù)學(xué)物理方法將瞬變電磁二次場信號轉(zhuǎn)換為地層電阻率和深度信息是一個復(fù)雜的過程[8-12]。早期地球物理研究人員采用“煙圈”快速成像方法[13-14]。一些線性迭代方法也被用于地球物理數(shù)據(jù)反演計算中，如高斯-牛頓法、阻尼最小二乘法和自適應(yīng)正則化反演方法等等[15-19]。上述方法的優(yōu)點(diǎn)是收斂速度快，計算量較小。與此同時也存在以下不足。例如：“煙圈”快速成像方法近似粗糙，數(shù)據(jù)處理精度低；在反演過程中，線性迭代方法存在對初始模型依賴性大、反演結(jié)果容易陷入局部極小值、靈敏度矩陣計算量大以及對計算機(jī)處理器要求高等缺點(diǎn)。

因此，研究人員選擇將非線性優(yōu)化反演算法應(yīng)用到地球物理反演與成像中并獲得成功。1995年，Eberhart等[20]首次提出粒子群優(yōu)化(PSO)算法，該算法在全局搜索能力和收斂速度方面表現(xiàn)出良好的性能，已被成功用于解決許多領(lǐng)域的問題。Godio、Yuan和Pallero等[21-23]首先將PSO算法應(yīng)用于電測深數(shù)據(jù)、大地電磁數(shù)據(jù)以及其他地球物理方法數(shù)據(jù)處理。除此之外，其他非線性優(yōu)化算法也被用于地球物理數(shù)據(jù)反演。例如：Boschetti等[24]使用遺傳算法進(jìn)行位場數(shù)據(jù)處理反演；Roy等[25]將蒙特卡羅算法應(yīng)用于引力場數(shù)據(jù)和磁場數(shù)據(jù)反演；Yin等[26]將模擬退火算法應(yīng)用于航空電磁數(shù)據(jù)反演。在瞬變電磁數(shù)據(jù)反演中，Li等[27-28]應(yīng)用神經(jīng)網(wǎng)絡(luò)方法、差分進(jìn)化算法和PSO算法進(jìn)行近似一維反演；Xu等[29]將PSO算法應(yīng)用于TEM數(shù)據(jù)反演，研究結(jié)果表明該算法具有高效的優(yōu)點(diǎn)。但上述研究內(nèi)容未考慮隨機(jī)噪聲對PSO算法精度的影響。

2008年，英國劍橋?qū)W者Yang[30]提出一種新的智能算法，即螢火蟲算法。該算法是一種功能強(qiáng)大的跨學(xué)科智能隨機(jī)算法，已成功應(yīng)用于許多領(lǐng)域。算法的核心思想來源于螢火蟲利用發(fā)光的生物學(xué)特性進(jìn)行覓食和求偶的行為。螢火蟲算法被提出后，受到眾多研究者的關(guān)注和研究，并不斷得到完善。Abdullah[31]提出一種混合進(jìn)化螢火蟲算法；郝曉瑩[32]提出一種自適應(yīng)步長的螢火蟲算法；陳文平[33]提出一種多策略分層學(xué)習(xí)的螢火蟲算法。研究表明螢火蟲算法具有良好的搜索能力和優(yōu)化能力，易于實(shí)現(xiàn)，實(shí)用性強(qiáng)，參數(shù)少，過程簡單，被廣泛應(yīng)用于工程和金融等領(lǐng)域。而在地球物理反演與成像中，螢火蟲算法的發(fā)展非常緩慢。當(dāng)前，只有Zhou等[34]將螢火蟲算法應(yīng)用到瑞利波數(shù)據(jù)處理中；王鵬飛[35]將該算法應(yīng)用于大地電磁數(shù)據(jù)處理中。

綜上所述，PSO算法提出時間早，在地球物理反演中已經(jīng)發(fā)展成熟；而螢火蟲算法是近些年來提出的一種新方法，在地球物理資料反演中應(yīng)用較少。因此，為了比較這兩種非線性優(yōu)化算法的應(yīng)用效果，本文將PSO算法和螢火蟲算法應(yīng)用于TEM數(shù)據(jù)反演中。先詳細(xì)闡述PSO算法和螢火蟲算法的原理和求解流程，然后分析了隨機(jī)噪聲對兩種算法精度的影響并結(jié)合野外瞬變電磁實(shí)測數(shù)據(jù)進(jìn)行處理，以對比說明兩種算法的實(shí)際應(yīng)用效果。

1 非線性反演方法

1.1 粒子群優(yōu)化算法

PSO算法的核心思想是：將鳥或魚簡化為粒子，粒子的位置代表最優(yōu)化問題中的可能解，食物的位置代表最優(yōu)解，所有粒子在一定的規(guī)則下，向著最優(yōu)解位置運(yùn)動[20]。

在尋優(yōu)開始時，粒子群優(yōu)化算法首先在搜索區(qū)域范圍內(nèi)隨機(jī)初始化m個粒子，作為迭代初始值，然后實(shí)時更新自己的速度及位置：

(1)

(2)

PSO算法的求解過程(圖1)如下：

1)PSO算法參數(shù)設(shè)置。設(shè)定粒子數(shù)目和最大迭代次數(shù)，確定c1、c2、ω的取值。隨機(jī)產(chǎn)生每個粒子的初始位置、初始速度以及速度迭代變化的取值范圍。

圖1 PSO算法計算流程圖

2)根據(jù)目標(biāo)函數(shù)計算每個粒子的適應(yīng)度，并根據(jù)適應(yīng)度更新每個粒子的個體最優(yōu)位置和群體最優(yōu)位置。算法迭代計算第一次時，個體最優(yōu)位置就是初始位置，群體最優(yōu)位置是適應(yīng)度最好的粒子位置，當(dāng)最優(yōu)位置適應(yīng)度優(yōu)于迭代前時，更新最優(yōu)位置，否則保留原位置。

3)判斷誤差精度或達(dá)到最大迭代次數(shù)。判斷計算結(jié)果是否達(dá)到誤差精度要求或達(dá)到最大迭代次數(shù)。如否，進(jìn)行下次迭代，更新粒子速度和位置，并計算新的適應(yīng)度，再進(jìn)行步驟2)，直到滿足終止條件停止迭代，輸出結(jié)果并保存。

1.2 螢火蟲算法

螢火蟲算法是模擬自然界中螢火蟲發(fā)光的生物學(xué)特征發(fā)展而來的，是基于群體搜索的隨機(jī)優(yōu)化算法[8]。標(biāo)準(zhǔn)螢火蟲算法的核心思想是，每個螢火蟲的位置代表待求解問題的一個解，螢火蟲的發(fā)光亮度代表待求解問題的目標(biāo)函數(shù)值，目標(biāo)函數(shù)值越好，亮度越強(qiáng)。亮度弱的螢火蟲被亮度強(qiáng)的螢火蟲吸引，并向其移動。隨著迭代的進(jìn)行，種群中亮度弱的螢火蟲不斷向比自己更亮的螢火蟲移動，最終大多數(shù)螢火蟲會聚集在最亮的螢火蟲附近，最亮的螢火蟲的位置就是問題的最優(yōu)解。

如上所述，螢火蟲算法有兩個因素，即亮度和吸引度。定義螢火蟲的相對亮度為

I=I0·e-γrij。

(3)

式中：I0為螢火蟲絕對亮度，與目標(biāo)函數(shù)相關(guān)，目標(biāo)函數(shù)越優(yōu)絕對亮度越大；γ為光吸收系數(shù)；rij為螢火蟲i和j之間的距離。rij計算公式為

(4)

式中：xi和xj分別為螢火蟲i和j在空間中的位置；d為模型參數(shù)的數(shù)目。螢火蟲的吸引度為

(5)

式中：βij為螢火蟲i和j之間的吸引度；β0為最大吸引力，通常取1。當(dāng)微弱的螢火蟲移動到最亮的螢火蟲時，其位置更新為

(6)

式中：α為步長因子；R為[0, 1]內(nèi)的隨機(jī)數(shù)；ε為在參數(shù)范圍內(nèi)隨機(jī)生成的位置矢量。

螢火蟲算法的求解過程[35](圖2)如下：

1)螢火蟲算法參數(shù)設(shè)置。設(shè)置螢火蟲數(shù)目m、β0、γ、α，以及最大迭代次數(shù)或搜索精度。

2)隨機(jī)初始化螢火蟲位置，計算目標(biāo)函數(shù)，并將其倒數(shù)作為各自絕對亮度。

3)計算群體中螢火蟲的相對亮度和吸引度，根據(jù)相對亮度決定螢火蟲的移動方向。

4)更新螢火蟲的空間位置，根據(jù)更新后的螢火蟲位置重新計算螢火蟲的亮度。

5)當(dāng)滿足誤差精度或達(dá)到最大迭代次數(shù)時轉(zhuǎn)到下一步；否則增加迭代次數(shù)，轉(zhuǎn)到步驟3)進(jìn)行下一次迭代搜索。

6)輸出全局極值和最優(yōu)個體解，保存結(jié)果。

圖2 螢火蟲算法計算流程圖

2 適應(yīng)度目標(biāo)函數(shù)構(gòu)造

通過漢克爾變換和正、余弦變換或G-S逆拉普拉斯變換可以求解水平層狀介質(zhì)磁場強(qiáng)度垂直分量正演響應(yīng)[29]。構(gòu)造適應(yīng)度目標(biāo)函數(shù)為

(7)

式中：N為時間道數(shù)目；hcn和hsn為第n個時間道垂直磁場正演計算響應(yīng)和實(shí)測數(shù)據(jù)響應(yīng)。定義誤差為目標(biāo)函數(shù)，求解最優(yōu)解就是求目標(biāo)函數(shù)f為最小時的電阻率值。目標(biāo)函數(shù)越小即誤差越小，表明正演計算響應(yīng)與實(shí)測數(shù)據(jù)響應(yīng)越接近。

3 理論數(shù)據(jù)處理

3.1 無噪聲理論數(shù)據(jù)反演

首先，選擇無噪聲理論數(shù)據(jù)測試PSO算法和螢火蟲算法的性能。建立三層和四層理論地電模型。三層地電模型第一、二、三層電阻率(ρ1、ρ2、ρ3)分別為20、100、20 Ω·m，厚度(h1、h2、h3)分別為50、100 m和無限大；四層地電模型第一、二、三、四層電阻率(ρ1、ρ2、ρ3、ρ4)分別為100、20、100、50 Ω·m，厚度(h1、h2、h3、h4)分別為50、100、100 m和無限大。利用余弦變換和漢克爾變換進(jìn)行三層和四層理論地電模型的正演計算[29]，然后采用PSO算法和螢火蟲算法對其進(jìn)行處理。此外，還選擇“煙圈”快速成像方法進(jìn)行比較分析。在PSO算法中，c1=c2=2，ω=0.99k·r3/2+0.1(r3為[0, 1]內(nèi)的隨機(jī)數(shù))，m=50；在螢火蟲算法中，α=0.1，γ=0.8，m=50。兩種算法均迭代計算60次。

圖3為三層地電模型的反演結(jié)果。從圖3a可以看出，第二層高阻體的“煙圈”快速成像方法反演結(jié)果與理論模型相比誤差較大；而PSO算法和螢火蟲算法的反演結(jié)果與理論模型一致。從圖3b可以看出，在前5次迭代中，PSO算法的誤差遠(yuǎn)大于螢火蟲算法的誤差；但隨著迭代計算進(jìn)行，兩種算法的誤差均收斂于穩(wěn)定值且螢火蟲算法誤差一直小于PSO算法。

圖4為四層地電模型的反演結(jié)果。從圖4a中可以看出，“煙圈”快速成像方法反演結(jié)果與理論模型相比具有較大誤差，并且只有第一層和第二層的電阻率信息得到客觀反映，第三層和第四層效果不明顯，深度信息誤差大且精度低；而PSO算法和螢火蟲算法能較好地反映出理論模型的電阻率和深度信息，但是在電阻率較高的第三層中，兩種算法獲得的電阻率信息存在一定的誤差，且PSO算法誤差大于螢火蟲算法，而在電阻率較低的第二層和第四層中，兩種算法誤差較小。從圖4b可以看出，在前6次迭代中，螢火蟲算法的誤差大于PSO算法的誤差，但是螢火蟲算法的收斂速度比PSO算法快；通過迭代計算，螢火蟲算法的誤差逐漸低于PSO算法，并趨于穩(wěn)定值。

a. 反演曲線；b. 目標(biāo)函數(shù)隨迭代次數(shù)變化曲線。

a. 反演曲線；b. 目標(biāo)函數(shù)隨迭代次數(shù)變化曲線。

表1和表2給出了兩種非線性優(yōu)化算法的反演結(jié)果和相對誤差。從表1可以看出：在三層地電模型中，PSO算法對第一層電阻率的反演精度好于螢火蟲算法，但隨著層數(shù)增加，PSO算法誤差逐漸大于螢火蟲算法；而在四層地電模型中，除了第三層之外，PSO算法的反演精度均優(yōu)于螢火蟲算法。從表2可以看出，PSO算法反演深度的誤差整體上大于螢火蟲算法。因此，可以說明PSO算法對地層深度信息的反演精度較差。

在算法的計算時間方面，當(dāng)?shù)螖?shù)為60次時，PSO算法處理三層地電模型需要6.50 min，螢火蟲算法需要48.00 min。在處理四層地電模型時，PSO算法耗時7.60 min，而螢火蟲算法耗時約113.20 min。與上述兩種算法的計算時間相比，“煙圈”快速成像方法的計算時間幾乎可以忽略不計(表3)。

3.2 含噪聲理論數(shù)據(jù)反演

在實(shí)際工作中，采集到的信號可能包含各種噪聲信號，會導(dǎo)致處理結(jié)果不穩(wěn)定。因此，有必要分析算法的抗噪能力。將5%和10%的高斯隨機(jī)噪聲(圖5a)添加到四層理論地電模型中，并使用PSO算法和螢火蟲算法進(jìn)行處理。圖5b、c為含噪四層地電模型的反演結(jié)果。從圖5b可以看出，當(dāng)添加隨機(jī)噪聲時，PSO算法和螢火蟲算法都可以反映理論模型的電阻率和深度信息，但結(jié)果逐漸偏離理論模型參數(shù)。從圖5c可以看出，無論在哪種噪聲下，螢火蟲算法的目標(biāo)函數(shù)總是可以快速收斂到一個較小的穩(wěn)定常數(shù)，而PSO算法的目標(biāo)函數(shù)值在迭代計算中遠(yuǎn)遠(yuǎn)大于螢火蟲算法，且噪聲越大，目標(biāo)函數(shù)值越大。因此，可以分析得出螢火蟲算法的收斂性優(yōu)于PSO算法。

表1 兩種非線性優(yōu)化算法電阻率反演結(jié)果

表2 兩種非線性優(yōu)化算法深度反演結(jié)果

表3 算法計算時間對比表

a. 含噪磁場強(qiáng)度垂直分量；b. 含噪反演曲線；c. 目標(biāo)函數(shù)隨迭代次數(shù)變化曲線。

表4和表5為兩種非線性優(yōu)化算法處理含噪數(shù)據(jù)的結(jié)果和相對誤差，可以看出，無論在哪種隨機(jī)噪聲下，螢火蟲算法的誤差整體上小于PSO算法，這表明螢火蟲算法的抗噪能力優(yōu)于PSO算法。

4 實(shí)測數(shù)據(jù)處理

將PSO算法和螢火蟲算法應(yīng)用于某城市地鐵線路巖溶勘察。該項目調(diào)查的目的是尋找隧道下方和周圍潛在的巖溶分布。 調(diào)查區(qū)淺表層為第四系覆蓋層，地層深部巖性為完整性較好的灰?guī)r。TEM測線約40 m，各測點(diǎn)間距2 m。

圖6為瞬變電磁多測道圖。從圖6可以看出，瞬態(tài)電磁信號在1.498 4～6.157 6 ms之間呈現(xiàn)增加趨勢，說明是低阻異常體的反映。圖7為“煙圈”快速成像方法、PSO算法、螢火蟲算法的反演結(jié)果。PSO算法和螢火蟲算法處理每個測點(diǎn)的數(shù)據(jù)分別需要5.00和8.50 min。從圖7可以看出，與PSO算法和“煙圈”快速成像方法的處理結(jié)果相比，螢火蟲算法反演得到的低阻異常范圍更小，異常細(xì)節(jié)更加突出，具有較好的分辨率。

圖6 瞬變電磁多測道圖

表4 含噪四層地電模型電阻率反演結(jié)果

表5 含噪四層地電模型深度反演結(jié)果

a. “煙圈”快速成像方法；b. PSO算法；c. 螢火蟲算法。

5 結(jié)論

本文采用PSO算法和螢火蟲算法進(jìn)行瞬變電磁數(shù)據(jù)反演，取得以下結(jié)論。

1)建立典型的三層和四層地電模型，分析兩種算法的收斂特性和抗噪特性。研究結(jié)果表明：PSO算法計算效率高，但數(shù)據(jù)反演精度低，抗噪能力差；螢火蟲算法在抗噪性、數(shù)據(jù)反演精度等方面均優(yōu)于PSO算法。

2)利用PSO算法和螢火蟲算法對巖溶調(diào)查實(shí)測數(shù)據(jù)進(jìn)行反演處理。結(jié)果表明：螢火蟲算法反演得到的低阻異常體細(xì)節(jié)較突出，分辨率較高。研究內(nèi)容為TEM數(shù)據(jù)反演提供了新的方法和手段。