滕明，王典,2，于涵，張敬權(quán)，徐雨歆

1.吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026；2.吉林大學(xué) 地質(zhì)資源立體探測虛擬仿真實驗教學(xué)中心，長春 130026

0 引言

螢火蟲算法(firefly algorithm,簡稱FA)是由Yang[1]于2008年提出的一種仿生隨機(jī)優(yōu)化算法。該算法設(shè)計思路源于對螢火蟲覓食、求偶等行為模擬，通過螢火蟲個體之間相互吸引、移動從而到達(dá)尋優(yōu)目的。與其他群智能算法相比，螢火蟲算法的優(yōu)點表現(xiàn)為：算法結(jié)構(gòu)簡單、所需調(diào)整參數(shù)少以及流程清晰易懂，已在數(shù)值優(yōu)化、工程技術(shù)和聚類分析等[2--6]方面得到成功應(yīng)用。和多數(shù)仿生隨機(jī)算法類似，F(xiàn)A算法也存在一些不足，例如收斂“早熟”現(xiàn)象、后期收斂速度緩慢和易陷入局部極值。針對FA算法存在的不足，許多學(xué)者進(jìn)行了研究與改進(jìn)。主要分為以下兩大類：一類是改進(jìn)FA算法參數(shù)設(shè)置方式。如Lukasik et al.[7]對吸收系數(shù)和移動步長進(jìn)行修改，改進(jìn)后FA算法的求解精度得以提高，但求解速度較慢；歐陽喆等[8]采用動態(tài)方式調(diào)整步長；Gandomi et al.[9]采用混沌序列修改FA算法參數(shù)設(shè)置，提升FA算法的求解精度和收斂速度；徐麗華等[10]針對FA算法“早熟”收斂現(xiàn)象，提出了一種變尺度混沌光強(qiáng)吸收系數(shù)調(diào)整策略的螢火蟲優(yōu)化算法；劉暢等[11]利用指數(shù)分布等改進(jìn)吸引項。另一類是將FA與粒子群算法(PSO)、遺傳算法(GA)和差分演化算法(DE)等優(yōu)化算法相結(jié)合，提高FA算法的性能。曹秀爽[12]提出當(dāng)FA算法處于全局搜索時，將模擬退火搜索機(jī)制融入FA算法，當(dāng)FA算法處于局部搜索時，把回火策略引入FA算法，以此解決螢火蟲易陷入局部極值和后期收斂速度慢等問題； Xia et al.[13]基于螢火蟲算法和粒子群算法，提出一種混合優(yōu)化算法。

筆者提出基于混沌搜索和動態(tài)步長修改的螢火蟲算法(CAFA)。在種群初值選取時，利用混沌序列代替隨機(jī)方式對螢火蟲種群進(jìn)行初始化?；煦缡谴_定性非線性系統(tǒng)內(nèi)存在的特殊隨機(jī)現(xiàn)象，具有以下特性：有界性、隨機(jī)性和遍歷性。有界性是指混沌變量運(yùn)動軌跡始終局限于一個確定區(qū)域內(nèi)，此區(qū)域被稱為混沌吸引域；隨機(jī)性是指某一時刻的混沌變量會受到前一狀態(tài)影響而確定出現(xiàn)；遍歷性是指混沌變量運(yùn)動軌跡不會停留于某一狀態(tài)而是不重復(fù)地遍歷空間區(qū)域中所有可能存在狀態(tài)?；煦绯跏蓟岣吡朔N群多樣性和搜索空間分布質(zhì)量，避免算法陷入局部最優(yōu)。另外，采用一種新型的雙曲線遞減步長代替?zhèn)鹘y(tǒng)的固定步長，實現(xiàn)了提高算法后期收斂速度、尋優(yōu)精度以及算法穩(wěn)定性的目的。

波阻抗參數(shù)是聯(lián)系地震、測井及地質(zhì)信息的紐帶。通過波阻抗反演可以把界面型的地震剖面轉(zhuǎn)換成巖層型的波阻抗剖面，可直接將地震資料與測井對比進(jìn)行儲層巖性解釋和物性分析，因此波阻抗反演技術(shù)在儲層參數(shù)估算[14]、油藏預(yù)測[15]等研究工作中發(fā)揮關(guān)鍵的作用。地震資料波阻抗反演是一個非線性反演問題，傳統(tǒng)的地震波阻抗反演常用線性化算法，如最速下降法、牛頓法和共軛梯度法等，此類方法計算效率高、穩(wěn)定性好，但過于依賴初始模型且算法收斂速度慢、易陷入局部極值等缺點。近些年，遺傳算法、粒子群算法和模擬退火算法等智能優(yōu)化算法已經(jīng)在波阻抗反演中得到廣泛應(yīng)用。 聶茹等[16]將免疫算法和遺傳算法相結(jié)合應(yīng)用于波阻抗反演；許永忠等[17]將粒子群算法應(yīng)用到火山巖的波阻抗反演，將反演結(jié)果與概率神經(jīng)網(wǎng)絡(luò)方法相對比，證明了粒子群算法在波阻抗反演的有效性。螢火蟲算法在波阻抗反演方向研究較少，筆者將CAFA算法應(yīng)用于波阻抗反演，通過理論模型試算分析CAFA算法的可行性。

1 傳統(tǒng)螢火蟲算法(FA)

1.1 生物學(xué)機(jī)制

自然界中大多數(shù)螢火蟲具有發(fā)光行為，發(fā)出的閃光信號可以被一定范圍內(nèi)其他螢火蟲所感知，螢火蟲借助此獨特的發(fā)光行為進(jìn)行覓食、求偶等群體行為。如圖1所示，將螢火蟲個體用搜索空間范圍內(nèi)紅點表示，螢火蟲的熒光亮度取決于搜索范圍內(nèi)紅點所在位置的目標(biāo)函數(shù)值。吸引度與亮度有關(guān)，越亮的螢火蟲擁有越強(qiáng)的吸引力，吸引視線范圍內(nèi)亮度弱的螢火蟲向其所在位置方向移動，若發(fā)光亮度相同，則螢火蟲各自隨機(jī)移動。把螢火蟲個體之間的不斷飛行和吸引的行為模擬成算法迭代和優(yōu)化過程，舍棄螢火蟲發(fā)光的某些生物學(xué)上意義，對閃光特性做出理想化處理[17]。

圖1 螢火蟲算法仿生模擬圖Fig.1 Bionic simulation diagram of firefly algorithm

①假設(shè)螢火蟲沒有雌雄之分。螢火蟲相互吸引與螢火蟲的性別無關(guān)。

②螢火蟲之間的吸引度與熒光亮度和距離有關(guān)。亮度高的螢火蟲會吸引周圍亮度弱的螢火蟲，隨著距離增大吸引度逐漸減小，亮度相同或是亮度比鄰域其他螢火蟲高的個體在自身周圍隨機(jī)運(yùn)動。

③螢火蟲的亮度與具體研究問題的目標(biāo)函數(shù)有關(guān)。

1.2 FA算法描述

螢火蟲算法的核心思想是熒光亮度弱的螢火蟲會被亮度強(qiáng)的螢火蟲吸引，并根據(jù)位置更新公式進(jìn)行位置更新，經(jīng)過多次迭代最終收斂到最優(yōu)值附近。下面定義相對亮度、吸引力、位置更新等核心公式。

相對熒光亮度[17]:

(1)

式中：I0表示自身(r=0處) 的熒光亮度，與目標(biāo)函數(shù)值有關(guān)，熒光亮度越高目標(biāo)值越優(yōu)；γ為光強(qiáng)吸收系數(shù)，衡量熒光亮度在傳播路徑中的減弱程度；rij為螢火蟲i和j之間的歐氏距離，對于其他優(yōu)化問題，距離公式依據(jù)具體問題進(jìn)行定義。

吸引度[17]：

(2)

式中：β0表示光源處(r=0處)的吸引度；γ為光強(qiáng)吸收系數(shù)；rij為螢火蟲i到螢火蟲j的歐氏距離[17]：

(3)

位置更新公式[17]：

Xi(t+1)=Xi(t)+β(rij){Xj(t)-Xi(t)}

+α(rand-1/2)

(4)

Xi(t+1)=Xi(t)+α(rand-1/2)

(5)

式中：Xi,Xj是t時刻螢火蟲i和j的位置；α為移動步長；β為吸引度系數(shù)；rand為0-1的隨機(jī)數(shù)。把t時刻螢火蟲Xi,Xj代入公式(1)并比較二者熒光亮度的大小，如果I(Xi)

在公式(4)中，Xi(t)表示螢火蟲當(dāng)前位置對更新位置的影響，起到平衡全局和局部搜索的能力；β(rij){Xj(t)-Xi(t))}反映群體信息交流的影響，是螢火蟲之間的信息共享，取決于吸引力的大??；α(rand-1/2)是帶有特定系數(shù)的隨機(jī)擾動項，避免螢火蟲陷入局部極值點。通過查閱文獻(xiàn)以及測試大量復(fù)雜函數(shù)發(fā)現(xiàn)FA算法主要存在兩點不足之處：①傳統(tǒng)FA算法的位置更新公式采用固定步長向較優(yōu)個體移動，如果被吸引螢火蟲與較優(yōu)個體之間距離小于固定步長因子，被吸引螢火蟲在移動時會躍過較優(yōu)個體，并隨著迭代進(jìn)行反復(fù)在較優(yōu)個體附近震蕩；如果步長較小，則需要多次的迭代計算，導(dǎo)致收斂速度降低。②FA采用隨機(jī)方式對螢火蟲個體位置進(jìn)行初始化，這樣會造成螢火蟲位置分布不均勻，導(dǎo)致個體間差異較大，降低種群的多樣性，算法易陷入局部最優(yōu)點。

2 改進(jìn)的螢火蟲算法(CAFA)

2.1 動態(tài)步長

在FA 算法搜索過程中，如果固定步長較大會使搜索精度降低且收斂后期出現(xiàn)震蕩現(xiàn)象。如果固定步長較小導(dǎo)致搜索速度降低，但求解精度提高。CAFA算法利用動態(tài)步長代替固定步長，動態(tài)步長是一種隨迭代次數(shù)增加，步長不斷變化的函數(shù)。在算法迭代初期，利用較大步長促使螢火蟲個體加速向較優(yōu)值附近移動，隨著迭代次數(shù)增加步長會逐漸減小直至零值，避免了算法后期被吸引螢火蟲個體與較優(yōu)個體之間距離小于固定步長值而出現(xiàn)振蕩現(xiàn)象。以四峰函數(shù)為例，當(dāng)坐標(biāo)為(0, 0)、 (0， -4)時，函數(shù)取得極大值fmax=2。圖2經(jīng)過15次迭代，對30只螢火蟲進(jìn)行了450次估算得到收斂效果圖，圖2a自適應(yīng)步長螢火蟲成功搜索到全局最優(yōu)值fmax=2，然而圖2b固定步長FA算法收斂精度為1.957，增加迭代次數(shù)固定步長FA算法收斂精度未能提高。

圖2 動態(tài)步長(a)及固定步長(b) Fig.2 Dynamic step(a) and fixed step(b)

在理想情況下，群體中所有個體隨著飛行次數(shù)的增加逐漸向最優(yōu)值附近移動，最終收斂到一點。

因此，對于任意的兩個螢火蟲i和j可以得到[18]：

(6)

(7)

式(6)表示搜索空間內(nèi)螢火蟲個體收斂到一點(最優(yōu)值)，式(7)表示收斂結(jié)束后個體(解)不再發(fā)生移動(變化)。結(jié)合上述公式可以得到[18]：

(8)

從式(8)可以得到，隨著算法進(jìn)行多次迭代，移動步長α趨于零。在傳統(tǒng)FA中，設(shè)置步長α是靜態(tài)的，它不能真正模擬算法搜索過程。通常對FA算法來說，α值較大有利于探索新的搜索空對全局最優(yōu)收斂沒有幫助。如果步長α值較小，則結(jié)果相反，因此移動步長α對算法的探索和收斂有著較大影響。本文使用雙曲線方程設(shè)計動態(tài)步長調(diào)整方案。雙曲線遞減動態(tài)步長：

(9)

式中：t為當(dāng)前迭代次數(shù)；Tmax為最大迭代次數(shù)；k為控制動態(tài)步長減小速率的參數(shù)。雙曲線遞減步長的動態(tài)遞減性質(zhì)如圖3所示。從圖3可以看出，步長在初期階段較大，然后隨著迭代次數(shù)的增加逐漸減小，這有助于算法平衡全局探索和局部搜索的能力。表1列舉其他學(xué)者設(shè)計的各種動態(tài)步長[19--22]和筆者提出的雙曲線遞減步長。以上述四峰函數(shù)作為目標(biāo)函數(shù)，測試各種動態(tài)步長的收斂效果。四峰函數(shù)全局最優(yōu)值為fmax=2。這里定義差值參數(shù)DP，它是全局最優(yōu)值與自適應(yīng)步長螢火蟲算法求解的差值，DP越接近于0值，表示算法求解精度越高。在表1中，步長α=0.4是FA算法的固定步長，其余的步長均為動態(tài)遞減步長且滿足收斂公式(8)。測試結(jié)果表明動態(tài)步長在求解精度和收斂速度方面明顯優(yōu)于固定步長。本文提出的雙曲線遞減步長螢火蟲算法迭代次數(shù)最少、收斂速度最快，在收斂性能方面明顯優(yōu)于上述文獻(xiàn)改進(jìn)的動態(tài)步長。

圖3 雙曲線動態(tài)步長收斂效果圖Fig.3 Convergence effect diagram of hyperbolic dynamic step

表1 測試自適應(yīng)步長收斂效果

2.2 混沌序列優(yōu)化FA

傳統(tǒng)FA算法采用隨機(jī)方式對螢火蟲種群個體位置加以初始化，導(dǎo)致螢火蟲個體位置分布不均勻、個體之間差異性較大及種群的多樣性降低，僅有少數(shù)優(yōu)秀個體在移動和迭代過程中被保留下來。為了解決這個問題，筆者將Logistic混沌序列引入傳統(tǒng)FA算法初始化階段，并結(jié)合具體問題相應(yīng)地修改從混沌空間到搜索空間的映射方式?；煦缧蛄芯哂歇毺匦再|(zhì)即有界性、隨機(jī)性、規(guī)律性和遍歷性，混沌變量類似于隨機(jī)變量，但能夠按照指定規(guī)律在搜索范圍內(nèi)不重復(fù)地遍歷所有狀態(tài)以確保生成搜索能力較強(qiáng)、分布較廣且均勻的初始種群，增強(qiáng)初始種群多樣性及其在搜索空間分布質(zhì)量，達(dá)到避免算法陷入局部極值點，提高算法求解精度。

采用Logistic映射函數(shù)生成混沌序列為

(10)

混沌序列初始化螢火蟲種群的基本步驟：

①隨機(jī)對D維空間內(nèi)螢火蟲進(jìn)行初始化得到一個混沌變量；

②按照式(10)遞推M-1次，最終得到的M個混沌變量；

③將產(chǎn)生混沌變量按式(11)映射到目標(biāo)函數(shù)的搜索空間，從而得到M只螢火蟲初始位置，公式如下：

xi,d=min{|LB|,|UB|}·yi,d

(11)

式中：LB和UB分別表示搜索空間第d維的下限和上限，yi,d是第i只螢火蟲相對應(yīng)的第d維混沌變量，xi,d是第i只螢火蟲在搜索空間第d維的坐標(biāo)值。圖4混沌序列初始化螢火蟲算法的流程圖；圖5直觀得出采用混沌序列優(yōu)化的螢火蟲算法，其種群分布的均勻性要明顯優(yōu)于隨機(jī)方式初始化的螢火蟲種群。引入霍普金斯統(tǒng)計量[23](Hopkins statistic)定量衡量初始種群分布均勻性，如果霍普金斯統(tǒng)計量H接近于0.5表示數(shù)據(jù)分布均勻，如果H接近于1表示數(shù)據(jù)有聚類趨勢，定義霍普金斯統(tǒng)計量H：

(12)

分別計算兩種初始化螢火蟲的方式的霍普金斯統(tǒng)計量H，采用混沌序列初始化的H1=0.485 6，而隨機(jī)方式初始化的H2=0.889 4，證明在種群分布均勻性方面上混沌序列初始化優(yōu)于傳統(tǒng)隨機(jī)方式。

圖4 混沌序列初始化螢火蟲算法的流程圖Fig.4 Flowchart of chaotic sequence initialization firefly algorithm

圖5 隨機(jī)方式與混沌搜索初始化螢火蟲位置對比圖Fig.5 Comparison chart of firefly position initialized by random method and chaotic search

圖6 測試函數(shù)的三維圖Fig.6 Three dimensional diagram of test function

2.3 改進(jìn)螢火蟲算法和傳統(tǒng)螢火蟲算法比較

3 CAFA反演波阻抗

3.1 波阻抗反演

以褶積模型為基礎(chǔ)進(jìn)行波阻抗反演，地震記錄的數(shù)學(xué)模型為：

S(t)=W(t)×R(t)

(13)

式中:S為地震記錄序列；W為地震子波序列是地震記錄的基本單位；R為地震反射系數(shù)序列，第i層的反射系數(shù)如下：

(14)

式中：ρi、vi是第i層的密度和速度，zi是第i層的波阻抗。由公式(13)和(14)可知，地震記錄與波阻抗為非線性關(guān)系，因此求取波阻抗是一個求解非線性組合優(yōu)化問題。目前模擬退火算法[24]、遺傳算法[25]、粒子算法[26]和差分進(jìn)化算法[27]等智能優(yōu)化算法已經(jīng)在反演波阻抗方面得到了成功應(yīng)用，為螢火蟲算法在波阻抗反演研究提供一些理論支持。筆者以改進(jìn)螢火蟲算法(CAFA)研究為基礎(chǔ)，將CAFA算法和地震記錄波阻抗反演相結(jié)合，驗證CAFA算法在地球物理反演方面的可行性。使用最小二乘法建立螢火蟲算法反演波阻抗的目標(biāo)函數(shù)，建立的目標(biāo)函數(shù)如下：

(15)

式中：M為空間采樣點；W為Ricker子波；R(t)為所求的反射序列；S(t)為實際地震記錄。CAFA經(jīng)過多次迭代后得到的反射系數(shù)和Ricker子波褶積與實際地震記錄的殘差平方和達(dá)到全局極小值，再利用CAFA算法得到反射系數(shù)來遞推各層的波阻抗。在反射系數(shù)遞推波阻抗過程中誤差會不斷積累，若反射系數(shù)存在較大誤差會導(dǎo)致遞推的波阻抗曲線越來越偏離實際的波阻抗曲線。

3.2 理論模型試算

為檢驗CAFA算法反演波阻抗的可行性，設(shè)計一個五層水平地層理論模型(圖 8)，螢火蟲種群規(guī)模n=50,光強(qiáng)吸收系數(shù)β=1,初始步長α=0.2，自適應(yīng)步長選取雙曲線遞減步長。圖8 理論模型各層介質(zhì)速度為1 600、2 000、2 300、1 900、2 600 m/s，各層介質(zhì)密度為1 500、1 900、2 000、1 800、2 400 kg/m3。圖9水平地層介質(zhì)波阻抗模型，該模型各層介質(zhì)波阻抗為2.40、4.60、3.06、1.82、3.42，單位為106kg·m-2·s-1。圖10是水平地層模型正演得到單道合成地震記錄，其中子波選取Ricker子波，子波主頻為50 Hz，采樣率為2 ms。從圖11可知，CAFA反演波阻抗與模型波阻抗吻合程度更好，能夠反演理論模型的各個層位波阻抗變化情況。圖12 中CAFA算法和FA算法分別反演合成地震記錄與模型合成地震記錄都具有較好吻合，且CAFA算法的求解精度明顯高于FA算法(CAFA算法求解精度為0.000 284、FA算法求解精度為0.005 320)。分別在原始合成地震記錄中加入3%、5%和10%隨機(jī)噪聲驗證CAFA算法的抗噪能力。圖13、圖14分別是波阻抗和合成地震記錄反演的結(jié)果，含3%、5%的隨機(jī)噪聲時，反演結(jié)果誤差較小能夠真實的反映有效地層信息。含10%的隨機(jī)噪聲時，誤差較大，但反演的總體變化趨勢仍然正確，這表明CAFA具有較強(qiáng)的抗噪能力。

圖8 五層水平地層Fig.8 Five horizontal strata

圖9 模型的波阻抗序列Fig.9 Wave impedance sequence of model

圖10模型的合成地震記錄Fig.10 Synthetic seismic record of model

圖11 CAFA和FA反演波阻抗Fig.11 Inversion wave impedance of CAFA and FA

圖12 CAFA和FA反演合成地震記錄Fig.12 Inversion synthetic seismic records of CAFA and FA

圖13 反演波阻抗效果(含3%，5%，10%隨機(jī)噪聲)Fig.13 Inversion of impedance effect(including 3%, 5%, and 10% random noise)

圖14 反演合成地震記錄(含3%，5%，10%隨機(jī)噪聲)Fig.14 Inverted synthetic seismic records(including 3%, 5%, and 10% random noise)

4 結(jié)論

(1)雙曲線遞減步長相比于傳統(tǒng)固定步長，不但提高算法求解精度(7.684 9e-07)，而且大幅度地提升收斂速度(1.3 s)，避免了算法收斂后期反復(fù)振蕩現(xiàn)象。

(2)采用混沌序列替換隨機(jī)方式初始化螢火蟲種群能夠提高初始種群分布質(zhì)量，有效地解決算法易陷入局部極值問題。

(3)利用CAFA算法反演水平地層的波阻抗模型，得到誤差較小的反演波阻抗和反演合成地震記錄。在理論模型中分別加入不同強(qiáng)度的隨機(jī)噪聲，驗證本文提出的CAFA算法具有較好的抗噪能力。