龐 珂,李 劍,蘇新彥,劉曉佳,魏曉曼,孫袖山

(1.中北大學(xué)省部共建動(dòng)態(tài)測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室,山西 太原 030051;2.中北大學(xué)信息探測(cè)與處理山西省重點(diǎn)實(shí)驗(yàn)室,山西 太原 030051)

0 引言

地下淺層震源定位技術(shù)是實(shí)現(xiàn)地下淺層爆炸、爆破、礦產(chǎn)資源勘探的關(guān)鍵技術(shù)手段,這種技術(shù)也經(jīng)常應(yīng)用在地震災(zāi)害評(píng)估以及應(yīng)急救援響應(yīng)的過程中[1-2]。與地下深層的震源定位相比,淺層空間因?yàn)榈刭|(zhì)結(jié)構(gòu)復(fù)雜、速度信息未知等原因,造成了在大區(qū)域淺層震源搜索過程中的難點(diǎn)和熱點(diǎn)問題。

目前常用的震源搜索方法包括模擬退火法、網(wǎng)格搜索法、粒子群優(yōu)化算法等[3]。文獻(xiàn)[4]結(jié)合模擬退火算法與粒子群算法提出了改進(jìn)的全局搜索算法。文獻(xiàn)[5]提出了一種基于粒子群優(yōu)化(PSO)算法的微震源智能定位方法,消除了傳統(tǒng)微震源定位時(shí)因縱波速度誤差帶來的不利影響。文獻(xiàn)[6]將粒子群算法用于提高虛擬場(chǎng)優(yōu)化(VFOM)的震源定位效率中,并與遺傳算法(GA)對(duì)比得出PSO可以為VFOM提供比GA更好的定位精度和計(jì)算效率。文獻(xiàn)[7]利用PSO算法估計(jì)震源的初始位置,然后通過牛頓迭代搜索方法確定震源位置,該方法能有效解決牛頓迭代法的初值選取問題,提高定位精度。文獻(xiàn)[8]提出了一種基于粒子群算法的巖體微震源分層定位方法,用于解決經(jīng)典法定位時(shí)速度模型反演不準(zhǔn)確和采用聯(lián)合法定位時(shí)解不唯一的問題,提高了收斂速度和解的穩(wěn)定性。文獻(xiàn)[9]提出一種基于遺傳算法與模擬退火算法的TDOA定位算法。文獻(xiàn)[10]提出了基于直接搜索算法的微地震震源快速網(wǎng)格搜索定位法,該方法無需劃分網(wǎng)格大小,在減少搜索次數(shù)的同時(shí),提高了計(jì)算效率和定位精度。文獻(xiàn)[11]用振幅疊加網(wǎng)格搜索法搜索區(qū)域中能量最大的位置,這種方法有效提高了計(jì)算效率,在定位精度上比單純的網(wǎng)格搜索法有了很大提高。文獻(xiàn)[12]在低信噪比下利用QPSO在區(qū)域中搜索陣列的最佳位置。文獻(xiàn)[13]提出了一種基于QPSO聯(lián)合走時(shí)-偏振角度信息的定位算法,有效降低了空間搜索范圍的局限性,適合大區(qū)域范圍的單震源定位。但是,上述方法在提高震源搜索精度的同時(shí)也存在一定的不足,例如:粒子群算法的搜索性能對(duì)初始位置的選取有一定的依賴性,無法對(duì)定位區(qū)域形成有效的覆蓋,導(dǎo)致算法不易收斂;而且粒子群算法本身也容易出現(xiàn)“早熟”現(xiàn)象,從而導(dǎo)致群體陷入局部收斂,造成定位誤差。

針對(duì)上述問題,本文提出一種利用DIR-QPSO聯(lián)合SRP的震源定位方法。采用Logistic混沌模型對(duì)初始種群進(jìn)行優(yōu)化,使得初始種群的分布相對(duì)均勻;然后通過DIR-QPSO算法引入輔助粒子群對(duì)主粒子群進(jìn)行約束,防止群體陷入局部收斂;同時(shí)以SRP為基礎(chǔ)構(gòu)建震源區(qū)域中的能量信息,并將其作為適應(yīng)度函數(shù),以此避免傳感器布設(shè)和拾取初至?xí)r刻對(duì)定位精度的影響。

1 基于SPR的震源定位模型

波束合成的震源定位算法是對(duì)于傳感器接收到的信號(hào)進(jìn)行時(shí)延補(bǔ)償-求和,得到該點(diǎn)的輸出功率,通過波束合成得到的輸出功率又稱為可控響應(yīng)功率(SRP, steered-response power),對(duì)區(qū)域內(nèi)的每個(gè)點(diǎn)進(jìn)行波束合成計(jì)算,輸出功率最大的點(diǎn)即震源所在的位置。波束形成示意圖如圖1所示。

圖1 波束形成示意圖Fig.1 Schematic diagram of beamforming

一個(gè)傳感器陣列有M個(gè)陣元,第m個(gè)傳感器接收到的信號(hào)可以表示為

xm(t)=hm(t)*s(t)+wm(t),

(1)

式(1)中,s(t)是震源信號(hào),hm(t)是震源到第m個(gè)傳感器的沖擊響應(yīng),“*”表示卷積,wm(t)為第m個(gè)傳感器接收到的噪聲,m=1,2,…,M。

對(duì)于空間中某一點(diǎn)q,該陣列的波束形成的定義為

(2)

式(2)中,τ(q,m)為第m個(gè)傳感器接收到震源信號(hào)的延時(shí)。波束形成的輸出功率稱為可控響應(yīng)功率(SRP),其定義式為

式(3)中,Y(ω,q)為Y(q)的傅里葉變換,SRP可以表示為

(4)

在實(shí)際情況中,傳感器接收到的信號(hào)和濾波器都是能量有限的,所以式(4)的積分是收斂的,求和與積分運(yùn)算的順序可以交換,得

又因?yàn)棣觤l(q)=τm(q)-τl(q),表示第m個(gè)傳感器與第l個(gè)傳感器接收位置q處信號(hào)的時(shí)間差,故式(4)可表示為

第m個(gè)傳感器與第l個(gè)傳感器接收到的信號(hào)的互相關(guān)(GCC)定義表示為

故,SRP可以表示為所有傳感器對(duì)的GCC之和可表示為

(8)

去掉其中的自相關(guān)值,可表示為

(9)

震源q所在位置為區(qū)域中可控響應(yīng)功率最大值處,即

2 震源搜索算法

2.1 QPSO算法

量子粒子群算法在粒子群優(yōu)化算法(PSO, particle swarm optimization)的基礎(chǔ)上取消了粒子的移動(dòng)方向?qū)傩?這樣在增加了粒子位置隨機(jī)性的同時(shí)消除了與粒子之前運(yùn)動(dòng)的關(guān)系。因?yàn)樵赒PSO算法中粒子的位置和速度在量子空間中不能一起確定,所以用波函數(shù)代替粒子位置和速度,通過蒙特卡羅方法求解波函數(shù),從而得出粒子位置和速度,QPSO算法模型如圖2所示。

圖2 QPSO算法模型示意圖Fig.2 Schematic diagram of the QPSO algorithm model

種群的局部吸引點(diǎn)pi,j=(pix,piy,piz)是以粒子的個(gè)體極值pbest和全局極值gbest來決定的,如式(11)所示:

pi,j(t)=(φ1·pbesti,j(t)+φ2·gbesti,j(t))/(φ1+φ2),

(11)

式(11)中,φ1和φ2為(0,1)內(nèi)的隨機(jī)數(shù)。

QPSO的粒子位置的更新公式為

Xi,j(t+1)=pi,j(t)±β|mbest-Xi,j(t)|ln[1/ui,j(t)],

(12)

式(12)中,ui,j(t)為(0,1)內(nèi)的隨機(jī)數(shù);Xi,j(t+1)為第t+1代粒子群的位置;β為收縮擴(kuò)張系數(shù),一般在(0,1)內(nèi)取值;mbest為種群平均最好位置,計(jì)算公式為

(13)

2.2 DIR-QPSO算法

利用DIR-QPSO算法進(jìn)行迭代時(shí),首先在給定的區(qū)域中初始化一個(gè)粒子種群S,將這個(gè)種群劃分為主粒子群S1和輔粒子群S2(S=S1∪S2)兩個(gè)部分,其中,主粒子群的局部吸引點(diǎn)pi,j是以主粒子的個(gè)體極值pbest和全局極值gbest來決定的,如式(14)所示:

pi,j(t)=(φ1·pbesti,j(t)+φ2·gbestj(t))/(φ1+φ2),

(14)

式(14)中,φ1和φ2為(0,1)內(nèi)的隨機(jī)數(shù)。

輔粒子群的局部吸引點(diǎn)qi,j是以主粒子群局部吸引點(diǎn)pi,j和輔粒子群全局最優(yōu)位置gbestj決定的,其具體表現(xiàn)形式為

qi,j(t)=[1-φi,j(t)]·pi,j(t)+φi,j(t)·gbestj(t),

(15)

式(15)中,φi,j(t)為[0,1]范圍內(nèi)的隨機(jī)數(shù)。

DIR-QPSO的粒子位置的更新公式如下:

式(16)中,ui,j(t)為(0,1)內(nèi)的隨機(jī)數(shù),Xi,j(t+1)為第t+1代主群粒子的位置,Yi,j(t+1)為第t+1代輔群粒子的位置,α為收縮擴(kuò)張系數(shù),一般取值在0～1之間,Rj(t)為種群隨機(jī)最好位置,它是在種群最好位置gbest與平均最好位置Mbest之間的一個(gè)參數(shù),能夠在兩者之間作出有效的協(xié)調(diào),具體形式為

Rj(t)=rj(t)gbest+(1-rj(t))Mbest,

(17)

式(17)中,rj(t)為[0,1]范圍內(nèi)的隨機(jī)數(shù)。

2.3 基于Logistic混沌模型的粒子初始化

無論是DIR-QPSO還是QPSO算法在初始化種群時(shí),常用隨機(jī)均勻分布的方法,但在這個(gè)過程中,可能會(huì)出現(xiàn)粒子重復(fù),空間覆蓋的范圍低,容易造成隨機(jī)生成的粒子質(zhì)量較差,從而導(dǎo)致算法定位的效果差,收斂區(qū)域錯(cuò)誤[14]。

故引入混沌模型為L(zhǎng)ogistic模型搜索的方式,如式(18)所示:

y(k+1)=α·y(k)·(1-y(k)),

(18)

其中,y是粒子x經(jīng)過式(19)所示的映射關(guān)系變換后的變量,k為迭代次數(shù),α為混沌映射程度。

2.4 構(gòu)建基于SRP的適應(yīng)度函數(shù)

通過式(9)和式(10)構(gòu)建的適應(yīng)度函數(shù)為

將式(20)作為DIR-QPSO算法的適應(yīng)度函數(shù),適應(yīng)度最大的粒子位置就是當(dāng)前種群的最優(yōu)位置gbestj,通過適應(yīng)度函數(shù)評(píng)價(jià)當(dāng)前粒子位置是否是實(shí)際震源位置的重要指標(biāo)。

2.5 定位模型解算流程

1) 初始化種群

首先設(shè)定震源的搜索范圍,設(shè)置種群規(guī)模為120、空間維數(shù)為3、迭代測(cè)試為500,根據(jù)式(18)、式(19)所示的Logistic混沌模型隨機(jī)生成初始粒子群,將這個(gè)初始種群劃分為兩個(gè)種群規(guī)模為60的主粒子群S1和輔粒子群S2。

2) 更新兩個(gè)種群的最優(yōu)震源位置

分別根據(jù)式(20)構(gòu)建的適應(yīng)度函數(shù)計(jì)算得到震源群中適應(yīng)度值最大的位置,設(shè)為當(dāng)前震源種群最優(yōu)震源gbest1和gbest2。

3) 更新全局最優(yōu)位置

判斷當(dāng)前計(jì)算得到的兩個(gè)最優(yōu)震源gbest1和gbest2哪個(gè)震源處的適應(yīng)度最大,則將這個(gè)適應(yīng)度最大的位置作為當(dāng)前的全局最優(yōu)Gbest。

4) 重復(fù)迭代

當(dāng)未達(dá)到迭代次數(shù)時(shí),通過式(16)再次分別更新兩個(gè)種群中最優(yōu)的震源位置,重復(fù)進(jìn)行步驟2),步驟3)。

5) 結(jié)果輸出

當(dāng)達(dá)到迭代次數(shù)時(shí),將全局震源最優(yōu)位置Gbest輸出,并設(shè)置為最優(yōu)化的震源位置。

3 實(shí)驗(yàn)驗(yàn)證

3.1 實(shí)驗(yàn)方案

實(shí)驗(yàn)采用本單位研發(fā)的分布式地下震源定位系統(tǒng)進(jìn)行定位,實(shí)驗(yàn)區(qū)域?yàn)樽匀煌两橘|(zhì),實(shí)驗(yàn)現(xiàn)場(chǎng)如圖3所示。通過與文獻(xiàn)[13]中提出的QPSO聯(lián)合走時(shí)偏振角度信息的算法進(jìn)行對(duì)比,驗(yàn)證本算法的可行性。

將震源搜索范圍設(shè)置為x[-50 m,50 m]、y[-50 m,50 m]、z[-50 m,0 m],同時(shí)在實(shí)驗(yàn)現(xiàn)場(chǎng)布設(shè)12個(gè)三軸振動(dòng)傳感器,如表1所示。

圖3 實(shí)驗(yàn)現(xiàn)場(chǎng)圖Fig.3 Experimental site map

表1 傳感器布設(shè)位置Tab.1 Sensor placement location

在區(qū)域中預(yù)設(shè)三組震源坐標(biāo)分別為(10,40,-20)、(40,-20,35)、(-30,20,-45),預(yù)設(shè)震源位置和傳感器在三維空間中的布設(shè)如圖4所示。按照預(yù)設(shè)震源依次進(jìn)行三發(fā)爆炸實(shí)驗(yàn),圖5為第一組實(shí)驗(yàn)中傳感器12接收到的信號(hào)。

圖4 預(yù)設(shè)震源位置和傳感器布設(shè)圖Fig.4 Preset hypocenter location and sensor layout

圖5 部分傳感器輸出信號(hào)圖Fig.5 Partial sensor output signal diagram

3.2 實(shí)驗(yàn)結(jié)果及評(píng)價(jià)方法

3.2.1實(shí)驗(yàn)結(jié)果分析

按照式(9)計(jì)算整個(gè)區(qū)域的能量,并繪制在預(yù)設(shè)震源處的三維能量場(chǎng)切片圖,如圖6所示。圖像可以直觀地反映出空間中能量的分布狀況,同時(shí)也可以在一定程度上反映震源的大致位置,但無法準(zhǔn)確得到震源的具體坐標(biāo)信息。

圖6 三維能量場(chǎng)切片圖Fig.6 3D energy field slice map

設(shè)置混沌映射程度α為4,混沌模型的迭代次數(shù)為300次,隨機(jī)生成120個(gè)粒子,作為文獻(xiàn)[13]算法迭代的初始粒子群,并將其分成2組各60個(gè)粒子,作為本算法的主、輔粒子群的初始粒子,將這些粒子按照Logistic模型分布在預(yù)設(shè)區(qū)域中,繪制粒子分布如圖7所示。通過對(duì)比可以看出經(jīng)過混沌搜索模型后粒子在空間中的分布相對(duì)于隨機(jī)均勻分布的方式,分布范圍更加廣泛。

設(shè)置兩種算法的迭代次數(shù)為500次,以預(yù)設(shè)震源1為例繪制兩種算法的粒子搜索軌跡,如圖8所示。通過圖8可以看出,本算法相比文獻(xiàn)[13]的算法,擴(kuò)大了搜索路徑的范圍,使得算法能搜索到更合適的震源位置。

圖7 初始粒子位置分布圖Fig.7 Initial particle position distribution map

圖8 算法軌跡對(duì)比圖Fig.8 Algorithm trajectory comparison chart

適應(yīng)度曲線對(duì)比圖如圖9所示。分析圖9,不難看出文獻(xiàn)[13]算法的適應(yīng)度曲線在迭代次數(shù)20次左右過早地進(jìn)入了收斂,30次左右群體陷入了局部收斂, 這導(dǎo)致了收斂區(qū)域錯(cuò)誤;而本算法則在10次左右開始收斂,45～55次左右再次收斂,然后在60次左右繼續(xù)收斂,100～140次左右后收斂速度逐漸下降,在迭代150次后最終達(dá)到收斂狀態(tài)。圖9直觀地反映出本算法具有擺脫局部收斂的優(yōu)勢(shì),這使得收斂區(qū)域更加準(zhǔn)確,從而提高了定位的精度。

圖9 適應(yīng)度曲線對(duì)比圖Fig.9 Fitness curve comparison chart

3.2.2算法評(píng)價(jià)方法

在通過DIR-QPSO聯(lián)合SRP對(duì)整個(gè)空間進(jìn)行精細(xì)化定位,得到最終的震源坐標(biāo)之后,對(duì)震源定位誤差進(jìn)行評(píng)價(jià)。本文采用均方誤差(MSE)來判斷震源定位結(jié)果的準(zhǔn)確性。針對(duì)預(yù)設(shè)的3組震源位置,每組重復(fù)測(cè)試算法10次取平均值,得到震源定位結(jié)果,同時(shí)分別測(cè)量?jī)煞N算法平均總收斂時(shí)間,如表2所示。

表2 震源定位結(jié)果Tab.2 Source location results

從表2中可以直觀地看出文獻(xiàn)[13]算法針對(duì)3組預(yù)設(shè)震源的定位結(jié)果,平均誤差均小于0.49 m,均方誤差MSE最小可達(dá)0.001 0 m,最大則有0.004 74 m;而DIR-QPSO聯(lián)合SRP的算法定位結(jié)果,平均誤差均小于0.235 m,MSE最小可達(dá)0.000 3 m,最大則只有0.001 32 m。同時(shí)DIR-QPSO聯(lián)合SRP算法的總收斂時(shí)間最大減少了0.539 2 s,最小也減少了0.106 7 s。通過對(duì)比可以得出:文獻(xiàn)[13]的算法與本文算法在實(shí)際定位結(jié)果雖然基本一致,但本文算法總體精度優(yōu)于文獻(xiàn)[13]的算法,誤差和總收斂時(shí)間也相對(duì)較少。

4 結(jié)論

本文提出了一種基于DIR-QPSO聯(lián)合SRP的地下爆炸震源定位方法。該方法通過Logistic混沌模型構(gòu)建初始化種群,利用DIR-QPSO相互約束的兩個(gè)主輔粒子群,消除了震源的搜索過程中容易出現(xiàn)的局部收斂的現(xiàn)象,同時(shí)結(jié)合SRP的能量聚焦特性,對(duì)震源位置進(jìn)行快速定位。

通過實(shí)驗(yàn)驗(yàn)證,及與文獻(xiàn)[13]中的QPSO聯(lián)合走時(shí)偏振角度信息定位的算法進(jìn)行對(duì)比,得出了本文算法具有定位精度較高、誤差小、搜索路徑范圍大、定位速度快和收斂范圍準(zhǔn)確的優(yōu)勢(shì),在地下淺層復(fù)雜空間定位研究領(lǐng)域具有一定的工程應(yīng)用價(jià)值。