劉李楠
(海軍士官學(xué)校,安徽 蚌埠 233012)
計(jì)算機(jī)層析成像就是利用地震波、電容值或電阻率變化等測量數(shù)據(jù)來反演地質(zhì)或物體內(nèi)部結(jié)構(gòu)的物質(zhì)屬性,并逐層剖析繪制其圖像的技術(shù),是目前應(yīng)用最廣泛的工程地質(zhì)探測和工業(yè)成像方法。成像問題的核心最終歸結(jié)于一個(gè)稀疏病態(tài)、混定、不相容方程組的求解[1-6],因此反演算法在整個(gè)層析成像技術(shù)中處于核心地位,反演算法的優(yōu)劣將直接關(guān)系層析成像的成敗。
層析成像中的反演方法可分為線性方法和非線性方法兩種。有效的線性反演方法主要有:Van 和Vorst[7]提出了雙穩(wěn)定共軛梯度法(BICGSTAB),它可以用于求解系數(shù)矩陣不對稱的線性方程組,它利用短遞歸的方法逐步減小殘量,占用內(nèi)存少,但收斂不規(guī)則,在有限精度運(yùn)算時(shí),這種不規(guī)則可能減慢收斂的速度;楊文采[8]提出的帶阻尼因子的LSQR 算法(DLSQR),提高了反演計(jì)算精度,避免了LSQR 算法在數(shù)據(jù)誤差大時(shí)造成的數(shù)值不穩(wěn)定,特別適于求解系數(shù)為大型稀疏矩陣的方程組,與其他迭代方法相比,在解奇異或病態(tài)問題時(shí),顯示出更快的收斂性及更好的可接受的結(jié)果,是目前最常用的線性反演方法;此外還有諸如精細(xì)積分迭代反演算法等方法,都從某些方面提升了反演計(jì)算的性能。非線性反演方法主要有:遺傳算法、模擬退火、粒子群和神經(jīng)網(wǎng)絡(luò)方法[9]等。
由于層析反演方程組是大型稀疏的非線性方程組,采用線性反演算法求解非線性問題不可避免地帶來了解的不適定性等問題,雖然計(jì)算速度快但往往精度較低,穩(wěn)定性也較差。近年來,隨著運(yùn)算平臺計(jì)算能力的大幅提升,以集群智能尋優(yōu)算法為代表的智能方法,因其良好的自適應(yīng)及全局搜索能力在工程技術(shù)中有了越來越多的應(yīng)用,它能更好地排除隨機(jī)干擾帶來的影響,對復(fù)雜介質(zhì)分布重建效果也更為優(yōu)越。
本文提出了一種利用量子行為粒子群算法(QPSO)優(yōu)化BP網(wǎng)絡(luò)解決非線性反演問題的方法,利用QPSO算法中粒子的量子隨機(jī)性優(yōu)化BP 神經(jīng)網(wǎng)絡(luò)的權(quán)值和閾值,提升網(wǎng)絡(luò)反演計(jì)算能力,仿真實(shí)驗(yàn)表明了方法的有效性。
粒子群算法(PSO)是一種群體智能優(yōu)化算法,通過對預(yù)設(shè)初始解的迭代運(yùn)算搜尋最優(yōu)值,具有易于實(shí)現(xiàn)、參數(shù)管理少等優(yōu)點(diǎn)。但從本質(zhì)上說,其不是一種全局收斂算法,受解空間局限、位置進(jìn)化公式等影響,算法中粒子的隨機(jī)性和智能性降低,易陷入局部極值導(dǎo)致魯棒性降低。孫俊等將量子空間引入PSO算法,利用量子理論中的測不準(zhǔn)原理來取代經(jīng)典力學(xué)理論中的粒子運(yùn)動(dòng)行為,極大地增強(qiáng)了解空間中粒子的隨機(jī)性,使得解落入局部極值的概率大大降低,全局搜索能力得到極大增強(qiáng)。
這一過程中,最主要的貢獻(xiàn)即為通過引入量子行為優(yōu)化了粒子位置更新公式[10]:
其中,M為粒子個(gè)數(shù),mbest、pi,d和pg,d分別為平均、個(gè)體和整個(gè)群體的迭代當(dāng)前最優(yōu)位置;α為擴(kuò)張參數(shù)。
BP 算法實(shí)際上就是含有隱層的多層前饋網(wǎng)絡(luò),核心是信息向前傳播而誤差向后傳播,當(dāng)輸出期望值達(dá)不到預(yù)設(shè)輸出時(shí),會利用誤差反向傳播來修改網(wǎng)絡(luò)的權(quán)值和閾值,使得網(wǎng)絡(luò)結(jié)構(gòu)更加優(yōu)化。然而單一的BP神經(jīng)網(wǎng)絡(luò)方法也有其明顯的困難與缺點(diǎn),如收斂性、局部極小以及隱層神經(jīng)元個(gè)數(shù)的選取等。本質(zhì)上BP 神經(jīng)網(wǎng)絡(luò)方法是隨機(jī)的,依賴于空間誤差表面的瞬時(shí)梯度值,因此收斂速度可能較慢甚至可能不收斂,易陷入局部極小,時(shí)間開銷較大,因此在應(yīng)用時(shí)可以綜合運(yùn)用其他方法來優(yōu)化,提升網(wǎng)絡(luò)的收斂性和準(zhǔn)確性。
QPSO算法優(yōu)化BP網(wǎng)絡(luò)反演計(jì)算分為兩部分,如圖1所示,主體右側(cè)是BP網(wǎng)絡(luò)訓(xùn)練部分,左側(cè)是QPSO算法優(yōu)化部分。QPSO算法中粒子代表BP網(wǎng)絡(luò)中的權(quán)值和閾值,利用每個(gè)粒子表示候選解空間內(nèi)的一個(gè)解,解的優(yōu)劣由適應(yīng)度函數(shù)衡量,這個(gè)適應(yīng)度函數(shù)是根據(jù)具體的優(yōu)化目標(biāo)而選取。如本文計(jì)算機(jī)層析成像反演計(jì)算中,粒子代表了BP 網(wǎng)絡(luò)的權(quán)值和閾值,適應(yīng)度函數(shù)可選為BP 網(wǎng)絡(luò)計(jì)算權(quán)值閾值與期望權(quán)值閾值的預(yù)測誤差,適應(yīng)度函數(shù)的優(yōu)化就是尋求誤差最小的過程。左右部分往返調(diào)用,迭代尋找最優(yōu)的粒子,即為BP網(wǎng)絡(luò)權(quán)值和閾值的最優(yōu)設(shè)置,建立了基于QPSO 優(yōu)化的BP 網(wǎng)絡(luò)反演算法(QPSO-BP反演算法)。
圖1 QPSO優(yōu)化BP網(wǎng)絡(luò)反演算法流程圖
采用QPSO-BP 算法進(jìn)行計(jì)算機(jī)層析成像反演計(jì)算,具體步驟及參數(shù)選擇如下:
(1)初始化:采用三層網(wǎng)絡(luò)結(jié)構(gòu),輸入層為監(jiān)測點(diǎn)測量數(shù)據(jù),節(jié)點(diǎn)數(shù)選擇為10;輸出層為反演介質(zhì)網(wǎng)格分布,節(jié)點(diǎn)數(shù)為10;隱層節(jié)點(diǎn)數(shù)根據(jù)經(jīng)驗(yàn)選擇為21。得出權(quán)值為420 個(gè),閾值為41個(gè)。
(2)計(jì)算適應(yīng)度函數(shù)值:選擇BP網(wǎng)絡(luò)訓(xùn)練實(shí)際輸出值和期望輸出值的差來計(jì)算適應(yīng)度函數(shù)值。
(3)集群更新:如果適應(yīng)度不滿足,按QPSO 算法更新,重復(fù)計(jì)算;如果滿足即停止,輸出最優(yōu)權(quán)值和閾值,完成網(wǎng)絡(luò)訓(xùn)練及優(yōu)化。
(4)反演計(jì)算:完成網(wǎng)絡(luò)訓(xùn)練后,將新測量炮點(diǎn)走時(shí)或電阻率數(shù)據(jù)輸入訓(xùn)練好的BP網(wǎng)絡(luò),經(jīng)計(jì)算得到最終反演結(jié)果,完成不同介質(zhì)計(jì)算機(jī)層析成像圖形重建。
為了驗(yàn)證反演算法的效果,建立二維層狀介質(zhì)模型,見圖2,模型尺寸均為100×300m,層狀介質(zhì)速度模型在0~150m深度范圍內(nèi)有5 層,每層厚度均為30m;150~300m 深度范圍內(nèi)有3 層,每層厚度均為50 米;從頂層至底層速度值分別為1800m/s、2000m/s、2200m/s、2400m/s、2600m/s、2800m/s、3000m/s及3200m/s。采用“井間”觀測方式,地表左側(cè)井位置坐標(biāo)為(0,0),發(fā)射源置于左側(cè)井中,個(gè)數(shù)為10 個(gè),第一個(gè)發(fā)射點(diǎn)位于深度10m處,均勻分布,間距為20m,接收點(diǎn)位于右側(cè)井,地表右側(cè)井位置坐標(biāo)為(100,0),接收點(diǎn)數(shù)目及深度分布同發(fā)射點(diǎn)。
圖2 二維層狀介質(zhì)模型圖
采用相同的離散規(guī)則網(wǎng)格,利用不同的反演方法求解計(jì)算。根據(jù)發(fā)射點(diǎn)至接收點(diǎn)射線建立反演方程,此時(shí)反演計(jì)算屬于超定問題無法得到精確解,分別采用PSO-BP 和QPSOBP 算法進(jìn)行求解。參數(shù)設(shè)置中,PSO 算法中粒子種群規(guī)模設(shè)置為80,QPSO 算法中粒子種群規(guī)模為20,最大迭代次數(shù)均為400,粒子速度取值區(qū)間為[-1,1],擴(kuò)張參數(shù)α為0.7;BP網(wǎng)絡(luò)部分,最大迭代代數(shù)設(shè)置為100,不設(shè)置訓(xùn)練終止條件。采用MATLAB 7.0,在硬件Intel(R)Core(TM)i5-3230M CPU 2.6G Hz,8G內(nèi)存的機(jī)器上運(yùn)行,適應(yīng)度進(jìn)化曲線、反演計(jì)算相對誤差如圖3、4所示。
圖3 兩種算法適應(yīng)度進(jìn)化曲線圖
圖4 兩種算法反演重建結(jié)果相對誤差圖
從圖3中可以看出,兩種算法都完成了測量數(shù)據(jù)的反演重建工作。QPSO算法由于對粒子位置更新采用量子行為極大增強(qiáng)了粒子運(yùn)動(dòng)的隨機(jī)性,迅速提升了全局尋優(yōu)能力,在適應(yīng)度計(jì)算迭代到約50次時(shí)已經(jīng)基本到達(dá)了最優(yōu),具有很好的快速尋優(yōu)能力。從誤差圖中看出,除個(gè)別網(wǎng)格介質(zhì)分布誤差較大外,大部分相對誤差在0.03 以內(nèi),基本完成了分層檢測板重建工作,PSO-BP 算法的最大誤差為0.097,最小誤差為9.73e-005,平均誤差為0.0218;而QPSO-BP算法計(jì)算最大誤差為0.065,最小誤差為1.49e-004,平均誤差為0.0098。仿真結(jié)果表明,QPSO-BP算法反演結(jié)果能更準(zhǔn)確反映出模型中各層速度分布,具有很好的反演效果。
層析成像的核心是反演計(jì)算,反演計(jì)算的重點(diǎn)即是對大型、病態(tài)方程組的求解,本文將量子行為粒子群算法引入BP神經(jīng)網(wǎng)絡(luò)權(quán)值和閾值優(yōu)化,然后利用優(yōu)化訓(xùn)練后的網(wǎng)絡(luò)進(jìn)行反演計(jì)算,主要具有以下優(yōu)點(diǎn):
(1)計(jì)算效率高。適應(yīng)度值能快速收斂最優(yōu),通過適應(yīng)度進(jìn)化曲線圖可以看出經(jīng)過不到50次尋優(yōu),適應(yīng)度基本達(dá)到了最優(yōu),收斂速度幾乎是PSO方法的三分之一。
(2)應(yīng)用性強(qiáng)。仿真實(shí)驗(yàn)中,PSO 方法粒子種群數(shù)目為80,而QPSO 方法粒子種群數(shù)目為20 時(shí)即能更好完成任務(wù),計(jì)算過程中參數(shù)選取以及迭代條件不會顯著影響收斂性、收斂速度,具有較好的應(yīng)用前景。