基于量子行為粒子群算法的航空瞬變電磁擬二維反演技術(shù)

何一鳴，薛國(guó)強(qiáng),4*，趙 煬

(1.中國(guó)科學(xué)院地質(zhì)與地球物理研究所 中國(guó)科學(xué)院礦產(chǎn)資源研究重點(diǎn)實(shí)驗(yàn)室，北京 100029; 2.中國(guó)科學(xué)院地球科學(xué)研究院，北京 100029; 3.中國(guó)科學(xué)院大學(xué) 地球與行星科學(xué)學(xué)院，北京 100049; 4.長(zhǎng)安大學(xué) 地球物理場(chǎng)多參數(shù)綜合模擬實(shí)驗(yàn)室(中國(guó)地球物理學(xué)會(huì)重點(diǎn)實(shí)驗(yàn)室)，陜西 西安 710054)

0 引 言

航空瞬變電磁法發(fā)射和接收裝置皆位于空中，完全擺脫了地形的限制，具有高效率和適應(yīng)性強(qiáng)等特點(diǎn)，其廣泛應(yīng)用于礦產(chǎn)資源探測(cè)、地質(zhì)填圖和地?zé)豳Y源勘查等領(lǐng)域[1-7]。該方法在人工脈沖信號(hào)激發(fā)下產(chǎn)生感應(yīng)電磁場(chǎng)響應(yīng)，然后通過(guò)數(shù)據(jù)反演擬合重構(gòu)地電信息。為獲得更加接近地下真實(shí)模型的反演結(jié)果，研究精細(xì)反演算法顯得極為必要。反演算法可劃分為確定性反演算法和隨機(jī)性反演算法,目前主流反演算法為確定性反演算法，例如廣義逆矩陣[8]、共軛梯度[9]和擬牛頓[10]等確定性反演算法，具有計(jì)算速度快的特點(diǎn)，但是瞬變電磁法反問(wèn)題具有非線性和病態(tài)性，以及確定性反演算法反演結(jié)果依賴于初始模型的選擇等問(wèn)題,在求解上述反問(wèn)題時(shí)易受人為因素干擾陷入到局部極小值中，導(dǎo)致最終反演結(jié)果偏離真實(shí)參數(shù)模型[11]。隨后，改進(jìn)正則化方法被提出來(lái)，在確定性反演算法基礎(chǔ)上加入模型約束來(lái)減小反演結(jié)果的非唯一性，同時(shí)提高反演結(jié)果的穩(wěn)定性。其中Occam反演算法采用最大光滑穩(wěn)定泛函所得到光滑反演結(jié)果能夠保證電性變化上的連續(xù)性，但是對(duì)目標(biāo)體邊界刻畫(huà)不夠清晰[12]；而聚焦反演算法采用最小梯度支撐穩(wěn)定泛函得到尖銳反演結(jié)果能夠顯著提高目標(biāo)體邊界刻畫(huà)精度，但是仍然存在假異常干擾。

隨著近些年來(lái)計(jì)算科學(xué)的快速發(fā)展，隨機(jī)性反演算法逐漸成為研究熱點(diǎn)，主要包括粒子群(Particle Swarm Optimization,PSO)算法、模擬退火算法、遺傳算法、貝葉斯算法和神經(jīng)網(wǎng)絡(luò)算法等[6,13-16]。其中粒子群算法作為一種群優(yōu)化理論的隨機(jī)性反演算法，具有較好的跳出局部極小值的能力，自從被提出后廣泛應(yīng)用于基因工程[17]、復(fù)合材料[18]、數(shù)學(xué)建模[19]、生命科學(xué)[20]等領(lǐng)域中。2006年,Fernandez-Alvarez等率先將其引入到地球物理反演中，成功地實(shí)現(xiàn)了垂直電測(cè)深的反問(wèn)題求解[21]。最近十幾年，眾多學(xué)者先后開(kāi)展了基于粒子群算法的反演研究工作，極大推動(dòng)了粒子群算法在地球物理領(lǐng)域的發(fā)展，但是由于粒子群算法在求解非線性映射關(guān)系下的高維反問(wèn)題時(shí)，仍然存在早熟收斂和收斂速度慢等問(wèn)題，所以限制了粒子群算法在電磁法二、三維反演中的發(fā)展[15,22-24]。

本文首先分析了早熟收斂和收斂速度緩慢的原因，然后分別給出了相應(yīng)的解決方案。一方面，在傳統(tǒng)粒子群算法中，早熟收斂問(wèn)題主要由于粒子群分布集中或運(yùn)動(dòng)隨機(jī)性不足，使得各粒子快速向全局最優(yōu)粒子聚集。本文引入量子行為粒子群(Quantum-behaved Particle Swarm Optimization，QPSO)算法改善傳統(tǒng)的粒子群算法，其中引入量子在勢(shì)阱中運(yùn)動(dòng)規(guī)律指導(dǎo)粒子群在空間位置更新，使得粒子可以出現(xiàn)在勢(shì)阱范圍內(nèi)任何存在概率分布的位置上，進(jìn)而能夠克服群體聚集性導(dǎo)致的早熟收斂問(wèn)題[25-26]。另一方面，局部極小值個(gè)數(shù)直接影響反演收斂速度，考慮到粒子群算法中局部極小值個(gè)數(shù)與模型參數(shù)的空間維度成指數(shù)正相關(guān)關(guān)系，因此，本文提出采用擬二維反演算法代替?zhèn)鹘y(tǒng)二維反演算法，將二維模型拆分為一維，設(shè)法提高粒子群算法的收斂速度[27]。但是在傳統(tǒng)擬二維反演中，采用添加正則化項(xiàng)實(shí)現(xiàn)橫向約束，如果在群體中開(kāi)展基于L曲線法正則化參數(shù)尋優(yōu)過(guò)程將要耗費(fèi)大量計(jì)算資源[28]。結(jié)合量子行為粒子群算法中全局最優(yōu)粒子在粒子群進(jìn)化過(guò)程中的重要地位和積極影響，提出了采用非線性滑動(dòng)窗口濾波器α-Trimmed方法開(kāi)展全局最優(yōu)粒子的模型參數(shù)光滑約束技術(shù)，實(shí)現(xiàn)粒子群算法的快速擬二維反演[29]。最后將量子行為粒子群算法擬二維反演技術(shù)應(yīng)用到含噪全航空瞬變電磁數(shù)據(jù)處理中，反演結(jié)果更接近真實(shí)的地電模型參數(shù)，驗(yàn)證了算法有效性和魯棒性，為全航空瞬變電磁反演解釋提供一種合適的解釋方法。

1 航空瞬變電磁法反問(wèn)題

盡管航空瞬變電磁法不同裝置形式之間存在一定差異，但原理上正反演問(wèn)題都可以簡(jiǎn)要表述為

d=F(m)

(1)

式中：F是根據(jù)電磁傳播規(guī)律建立的數(shù)學(xué)物理方程(正演算子)；d為地表觀測(cè)的響應(yīng)數(shù)據(jù)；m為模型參數(shù)集。

在模型參數(shù)和正演算子已知的情況下，計(jì)算瞬變電磁響應(yīng)是一個(gè)正問(wèn)題。反之，如果觀測(cè)數(shù)據(jù)和正演算子已知，那么推斷地球內(nèi)部結(jié)構(gòu)或地下礦藏位置是一個(gè)典型的反問(wèn)題。

用于求解反問(wèn)題的反演算法可分為確定性反演算法[8-10]和隨機(jī)性反演算法[13-16]。確定性反演算法[圖1(a)]在模型空間中從某一個(gè)初始模型開(kāi)始，沿著靈敏度函數(shù)確定的方向移動(dòng)，當(dāng)正演響應(yīng)和記錄的數(shù)據(jù)集之間的擬合差滿足要求，即得到解，反演結(jié)果依賴于初始模型選取，易陷入局部最優(yōu)解。而隨機(jī)性反演算法[圖1(b)]在整個(gè)模型空間內(nèi)進(jìn)行隨機(jī)搜索，在模型空間中多個(gè)點(diǎn)同時(shí)開(kāi)始，依據(jù)個(gè)體經(jīng)驗(yàn)積累和群體的信息交流確定各個(gè)模型運(yùn)動(dòng)方向，且具有一定隨機(jī)性，最后將所有滿足擬合殘差要求的模型組成表征反問(wèn)題解的模型集，其中不僅包括一個(gè)全局最優(yōu)解，還會(huì)有多個(gè)次優(yōu)解，因此，隨機(jī)性反演算法具有更強(qiáng)的全局尋優(yōu)能力。

圖1 確定性反演算法與隨機(jī)性反演算法對(duì)比

2 方法原理

2.1 量子行為粒子群算法

粒子群算法是一種著名的群理論隨機(jī)性反演算法，在1995年由Kennedy和Eberhart首次提出，靈感源自社會(huì)群體在自然界中搜索食物的過(guò)程，依據(jù)個(gè)體知識(shí)積累和群體信息共享為尋優(yōu)運(yùn)動(dòng)過(guò)程提供指導(dǎo)，不斷更新粒子群參數(shù)信息，最終獲得最優(yōu)解[30]。粒子的速度V和位置矢量X計(jì)算公式分別為

(2)

(3)

計(jì)算粒子群的目標(biāo)函數(shù)ψ的表達(dá)式為

(4)

式中：N為數(shù)據(jù)個(gè)數(shù)；Dobs為擬合數(shù)據(jù)；d為仿真數(shù)據(jù)。

自從粒子群算法提出后，學(xué)者們先后提出了大量的改進(jìn)算法，量子行為粒子群算法是其中之一[31-34]。在量子行為粒子群算法中，采用量子在多維空間中運(yùn)動(dòng)理論為粒子群位置更新提供指導(dǎo)，依據(jù)量子理論中勢(shì)阱概念劃分粒子運(yùn)動(dòng)范圍，使得粒子可以出現(xiàn)在任何存在概率分布的位置上，進(jìn)而能夠克服傳統(tǒng)粒子群算法中早熟收斂問(wèn)題。量子行為粒子群算法的粒子群更新策略(圖2)描述為

圖2 量子行為粒子群算法示意圖

(5)

(6)

(7)

(8)

(9)

(10)

(11)

2.2 α-Trimmed方法橫向約束

依據(jù)各測(cè)點(diǎn)對(duì)應(yīng)的全局最優(yōu)粒子對(duì)于所有其他粒子的空間約束力，提出采用α-Trimmed方法開(kāi)展相鄰測(cè)點(diǎn)的全局最優(yōu)粒子間的橫向模型約束。首先將各測(cè)點(diǎn)對(duì)應(yīng)的全局最優(yōu)粒子整合為一個(gè)最優(yōu)模型矩陣mD,N，其中D為模型中變量數(shù)，N為測(cè)點(diǎn)個(gè)數(shù)。然后依據(jù)需要的橫向約束尺度選擇合適的固定窗寬n，滑動(dòng)窗口讀取矩陣數(shù)據(jù)，對(duì)選中的數(shù)據(jù)進(jìn)行由小到大的排序得到新的序列

m(1∶D,αn+1+j)≤m(1∶D,αn+2+j)≤

…≤m(1∶D,n-αn+j)

(12)

式中：α為裁剪參數(shù)。

選擇裁剪參數(shù)并裁掉序列中前段和后段相同數(shù)目的數(shù)據(jù)，最后將剩下的數(shù)據(jù)取平均值代替窗中心點(diǎn)的值[29]，其表達(dá)式為

(13)

最后以N=5,α=0.25為例，給出了一個(gè)典型α-Trimmed方法處理流程，包括選區(qū)數(shù)據(jù)、排序、裁剪極值和求平均值等步驟。然后，通過(guò)滑動(dòng)窗口對(duì)整體數(shù)據(jù)進(jìn)行處理，實(shí)現(xiàn)橫向約束效果(圖3)。

圖3 α-Trimmed方法的橫向約束示意圖

2.3 量子行為粒子群擬二維反演算法計(jì)算步驟

量子行為粒子群擬二維反演算法計(jì)算步驟(圖4)如下：

圖4 量子行為粒子群算法擬二維反演流程

步驟1，獲取Nobs個(gè)測(cè)點(diǎn)的響應(yīng)數(shù)據(jù)和收發(fā)裝置參數(shù)。

步驟2，設(shè)置反演參數(shù)。

步驟4，計(jì)算第N個(gè)測(cè)點(diǎn)所有粒子的適應(yīng)度,依據(jù)式(5)～(9)更新各測(cè)點(diǎn)粒子群最優(yōu)解和全局最優(yōu)解。

步驟5，根據(jù)式(5)～(8)依次計(jì)算第N個(gè)測(cè)點(diǎn)能量勢(shì)阱的中心位置、能量勢(shì)阱范圍和平均最優(yōu)值，再計(jì)算創(chuàng)造力因子，并根據(jù)式(4)更新粒子群的位置信息。

步驟6，計(jì)算第N個(gè)測(cè)點(diǎn)所有粒子的適應(yīng)度,然后更新粒子群最優(yōu)解和全局最優(yōu)解。

步驟7，如果N

步驟8，利用各測(cè)點(diǎn)全局最優(yōu)解生成最優(yōu)模型矩陣，并設(shè)為外部存檔集。

步驟9，開(kāi)展α-Trimmed方法橫向約束,然后更新各測(cè)點(diǎn)粒子群最優(yōu)解和全局最優(yōu)解,t=t+1。

步驟10，如果t≤T,N=1,返回步驟4，否則輸出全局最優(yōu)解作為最終擬二維反演解。

3 數(shù)值模擬

3.1 一維反演

在一維反演中，仿真模擬的地質(zhì)模型由5層構(gòu)成，由上到下各層電阻率和層厚度參數(shù)分別為[200，20，300，600，100]和[100，80，120，80]。其中發(fā)射波形為25 Hz雙極性方波，電流大小為1 A,線圈半徑為25 m，離地高度為30 m，記錄時(shí)間道為0.001～10.000 ms，接收點(diǎn)位于發(fā)射線圈中心，觀測(cè)磁場(chǎng)強(qiáng)度垂直分量的時(shí)間導(dǎo)數(shù)(dH/dt)響應(yīng)曲線，添加隨機(jī)噪聲后信噪比達(dá)到26 dB。

為檢驗(yàn)量子行為粒子群算法跳出局部最優(yōu)解的能力，對(duì)比了在相同擬合差條件下，兩種主流的正則化約束反演算法(Occam反演算法和聚集反演算法)與傳統(tǒng)粒子群算法和量子行為粒子群算法尋優(yōu)效果對(duì)比。傳統(tǒng)粒子群算法和Occam反演算法的反演結(jié)果[圖5(a)]表明Occam反演算法中提取異常體的邊界信息較為困難，而傳統(tǒng)粒子群算法能夠更好地刻畫(huà)邊界信息。量子行為粒子群算法和聚焦反演算法的反演結(jié)果[圖5(c)]表明聚焦反演算法對(duì)于邊界信息更敏感，但是存在過(guò)多的假異常。在量子行為粒子群算法的一維反演結(jié)果中，在100～190 m深度觀察到低阻層，最小電阻率為8 Ω·m；在310～380 m深度觀察到高阻層，最大電阻率為616 Ω·m。量子行為粒子群算法反演結(jié)果更加接近真實(shí)地質(zhì)模型，有效壓制了反問(wèn)題的非唯一性。在相同擬合差條件下，量子行為粒子群算法搜索精度高于目前主流的反演算法。

圖5 基于5層地電模型的反演擬合圖

在模型參數(shù)設(shè)置中，層數(shù)為50層，厚度參數(shù)固定，粒子數(shù)為200，最大迭代次數(shù)為50。Occam反演算法和聚焦反演算法初始模型為100 Ω·m的半空間模型；傳統(tǒng)粒子群算法初始模型在搜索空間中隨機(jī)生成，自我認(rèn)知系數(shù)和社會(huì)認(rèn)知系數(shù)取值均為2；量子行為粒子群算法中的創(chuàng)造力參數(shù)取值為[0.5,1.0],隨迭代次數(shù)線性遞減變化，實(shí)現(xiàn)反演早期大范圍全局搜索和反演晚期局部精細(xì)搜索。

3.2 擬二維反演

在擬二維反演算法研究中，地質(zhì)模型背景電阻率為上、下兩層,分別為100和300 Ω·m。如圖6(a)所示，在淺部和深部分別設(shè)計(jì)一個(gè)60 m×60 m的正方形高阻異常體和一個(gè)100 m×60 m長(zhǎng)方形低阻異常體。其中，發(fā)射波形為25 Hz雙極性方波，電流為1 A,線圈半徑為15 m，發(fā)射線圈和接收點(diǎn)的離地高度分別為45和50 m，記錄時(shí)間道為0.1～10.0 ms，接收磁感應(yīng)強(qiáng)度垂直分量的時(shí)間導(dǎo)數(shù)(dB/dt)響應(yīng)曲線。在25 m點(diǎn)距的500 m剖面上共包括21個(gè)測(cè)點(diǎn)。響應(yīng)數(shù)據(jù)如圖7(a)所示，淺部異常體位置和深部異常體位置對(duì)應(yīng)的數(shù)據(jù)早期測(cè)道變化都明顯。為了測(cè)試量子行為粒子群算法的魯棒性,在得到的原始數(shù)據(jù)中加入高斯隨機(jī)噪聲,使信噪比達(dá)到26 dB[圖7(b)]。

圖6 電阻率剖面

圖7 多測(cè)道圖

模型參數(shù)反演計(jì)算中分別開(kāi)展了基于L曲線法擬二維聚焦反演和基于α-Trimmed方法的擬二維粒子群反演。α-Trimmed方法經(jīng)過(guò)50次迭代計(jì)算后，各算法的反演結(jié)果如圖6所示。圖6(b)和(c)反演結(jié)果表明，與聚焦反演算法相比，粒子群算法反演結(jié)果中包括更少的假異常，更接近真實(shí)地電模型。對(duì)比圖6(c)和(d)反演結(jié)果，量子行為粒子群算法顯著改善了反演結(jié)果的精度。上述實(shí)驗(yàn)驗(yàn)證了基于量子行為粒子群算法的航空瞬變電磁法擬二維反演技術(shù)的有效性和魯棒性。其中，α-Trimmed方法的參數(shù)N=5，5點(diǎn)平滑既能保證橫向連續(xù)又不會(huì)丟失真實(shí)地電信息；裁剪參數(shù)設(shè)置為α=0.25，去除模型參數(shù)中的極值，保證了反演結(jié)果穩(wěn)定性。

4 結(jié) 語(yǔ)

(1)首先從理論上分析了傳統(tǒng)粒子群算法在高維反演中存在早熟收斂和速度緩慢等問(wèn)題的原因，然后結(jié)合前人改進(jìn)粒子群算法的技術(shù)和航空瞬變電磁反演問(wèn)題特點(diǎn)，提出了基于量子行為粒子群算法與α-Trimmed方法相結(jié)合的擬二維反演技術(shù)，擬克服了早熟收斂和反演速度緩慢問(wèn)題。

(2)經(jīng)過(guò)一維和二維數(shù)值模擬分析，結(jié)果表明：無(wú)論是在相同擬合差條件下，還是相同迭代次數(shù)下，量子行為粒子群算法都能夠提供更接近真實(shí)的地電模型反演結(jié)果，算法的可靠性和魯棒性得到了驗(yàn)證。

(3)目前粒子群反演算法的實(shí)際計(jì)算速度仍大幅落后于傳統(tǒng)的確定性反演算法，下一步計(jì)劃開(kāi)展關(guān)于粒子群反演算法和確定性反演算法的混合算法，將已成熟應(yīng)用于電磁法反演中的確定性反演算法與粒子群算法相結(jié)合，希望能夠獲得更高精度、更快速度的混合算法，在更大的搜索空間中開(kāi)展擬三維航空瞬變電磁反演研究。