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

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

0 引 言

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

隨著近些年來計算科學(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等率先將其引入到地球物理反演中，成功地實現(xiàn)了垂直電測深的反問題求解[21]。最近十幾年，眾多學(xué)者先后開展了基于粒子群算法的反演研究工作，極大推動了粒子群算法在地球物理領(lǐng)域的發(fā)展，但是由于粒子群算法在求解非線性映射關(guān)系下的高維反問題時，仍然存在早熟收斂和收斂速度慢等問題，所以限制了粒子群算法在電磁法二、三維反演中的發(fā)展[15,22-24]。

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

1 航空瞬變電磁法反問題

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

d=F(m)

(1)

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

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

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

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

2 方法原理

2.1 量子行為粒子群算法

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

(2)

(3)

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

(4)

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

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

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

(5)

(6)

(7)

(8)

(9)

(10)

(11)

2.2 α-Trimmed方法橫向約束

依據(jù)各測點(diǎn)對應(yīng)的全局最優(yōu)粒子對于所有其他粒子的空間約束力，提出采用α-Trimmed方法開展相鄰測點(diǎn)的全局最優(yōu)粒子間的橫向模型約束。首先將各測點(diǎn)對應(yīng)的全局最優(yōu)粒子整合為一個最優(yōu)模型矩陣mD,N，其中D為模型中變量數(shù)，N為測點(diǎn)個數(shù)。然后依據(jù)需要的橫向約束尺度選擇合適的固定窗寬n，滑動窗口讀取矩陣數(shù)據(jù)，對選中的數(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為例，給出了一個典型α-Trimmed方法處理流程，包括選區(qū)數(shù)據(jù)、排序、裁剪極值和求平均值等步驟。然后，通過滑動窗口對整體數(shù)據(jù)進(jìn)行處理，實現(xiàn)橫向約束效果(圖3)。

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

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

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

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

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

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

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

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

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

步驟7，如果N

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

步驟9，開展α-Trimmed方法橫向約束,然后更新各測點(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，記錄時間道為0.001～10.000 ms，接收點(diǎn)位于發(fā)射線圈中心，觀測磁場強(qiáng)度垂直分量的時間導(dǎo)數(shù)(dH/dt)響應(yīng)曲線，添加隨機(jī)噪聲后信噪比達(dá)到26 dB。

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

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

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

3.2 擬二維反演

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

圖6 電阻率剖面

圖7 多測道圖

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

4 結(jié) 語

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

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

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