地下震源能量場快速掃描定位算法

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

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

0 引言

地下淺層震源定位技術是實現地下淺層爆炸、爆破、礦產資源勘探的關鍵技術手段,這種技術也經常應用在地震災害評估以及應急救援響應的過程中[1-2]。與地下深層的震源定位相比,淺層空間因為地質結構復雜、速度信息未知等原因,造成了在大區域淺層震源搜索過程中的難點和熱點問題。

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

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

1 基于SPR的震源定位模型

波束合成的震源定位算法是對于傳感器接收到的信號進行時延補償-求和,得到該點的輸出功率,通過波束合成得到的輸出功率又稱為可控響應功率(SRP, steered-response power),對區域內的每個點進行波束合成計算,輸出功率最大的點即震源所在的位置。波束形成示意圖如圖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個傳感器的沖擊響應,“*”表示卷積,wm(t)為第m個傳感器接收到的噪聲,m=1,2,…,M。

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

(2)

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

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

(4)

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

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

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

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

(8)

去掉其中的自相關值,可表示為

(9)

震源q所在位置為區域中可控響應功率最大值處,即

2 震源搜索算法

2.1 QPSO算法

量子粒子群算法在粒子群優化算法(PSO, particle swarm optimization)的基礎上取消了粒子的移動方向屬性,這樣在增加了粒子位置隨機性的同時消除了與粒子之前運動的關系。因為在QPSO算法中粒子的位置和速度在量子空間中不能一起確定,所以用波函數代替粒子位置和速度,通過蒙特卡羅方法求解波函數,從而得出粒子位置和速度,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)內的隨機數。

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

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

(12)

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

(13)

2.2 DIR-QPSO算法

利用DIR-QPSO算法進行迭代時,首先在給定的區域中初始化一個粒子種群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)內的隨機數。

輔粒子群的局部吸引點qi,j是以主粒子群局部吸引點pi,j和輔粒子群全局最優位置gbestj決定的,其具體表現形式為

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

(15)

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

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

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

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

(17)

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

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

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

故引入混沌模型為Logistic模型搜索的方式,如式(18)所示:

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

(18)

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

2.4 構建基于SRP的適應度函數

通過式(9)和式(10)構建的適應度函數為

將式(20)作為DIR-QPSO算法的適應度函數,適應度最大的粒子位置就是當前種群的最優位置gbestj,通過適應度函數評價當前粒子位置是否是實際震源位置的重要指標。

2.5 定位模型解算流程

1) 初始化種群

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

2) 更新兩個種群的最優震源位置

分別根據式(20)構建的適應度函數計算得到震源群中適應度值最大的位置,設為當前震源種群最優震源gbest1和gbest2。

3) 更新全局最優位置

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

4) 重復迭代

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

5) 結果輸出

當達到迭代次數時,將全局震源最優位置Gbest輸出,并設置為最優化的震源位置。

3 實驗驗證

3.1 實驗方案

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

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

圖3 實驗現場圖Fig.3 Experimental site map

表1 傳感器布設位置Tab.1 Sensor placement location

在區域中預設三組震源坐標分別為(10,40,-20)、(40,-20,35)、(-30,20,-45),預設震源位置和傳感器在三維空間中的布設如圖4所示。按照預設震源依次進行三發爆炸實驗,圖5為第一組實驗中傳感器12接收到的信號。

圖4 預設震源位置和傳感器布設圖Fig.4 Preset hypocenter location and sensor layout

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

3.2 實驗結果及評價方法

3.2.1實驗結果分析

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

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

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

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

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

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

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

圖9 適應度曲線對比圖Fig.9 Fitness curve comparison chart

3.2.2算法評價方法

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

表2 震源定位結果Tab.2 Source location results

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

4 結論

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

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