董隴軍，胡清純，童小潔，劉有芳

a School of Resources and Safety Engineering, Central South University, Changsha 410083, China

b Department of Geophysics, Colorado School of Mines, Golden, CO 80401, USA

1. 引言

隨著淺部礦產(chǎn)資源的逐漸枯竭和地下空間的利用，地下工程逐漸向深部進發(fā)。在深部區(qū)域，巨大的地應力加上構造應力，使得巖體積蓄了大量的能量。地下開挖如同卸載，使得巖體積蓄的能量瞬間猛烈釋放，形成巖爆等災害。這種高應力引起的災害往往對巖體工程造成很大破壞，且常造成人員的傷亡。巖爆災害涉及諸多復雜因素的巖體力學行為過程，很難用傳統(tǒng)的巖石力學理論來進行分析和解釋，是目前的世界性難題之一。

微震監(jiān)測技術作為一種有效手段，可以對巖爆進行監(jiān)測預警[1]。目前，微震監(jiān)測技術已經(jīng)廣泛應用于深、高應力礦山的地壓安全監(jiān)測[2-6]。此外，加拿大、南非等國家也建立了全國性的礦山微震監(jiān)測網(wǎng)絡系統(tǒng)[7,8]。并且，微震監(jiān)測技術還在航空、軍工、橋梁結構、隧道工程、精密制造等各個領域得到了迅速發(fā)展和應用 [9-13]。

震源定位是微震/聲發(fā)射監(jiān)測中最經(jīng)典和最基本的問題之一[13]。1912年，Geiger [14]首次提出了定位方法。此后，提高微震震源定位的精度和效率，一直是學者們在微震監(jiān)測技術領域的重要研究內(nèi)容。隨著大量定位方法的提出，震源定位的誤差和所需時間不斷降低。經(jīng)典的源定位方法有Inglada方法[15]、USBM方法[16,17]、Thueber方法[18]、Powell方法[19]、單純型法[20]、雙重殘差法[21]等。近些年，Dong和Li等[22,23]研究了波速對定位的影響，提出了無需預先測速的定位方法。該方法避免了測量速度區(qū)域與實際微震區(qū)域的差異，大幅提高了定位精度。該方法還節(jié)約了前期測量速度造成的人員、時間和經(jīng)濟成本，比傳統(tǒng)方法更加方便實用。

Dong等在此基礎上發(fā)展出了三維綜合解析定位[24,25]、多步定位[26]、數(shù)值解析協(xié)同定位[27]等方法，取得了良好的效果。隨著研究的深入，一些學者開始關注物體的幾何形狀對定位的影響。Park等[9]、Baxter等[28]、Eaton等[29]通過訓練數(shù)據(jù)建立數(shù)據(jù)集合來評估震源位置。這類方法原理上可以適用于任何發(fā)生在物體表面震源的定位，但對于塊體結構內(nèi)部震源的定位，其適應性較弱。目前關于此類研究的試驗對象均為平面結構，且定位方法本身非常耗時。Gollob等[30]考慮了對象的幾何情況對定位的影響，但沒有考慮預測波速與實際波速之間的誤差。

本文針對含孔洞的三維結構，提出了一種基于A*搜索算法的無需測速震源定位方法（VFH）。該方法通過等距網(wǎng)格點來搜索路徑的方式避免了人工重復訓練，引入A*算法[31,32]并采用網(wǎng)格點來適應含有空區(qū)的復雜結構，并利用了無需預先測速的定位方法的優(yōu)點。該方法能有效地適應于三維復雜結構，并顯著提高了震源的定位精度。

2. 理論與方法

VFH使用改進的A*搜索算法和無需預先測速定位方法的思想來確定聲發(fā)射源的位置。該方法可分為四個部分：首先，劃分需要定位物體的網(wǎng)格并用0和1表示物體的形狀；其次，采集傳感器接收到的聲發(fā)射事件產(chǎn)生的P波到時信號；然后，使用A*算法搜索每個傳感器與每個網(wǎng)格點之間的最快波形路徑；最后，根據(jù)引入最小偏差量D來確定聲發(fā)射源的位置。圖1為VFH的定位流程圖。本節(jié)將描述該方法的每一部分。

圖1. VFH的定位過程流程圖。

2.1. 確定初始環(huán)境

在定位區(qū)域上，確定空區(qū)的幾何形狀和具體位置。根據(jù)空區(qū)的情況和定位的精度要求，確定單位立方體網(wǎng)格的尺寸大小。一般來說，網(wǎng)格劃得越密，定位精度越高，計算量也會成倍地增加。但當網(wǎng)格足夠密時，繼續(xù)增加網(wǎng)格密度，定位精度將不會有明顯變化。建立和網(wǎng)格節(jié)點相同尺寸的零矩陣M，將矩陣索引位置(x, y, z)與網(wǎng)格節(jié)點位置一一對應，并將對應空區(qū)的M (x, y, z)值更改為1。網(wǎng)格節(jié)點形成一個集合，在后續(xù)節(jié)點間搜索最快波形路徑時它們被作為起始點。假定P波在周圍非空區(qū)域的傳播速度為一個定值，為未知數(shù)，用C來表示。對于未知震源P0，設其位置坐標為(x0, y0, z0)，激發(fā)的初始時間為t0。

2.2. 搜索最快波形路徑

將集合內(nèi)的每個網(wǎng)格節(jié)點Pxyz當作潛在震源的激發(fā)位置，追蹤最短路徑，得到網(wǎng)格節(jié)點到第k個臺站的理論最短路徑。若網(wǎng)格節(jié)點Pxyz位于空區(qū)內(nèi)，則認為= ∞。VFH采用改進的A*算法來追蹤最短路徑

2.2.1. A*搜索算法

A*搜索算法由斯坦福研究所（現(xiàn)為SRI國際）的Peter Hart、Nils Nilsson和Bertram Raphael于1968年首次提出[31,32]。在計算機科學中，A*算法作為Dijkstra算法的擴展，因其高效性而被廣泛應用于搜尋最短路徑及圖的遍歷。A*算法通過使用啟發(fā)式的方法來指導搜索路徑，從而實現(xiàn)更好的性能。在其主循環(huán)的每次迭代中，A*算法根據(jù)當前節(jié)點到終點的估價代價f(m)來確定擴展路徑。具體來說，A*算法始終選擇最小長度的路徑。

式中，m表示選擇路徑的當前節(jié)點；d(m)表示從起始節(jié)點到節(jié)點m的實際最快路徑長度；而h(m)表示一個啟發(fā)式函數(shù)，用于估算從節(jié)點m到目標的最快路徑長度[33]。在本文中估值函數(shù)h(m)采用歐氏距離。當從起始點拓展到目標點或者沒有路徑可以拓展時，A*算法停止搜索。

在實際三維應用中，傳統(tǒng)的A*算法采用中心點，一般只考慮相鄰層的26個節(jié)點來選擇下一個節(jié)點，如圖 2（a）、（b）所示。在L形的區(qū)域內(nèi)，利用傳統(tǒng)的A*算法在內(nèi)追蹤最短路徑，搜索到如圖2（c）所示的一條路徑。從圖中可以發(fā)現(xiàn)，追蹤的最短路徑存在兩處不合理的地方：①搜索得到的路徑有明顯的鋸齒狀，這是傳統(tǒng)的A*算法的自身限制造成的；②路徑的節(jié)點均為立方體網(wǎng)格的中心，這意味著傳感器也要安放在立方體網(wǎng)格的中心，與實際不符合。

2.2.2. 改進的A*算法

為了更有效地追蹤最短路徑，我們改進了A*算法，采用網(wǎng)格點進行搜索，見圖3（a）。這樣可以使得搜索得到的路徑的節(jié)點均在網(wǎng)格節(jié)點上。這也意味著傳感器可以貼在物體的表面節(jié)點上，從而更加符合實際情況。

為了讓搜索得到的路徑不具有明顯的鋸齒狀，我們讓節(jié)點與周圍更多層的節(jié)點建立有效聯(lián)系。傳統(tǒng)的A*算法中，一個節(jié)點向相鄰一層的26個節(jié)點進行拓展。這意味著當前節(jié)點向周圍拓展的方向只有26個可以選擇。

讓節(jié)點與周圍更多相鄰層的節(jié)點之間進行聯(lián)系，可以使得當前節(jié)點搜索路徑時選擇的方向更多，這使得搜索得到的路徑更精確。根據(jù)節(jié)點拓展模型的對稱性，只畫出模型的八分之一來進行解釋說明。圖3（b）、（c）分別顯示了當前網(wǎng)格節(jié)點與周圍兩層（124個節(jié)點）、三層（342個節(jié)點）建立了聯(lián)系。節(jié)點z(i)與層數(shù)i之間的關系表示為：

在向外拓展的過程中，部分節(jié)點間形成的方向重復，故可以不用考慮（圖1）。去除這些方向對應的節(jié)點間連接可以減少計算量。隨著每個節(jié)點與周圍更多的網(wǎng)格節(jié)點建立聯(lián)系，搜索得到的路徑的誤差越小，但同時帶來計算量的增加。

圖2. 傳統(tǒng)的A*搜索算法。（a）當前節(jié)點連接到相關的26個節(jié)點；（b）圖2（a）的主要視圖；（c）傳統(tǒng)A*搜索算法搜索的路徑。

圖3. 當前網(wǎng)格節(jié)點與關聯(lián)網(wǎng)格節(jié)點間形成的方向。（a）當前節(jié)點連接到相鄰的一層；（b）當前節(jié)點連接到相鄰的兩層；（c）當前節(jié)點連接到相鄰的三層。

設計一個塊體模型和一個長條狀模型來探討層數(shù)i的合理取值。將模型劃分網(wǎng)格，如圖4所示。假設在Source點觸發(fā)震源，波從Source點到達K點（K=A,B,C,D,E,F,G,H,A′,B′,C′,D′,E′,F′,G′,H′），形成路徑LK。將波的實際最快路徑距離DR記錄到表1中。分別使用i層（i=1,2,3,4,5,6,…）的模型追蹤路徑LK。搜索得到的路徑距離DSi與DR之間的相對誤差E可以表示為：

將得到的路徑的相對誤差E記錄到表1中，并選出每個模型的最大路徑誤差Ei-max。算法的計算量O與拓展的節(jié)點數(shù)Z呈正相關關系。因而，可以用節(jié)點數(shù)Z來近似地表示計算量O，即Oi∝Z(i)。圖（5）顯示層數(shù)i分別與Ei-max和關聯(lián)節(jié)點數(shù)Z(i)之間的關系。

從圖5中可以看出，當拓展的層數(shù)達到3層時，理論路徑的相對誤差小于3%，能基本滿足定位要求，且計算量增速相對較小。因此，本文中的A*算法考慮當前節(jié)點與周圍3層的節(jié)點建立聯(lián)系。為了確定路徑的具體位置，A*（或Dijkstra）算法在確定最短路徑的長度后，需再次進行反向搜索。本文在搜索路徑時增加一個數(shù)組，在當前節(jié)點對應的位置上記錄前一個節(jié)點的坐標。這樣就避免了反向搜索，提高了運算效率。

2.3. 采集到時數(shù)據(jù)

圖4. 節(jié)點間的理想路徑。

在待測物體的不同位置分別裝上m個傳感器。各個傳感器均位于網(wǎng)格點上，其位置均為已知。對于三維模型，未知數(shù)有5個[P波的波速、聲發(fā)射源坐標(x0,y0,z0)、激發(fā)的初始時間t0]，因而m需為大于或等于5的整數(shù)。對于接收到信號的第k個傳感器Sk，記錄其位置坐標為(xk,yk,zk)，接收到聲發(fā)射P波信號的初至到時為。計算兩個不同的傳感器Sl和Sm的實際到時差，用Δ表示。

2.4. 震源定位

對于點Pxyz激發(fā)的震源，理論旅行時間等于震源到第k個傳感器之間的最短傳播路徑除以波速C。兩個不同的臺站l和m的到時差等于旅行時間之差Δ。根據(jù)Δ與Δ之差的平方，引入Dxyz來描述點Pxyz與未知聲發(fā)射源P0的偏離程度，表示為：

式中，當采樣點落入空區(qū)內(nèi)，則令Dxyz= ∞。

每個網(wǎng)格點都將得到對應的Dxyz值。Dxyz的值越大，表示點Pxyz與未知震源P0的偏離程度越大。因此，最小的Dxyz值對應的坐標便可認為是微震/AE定位的坐標。

3. 結果與討論

為了評價新的微震震源定位算法性能，對空心柱體砂漿結構件進行了斷鉛試驗。將10 cm×10 cm×10 cm的立方體試件中間挖去了尺寸為Φ6 cm×10 cm的圓柱體，如圖6（a）所示。為了更快地計算，將該試件劃分成25×25×25個相同尺寸的立方體小網(wǎng)格方塊，形成定位模型，如圖6（b）所示。

表1 真實路徑距離與追蹤路徑距離之間的相對誤差

圖5. 相對誤差和相關節(jié)點數(shù)與模型層數(shù)的關系。

將空區(qū)所在方塊標記為不通過，其他方塊標記為通過。將6個聲發(fā)射傳感器固定在試件上，在傳感器與試件之間加入耦合劑來獲得更好的聲耦合。傳感器的位置均位于網(wǎng)格節(jié)點上[圖6（c）]，坐標位置見表2。

圖6. VFH定位過程。（a）試件；（b）試件建模和網(wǎng)格劃分；（c）確定傳感器的坐標（單位：cm）；（d）使用改進的A*搜索算法進行的從節(jié)點到傳感器的路徑搜索（單位：cm）。

為方便計算，將模型坐標進行轉化，來使得節(jié)點矩陣M的索引與坐標位置一一對應。AE數(shù)據(jù)采集時采用40 dB的閾值和5 MHz的采樣率。在試件不同位置進行斷鉛試驗，每個位置進行兩次，記錄傳感器接收聲發(fā)射事件產(chǎn)生的P波到時，記錄在表3中。搜索每個潛在震源到傳感器的路徑LS，并計算距離，如圖6（d）所示。再根據(jù)傳感器接收到的到時，計算每個潛在震源的偏差值Dxyz。確定最小Dxyz值對應在試件上的位置坐標(x,y,z)，并將位置坐標進行轉化，結果見表3。同時，用傳統(tǒng)的Geiger法進行定位，定位結果記錄在表3中。

表2 傳感器布置的坐標

分析定位結果，將新方法和傳統(tǒng)方法的定位結果與實際斷鉛點進行對比，將誤差記錄在表3內(nèi)。從表中可以看出，傳統(tǒng)的Geiger法的最大定位誤差為4.5 cm，遠大于VFH的定位誤差。為了便于觀察，圖7（a）為兩種方法的定位結果和定位誤差的可視化圖形。圓圈尺寸代表震源定位結果的誤差大小?？梢悦黠@地看出傳統(tǒng)定位方法的圓圈相較于VFH的要大很多。圖7（b）為兩種方法定位誤差的箱線圖。在圖7（b）中，使用VFH得到的震源定位誤差的中位數(shù)約為1.0 cm，而使用傳統(tǒng)方法得到的震源定位誤差的中位數(shù)約為1.9 cm。根據(jù)表中每個斷鉛試驗的定位誤差，可以容易地計算出新方法的平均值為1.20 cm，而傳統(tǒng)方法的定位誤差為2.02 cm。新方法的平均定位精度較傳統(tǒng)方法提高了近40%。因此，在復雜的三維結構中，新方法的定位精度較傳統(tǒng)方法有了較大的提高。

表3 斷鉛點的定位結果和誤差

（續(xù)表）

圖7. 兩種方法得到的震源定位結果及誤差。（a）定位結果及誤差可視化；（b）定位誤差箱線圖。

4. 結論

在本文中，我們提出了一種基于A*算法的無需測速定位方法（VFH）。我們通過等距網(wǎng)格點來搜索路徑的方式避免了人工重復訓練。對于含有不規(guī)則空區(qū)的三維復雜結構，我們改進了A*算法并采用網(wǎng)格點來搜索符合彈性波實際傳播規(guī)律的路徑。并且，我們還在計算中將彈性波的波速當作未知數(shù)，來減小因監(jiān)測過程中物體波速的變化帶來的定位誤差。

為了評價新方法的準確性和有效性，在空心立方體混凝土結構件上進行了斷鉛試驗。定位結果顯示，新方法得到的平均定位誤差為1.20 cm，它的平均定位精度較傳統(tǒng)的Geiger法提高了40%。這表明新方法能有效地適應幾何不規(guī)則的三維復雜結構。

該方法克服了傳統(tǒng)定位方法中波速難以確定的缺陷，考慮了幾何形狀對定位的影響，且適用于三維定位，完善了微震定位方法。本文的方法亦適用于涉及聲發(fā)射定位的無損檢測等其他領域。

致謝

感謝國家自然科學基金（51822407、51774327）、湖南省自然科學基金（2018JJ1037）、中南大學創(chuàng)新驅動計劃項目（2020CX014）的資助。

Compliance with ethics guidelines

Longjun Dong, Qingchun Hu, Xiaojie Tong, and Youfang Liu declare that they have no conflict of interest or financial conflicts to disclose.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2019.12.016.