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

(1.中北大學(xué)省部共建動態(tài)測試技術(shù)國家重點實驗室,山西 太原 030051;2.中北大學(xué)信息探測與處理山西省重點實驗室,山西 太原 030051)

0 引言

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

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

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

1 基于SPR的震源定位模型

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

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

一個傳感器陣列有M個陣元,第m個傳感器接收到的信號可以表示為

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

(1)

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

對于空間中某一點q,該陣列的波束形成的定義為

(2)

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

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

(4)

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

又因為τml(q)=τm(q)-τl(q),表示第m個傳感器與第l個傳感器接收位置q處信號的時間差,故式(4)可表示為

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

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

(8)

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

(9)

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

2 震源搜索算法

2.1 QPSO算法

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

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

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

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

(11)

式(11)中,φ1和φ2為(0,1)內(nèi)的隨機數(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)的隨機數(shù);Xi,j(t+1)為第t+1代粒子群的位置;β為收縮擴張系數(shù),一般在(0,1)內(nèi)取值;mbest為種群平均最好位置,計算公式為

(13)

2.2 DIR-QPSO算法

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

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

(14)

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

輔粒子群的局部吸引點qi,j是以主粒子群局部吸引點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)的隨機數(shù)。

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

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

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

(17)

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

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

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

故引入混沌模型為Logistic模型搜索的方式,如式(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ù)評價當(dāng)前粒子位置是否是實際震源位置的重要指標(biāo)。

2.5 定位模型解算流程

1) 初始化種群

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

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

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

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

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

4) 重復(fù)迭代

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

5) 結(jié)果輸出

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

3 實驗驗證

3.1 實驗方案

實驗采用本單位研發(fā)的分布式地下震源定位系統(tǒng)進行定位,實驗區(qū)域為自然土介質(zhì),實驗現(xiàn)場如圖3所示。通過與文獻[13]中提出的QPSO聯(lián)合走時偏振角度信息的算法進行對比,驗證本算法的可行性。

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

圖3 實驗現(xiàn)場圖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è)震源依次進行三發(fā)爆炸實驗,圖5為第一組實驗中傳感器12接收到的信號。

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

圖5 部分傳感器輸出信號圖Fig.5 Partial sensor output signal diagram

3.2 實驗結(jié)果及評價方法

3.2.1實驗結(jié)果分析

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

圖6 三維能量場切片圖Fig.6 3D energy field slice map

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

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

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

圖8 算法軌跡對比圖Fig.8 Algorithm trajectory comparison chart

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

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

3.2.2算法評價方法

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

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

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

4 結(jié)論

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

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