土體中巖石破壞次聲波的三維多測點振速矢量直線匯聚聲源定位方法

趙久彬， 劉元雪， 柏 準， 楊駿堂， 何少其

(1. 陸軍勤務(wù)學(xué)院 軍事設(shè)施系，重慶 401311； 2. 巖土力學(xué)與地質(zhì)環(huán)境保護重慶市重點實驗室，重慶 401311)

滑坡泥石流等地質(zhì)災(zāi)害發(fā)生前，巖石的爆裂、摩擦和斷裂等破壞會產(chǎn)生低頻的次聲波[1-2]，由于次聲波長較長頻率低的特點，使其能夠在傳播中衰減少、抗干擾能力和穿透能力較強，這些優(yōu)勢使得次聲波能夠遠距離傳輸，可以不考慮次聲傳感器和被測物體的耦合性。次聲波的監(jiān)測技術(shù)在地質(zhì)災(zāi)害領(lǐng)域的研究,主要采用聲發(fā)射技術(shù)、小波變換、小波包變換、時頻分析等信號處理方法對次聲信號的相對能量、頻率等頻譜和時頻特征進行分析，獲取出地質(zhì)災(zāi)害發(fā)生前次聲信號的重要特征，為地質(zhì)災(zāi)害監(jiān)測預(yù)警前兆特征所參考[3]。由于次聲信號傳播距離遠衰減少，用在地質(zhì)災(zāi)害發(fā)生地點的定位預(yù)測顯得格外有應(yīng)用價值，但國內(nèi)外對次聲波定位的研究還很少。

由于軍事應(yīng)用和海洋探測等領(lǐng)域的發(fā)展需求，聲源定位在海洋水聲定位的研究比較深入，聲源定位技術(shù)從傳統(tǒng)的時延算法[4-5]、波束形成算法[6-7]，向更加精確的匹配場[8-9]算法發(fā)展。但是在陸地的微震、泥石流、滑坡等聲波定位技術(shù)研究領(lǐng)域，還停留在室內(nèi)試驗階段[10]，黃曉紅等[11]采用在室內(nèi)巖土試塊內(nèi)激發(fā)超聲信號，通過部署在表面的傳感器收集信號，采用時延算法進行聲發(fā)射定位。對于巖石破壞的次聲信號研究，在信號處理領(lǐng)域的定位方法和聲場傳播理論方面還有很大的研究空間。

本文首先分析滑坡泥石流等地質(zhì)災(zāi)害中，災(zāi)害發(fā)生前夕由于巖石破壞產(chǎn)生次聲波機理，通過數(shù)值模擬發(fā)現(xiàn)土體中爆炸發(fā)出的次聲波的頻率范圍與巖石破壞的次聲波頻率范圍相當，可以通過該方法呈現(xiàn)出巖石破壞產(chǎn)生次聲波現(xiàn)象，使次聲波在土體及空氣中傳播。通過多個位于聲場中的次聲傳感器監(jiān)測點接收信號，基于測點振速矢量指向聲源的原理，提出一種三維多測點振速矢量直線匯聚次聲聲源定位方法，作為初步展示，采用數(shù)值模擬的方法建立土質(zhì)滑坡模型和土壤空氣分層介質(zhì)模型中發(fā)出次聲信號，在模型體中部署3個次聲傳感器收集信號，采用三維多測點振速矢量直線匯集次聲聲源定位方法，準確獲得了聲源坐標位置。

1 巖石破壞產(chǎn)生次聲波機理及數(shù)值模擬

1.1 巖石破壞產(chǎn)生次聲波機理

國內(nèi)外學(xué)者通過長期觀測和試驗發(fā)現(xiàn)，巖石結(jié)構(gòu)從裂紋產(chǎn)生、貫通和破壞過程中能夠輻射出聲波[12-17]。當巖石結(jié)構(gòu)處于非穩(wěn)定平衡狀態(tài)時，巖石內(nèi)部產(chǎn)生微觀下晶體之間的破裂和拓展，形成能量的積累。當能量積累到臨界值時巖石爆炸釋放能量產(chǎn)生應(yīng)力波，形成聲發(fā)射。巖石破壞過程中的聲發(fā)射會產(chǎn)生高頻和低頻聲波。高頻聲波由于頻率高衰減很快，無法進行遠距離傳播，監(jiān)測設(shè)備不能接收到信號。而低頻信號尤其是次聲波具有頻率低衰減少的特點，可以進行遠距離傳播。由于巖體的聲發(fā)射是滑坡崩塌等地質(zhì)災(zāi)害臨發(fā)的重要標志，可以對災(zāi)害體進行次聲波監(jiān)測，通過數(shù)學(xué)建模方法進行災(zāi)害發(fā)生地點進行定位，從而實現(xiàn)提前災(zāi)害預(yù)警[18]。

國內(nèi)外很多學(xué)者通過巖石破壞試驗收集次聲波特征，發(fā)現(xiàn)巖石破壞次聲波異常信號能量主要集中在3.19～7.81 Hz內(nèi)，而在巖土體爆炸與沖擊工程中，沖擊波引發(fā)和衰減形成次聲波在遠距離傳播，收集的次聲波信號的次聲能量集中在2～8 Hz。由于在地質(zhì)環(huán)境中實施土體巖石破壞試驗難度較大，可以在巖土體中設(shè)置小當量炸藥爆炸，模擬地質(zhì)災(zāi)害中巖石破壞產(chǎn)生的次聲現(xiàn)象，研究次聲定位方法。初期我們采用數(shù)值模擬方法，提出一種三維多測點振速矢量直線匯聚次聲聲源定位方法。

1.2 數(shù)值模型的建立

本文設(shè)計的模型一共有土體和炸藥兩部分。土體的尺寸為8 000 m×8 000 m×8 000 m，模擬巖石爆裂的小當量炸藥尺寸為邊長75 m×75 m×75 m的立方體，將其安放在土體中心，模擬地質(zhì)災(zāi)害發(fā)生前，巖石擠壓爆裂產(chǎn)生次聲波的現(xiàn)象[19]，如圖1所示。采用LS-DYNA程序中的多物質(zhì)ALE算法模擬巖石爆裂體爆炸后沖擊波傳播及土介質(zhì)的運動[20]。

圖1 立方體1/8網(wǎng)格模型

選用SOLID164實體進行網(wǎng)格劃分，由于模型是全對稱結(jié)構(gòu)，故可分析模型的1/8結(jié)構(gòu)，共有243 085個單元，250 047個節(jié)點。土體和巖石爆裂體采用映射方法劃分網(wǎng)格，其中土體采用非等距離劃分方法，距離中心近密遠疏，巖石爆裂體采用小劑量TNT炸藥材料替代。為了模擬出無限的土體區(qū)域，模型邊界均采用透射條件，選擇炸藥材料模型MAT_HIGH_EXPLOSIVE_BURN和EOS_JWL關(guān)鍵字為巖石爆裂體中炸藥的材料模型和狀態(tài)方程，土體選用土壤泡沫材料模型MAT_SOIL_AND_FOAM關(guān)鍵字。模型外界邊界執(zhí)行零位移全約束，數(shù)值模擬采用單位為m-kg-μs建模，計算時間為300 s,每3.6 s提取一次數(shù)據(jù)，計算生成K文件進行顯示動力運算。

1.3 標準計算結(jié)果

巖石爆裂體中爆炸發(fā)生后，產(chǎn)生巨大的能量釋放，形成球形沖擊波向外在土體中傳播，沖擊波的能量很大頻率很高，在土體介質(zhì)中衰減很快，衰減形成較低頻率的壓力波繼續(xù)向外傳播[21]。如圖2所示為爆炸過程中土體結(jié)構(gòu)在不同時刻的等值壓力圖:圖2(a)為發(fā)生爆炸時刻；圖2(b)說明爆炸后沖擊波以球面波形式向外擴大，球形沖擊波向外傳播形成更大的球面波，最外層能量最小，最里層能量最大，說明在傳播過程中發(fā)生能量衰減；如圖2(c)所示，在球面沖擊波向外傳播的同時，內(nèi)部出現(xiàn)了形狀猶如“水花”的壓力波，繼而新的沖擊波再次被激發(fā)，形成球形沖擊波線外傳播，但是第二輪沖擊波的能量較第一次低，形成了能量層次為“低-中-低”組成的球形波向外傳播；接著在中心再次產(chǎn)生了能量很低的壓力波，如圖2(d)所示，產(chǎn)生的沖擊波在中心形成漣漪，由于能量很低已經(jīng)無法向外傳播；最后如圖2(e)所示，之前形成的層次為“低-中-低”能量的球形波最終衰減為能量低的壓力波，可以傳播到很遠的距離而不衰減，這說明此種壓力波具有頻率低、衰減小、傳播遠的特征，說明產(chǎn)生了次聲波向外傳播；傳播最終達到模型邊界，由于邊界設(shè)置為透射條件，次聲波向外透射繼續(xù)傳播，如圖2(f)所示。

圖2 各時刻1/2面的等值壓力圖

1.4 各測點的信號壓力和振速矢量信號分析

為了了解土體中不同區(qū)域中受沖擊波影響的作用，設(shè)置3個不同位置的監(jiān)測點，其單元編號及坐標為A455(2 575.3,0,0)，A4024(1 588.42,0,1 782.23)，A13046(1 525.80,845.52,0)，如圖3所示為3個監(jiān)測點的單元位置,分別設(shè)置在爆裂點的正上方、左方較遠位置和右方較近位置，以便于比較各個不同特征位置接收到的應(yīng)力波信號異同。圖4為在爆裂期間3個監(jiān)測點接收到的應(yīng)力波變化曲線圖，從圖中我們可以發(fā)現(xiàn)，應(yīng)力波首先到達距離爆裂點最近的A13046監(jiān)測點，隨后分別達到距離較遠的A4024和A455監(jiān)測點，3個監(jiān)測點接收到的第一輪應(yīng)力波能量最大，隨后減小。

圖4 各監(jiān)測點壓力

圖5 A4024單元壓力與各向質(zhì)點振速對比圖

圖6 A4024單元壓力和質(zhì)點振速頻譜圖

2 次聲波傳播模型及定位方法

2.1 次聲波的多徑傳播模型

我們建立二維聲場模型，其中土層尺寸為8 000 m×4 000 m，土壤的參數(shù)為：ρ=1 500 kg/m3,cp=1 500 m/s。采用聲場傳播仿真軟件BELLHOP[22]，在MATLAB軟件進行在次聲波多徑傳播仿真，如圖7所示。由文獻[23]知土壤的剪切波傳播速度為cs≤200 m/s，遠遠小于縱波速度，且在土質(zhì)中衰減較快，當縱波到達時剪切波還未到達，故能在聲波信號中明確辨別，所以本文不考慮土層的剪切波作用。次聲源向接收點方向發(fā)出多條路徑的次聲波，其中直線傳播路徑為直達波，聲能損耗最小，其余路徑與上界面和下界面均發(fā)生了發(fā)射，其聲能將發(fā)生較大損耗。從圖中可以看到直達波與接收點的路徑最短，故直達波最先達到接收點，由于直達波聲能在所有聲線中最大，且最早到達接收點，所以在接收傳感器中的信號序列中，第一個聲壓波峰即為直達波的聲壓。

圖7 次聲波多徑傳播模型

根據(jù)文獻[24]，當在彈性介質(zhì)中不考慮剪切波時，其點聲源存在速度勢函數(shù)為

(1)

質(zhì)點的各向振速為

(2)

由于各質(zhì)點振速分量為速度勢函數(shù)的導(dǎo)數(shù)，與時間因子無關(guān)，所以各質(zhì)點振速分量應(yīng)同時達到最值。由于直達波聲線傳播路徑最短，傳播時間最小，且聲能損失最小，在接收點振速信號中第一個最大波峰即為直達波信號，本文采用該信號進行多測點振速矢量直線匯聚定位計算。

2.2 三維多測點振速矢量直線匯聚次聲聲源定位方法

巖石破壞產(chǎn)生的沖擊波可視為一種應(yīng)力彈性波，向外傳播后衰減形成次聲波向外繼續(xù)傳播[25]。由于應(yīng)力的方向與振速的方向具有一致性，在介質(zhì)中次聲波到達某質(zhì)點時，打破了該質(zhì)點的靜力平衡，其增加的最大的應(yīng)力方向即為振速最大的方向，也是縱波的傳播方向。本文提出一種三維多測點幾何定位聲源方法，以三測點為例。

在聲場中布置3個不共面不共線測點，其位置分別為P1(x1,y1,z1)，P2(x2,y2,z2)，P3(x3,y3,z3)，為了測得次聲波到達此點在各個方向引起的質(zhì)點運動情況，測得3個測點的最大振速為V1(p1,q1,r1)，V2(p2,q2,r2)，V3(p3,q3,r3)，故其振速的方向矢量分別為

(3)

設(shè)L1，L2，L3分別為經(jīng)過測點且其法線為方向矢量的直線方程，其形式為

(4)

由于所考慮聲波為直線傳播，在不考慮聲波的折射和反射等影響下，上述三條振速矢量直線方程理論上的交點則為聲源所在點，但實際中存在測量誤差、折射等影響因素，往往三條振速矢量直線不能交于一點，但一定會向一個區(qū)域內(nèi)聚集。如圖8所示，線段AB垂直于直線L1與直線L2，也是L1與L2兩直線間最短距離，若由直線L1與直線L2確定聲源預(yù)測點，則線段AB中點O′12(x12,y12,z12)為最優(yōu)預(yù)測點，證明如下：

圖8 三測點振速矢量直線匯聚于聲源點

由于聲源在線段AB上任意點出現(xiàn)的概率相等，那么需要在AB線段上尋找一點Q(xQ,yQ,zQ)，該點與其他點的距離總和最小時，則出現(xiàn)聲源的概率最大。采用函數(shù)最值法求解，設(shè)AB長度為L，取目標函數(shù)數(shù)D(xQ,yQ,zQ)，為Q點與其他點距離平方的總和，則其表達式為

(5)

對其偏微分得

(6)

令式(6)中方程式為0，得到

(7)

式(7)表明，當線段AB為均值密度時，Q(xQ,yQ,zQ)的坐標為線段AB的質(zhì)心，也就是AB的中點，那么對于由直線L1與直線L2確定的最優(yōu)聲源預(yù)測點為線段AB的中點，該點坐標為O′12(x12,y12,z12)。

同理，線段CD垂直于直線L2與直線L3，確定的最有聲源預(yù)測點坐標為O′23(x23,y23,z23)；線段EF垂直于直線L1與直線L3，確定的最優(yōu)聲源預(yù)測點坐標為O′13(x13,y13,z13)。所以，上述的聚集區(qū)域就為由O′12，O′23，O′13組成的三角形，真實聲源將出現(xiàn)在該三角區(qū)域中，按照上述函數(shù)最值法同樣可以證明，三角形的質(zhì)心O′與其他點的距離總和最小，故將質(zhì)心坐標作為聲源預(yù)測點的概率最大。則三測點振速矢量直線的聲源點定位點坐標為

(8)

波導(dǎo)環(huán)境中的密度變化會影響該三角區(qū)域的大小，例如：當接收器位于聲速較小，次聲源位于聲速較大的波導(dǎo)環(huán)境中，其定位的三角區(qū)域會增大；當接收器位于聲速較大，次聲源位于聲速較小的波導(dǎo)環(huán)境中，其定位的三角區(qū)域會減小。由于本文的定位算法基于概率和統(tǒng)計優(yōu)化的數(shù)學(xué)原理，可以有效減少由于界面密度變化產(chǎn)生的定位誤差，在應(yīng)對不同介質(zhì)之間傳播定位具有適用性。

以上為本文提出的三維多測點多振速矢量直線匯聚次聲聲源定位方法，下面以1.2節(jié)建立的模型具體說明。選取該模型的單元編號為A13046，A4024，A455三點進行驗證該算法。如圖9～圖11所示為各個質(zhì)點的振速變化圖，可知A455監(jiān)測點的振速最大值V1(0.001 174 62,06.35×10-6,6.37×10-6)m/s，A4024監(jiān)測點的振速最大值為V2(0.000 803 095,0.000 009 69,0.000 900 883)m/s，A13046監(jiān)測點的振速最大值為V3(0.001 488 24，0.000 836 661，0.000 018 1)m/s。根據(jù)式(4)求出三條經(jīng)過測點且指向聲源方向的直線方程如圖12所示，其形式為

圖9 A455監(jiān)測點各向振速

圖10 A4024監(jiān)測點各向振速

圖11 A13046監(jiān)測點各向振速

圖12 模型中三測點直線匯聚示意圖

(9)

根據(jù)上述方法，求出各坐標點如表1所示。其中：O點為爆裂點，為所建模型的原點；O′為本文提供的聲源定位方法所得的聲源坐標，發(fā)現(xiàn)該坐標位置在爆裂模型內(nèi)，結(jié)果表明該方法預(yù)測出的聲源位置位于巖爆體內(nèi)，且圓概率誤差僅有8.09 m，這個誤差在地質(zhì)災(zāi)害監(jiān)測預(yù)警工程中，是可以接受的。

表1 三測點定位坐標

3 數(shù)值模擬驗證

3.1 土質(zhì)滑坡

為了驗證該定位方法的有效性，本文建立實體土質(zhì)滑坡爆炸模型。某一處土質(zhì)滑坡，其尺寸如圖13所示，其爆裂體中炸藥尺寸為6 m×6 m×6 m，采用LS-DYNA程序中的多物質(zhì)ALE算法模擬爆裂后次聲波的傳播運動。其數(shù)值模擬方法和參數(shù)與1.2節(jié)方法一致。

圖13 土質(zhì)滑坡1/2網(wǎng)格模型

按照上述方法，為便于施工方便，在滑坡體的表面任意選擇3個地點安裝次聲波傳感器，收集當?shù)氐母飨蛸|(zhì)點振速信號。為便于敘述，本文選擇了如下地點如圖14所示，其單元編號及坐標為H15028(218.30,40.85,101.77)m，H20706(38.06,130.97,9.52)m，H16918(115.93,92.03,203.42)m，以便于比較各個不同特征位置接收到的應(yīng)力波信號異同。

圖14 土質(zhì)滑坡次聲傳感器所在位置

分析得到H15028監(jiān)測點的振速最大值V1(0.000 143 017, 2.93×10-5, 6.25×10-5)m/s，H20706監(jiān)測點的振速最大值為V2(0.000 123 284, 0.000 485 8, 3.07×10-5)m/s，H16918監(jiān)測點的振速最大值為V3(4.76×10-5, 3.86×10-5, 8.54×10-5)m/s。求出三條經(jīng)過測點且指向聲源方向的直線方程如下，其空間位置如圖15所示。

圖15 滑坡模型中三測點直線匯聚示意圖

(10)

根據(jù)本文提出的方法，求出各坐標點如表2所示。其中：R點為爆裂點，為所建模型的原點；R′為本文提供的聲源定位方法所得的聲源坐標，發(fā)現(xiàn)該坐標位置在爆裂模型內(nèi)，結(jié)果表明該方法預(yù)測出的聲源位置位于巖爆體內(nèi)，且圓概率誤差僅有5.11 m，這個誤差在地質(zhì)災(zāi)害監(jiān)測預(yù)警工程中，說明該算法能夠精確定位，應(yīng)用于災(zāi)害監(jiān)測預(yù)警具有廣闊前景。

表2 滑坡模型中三測點定位坐標

3.2 土壤和空氣不同介質(zhì)

聲波從介質(zhì)向另一種介質(zhì)傳播將發(fā)生反射和透射現(xiàn)象，其規(guī)律滿足斯涅爾定律。實際的巖土體的局部破壞發(fā)生在巖土層中，次聲波將從巖土層中透射到空氣中。為探究聲源定位方法在在不同介質(zhì)中應(yīng)用的定位效果和誤差，我們建立了土壤-空氣分層模型。其尺寸如圖16所示，其中爆裂體中炸藥尺寸為6 m×6 m×6 m，采用LS-DYNA程序中的多物質(zhì)ALE算法模擬爆裂后次聲波的傳播運動。其數(shù)值模擬方法和參數(shù)與土質(zhì)滑坡模型一致。

圖16 土壤和空氣的1/4模型

選擇測點如圖17所示，其單元編號及坐標為A958549(3,93,18)m，A961489(3,57,90)m，A529488(51,54,3)m，以便于比較各個不同特征位置接收到的應(yīng)力波信號異同。各監(jiān)測位置收集到傳感器的各向振速數(shù)據(jù)，分析得到A958549監(jiān)測點的振速最大值V1(5.47×10-17, 1.51×10-15, 3.7×10-16)m/s，A962351監(jiān)測點的振速最大值為V2(3.808 76×10-9,7.917 39×10-8,1.326 95×10-7)m/s，A529488監(jiān)測點的振速最大值為V3(0.000 049 8,0.000 056 7,0.000 003 53)m/s。求出三條經(jīng)過測點且指向聲源方向的直線方程如下，其空間位置如圖18所示。

圖18 不同介質(zhì)模型中三測點直線匯聚示意圖

(11)

根據(jù)本文提出的方法，求出各坐標點如表3所示，其中：P點為爆裂點，為所建模型的原點；P′為本文提供的聲源定位方法所得的聲源坐標，結(jié)果表明該方法預(yù)測出的聲源位置位于巖爆體內(nèi)，且圓概率誤差僅有4.04 m。為了比較本文定位方法在不同介質(zhì)中的定位效果，我們以最遠測點的距離為基準，比較土質(zhì)滑坡模型與不同介質(zhì)模型定位誤差。土質(zhì)滑坡的定位誤差為2.03%，而不同介質(zhì)模型定位誤差為4.26%，說明次聲波的界面?zhèn)鞑バ?yīng)會影響定位精度，不同介質(zhì)中的定位效果不如同一介質(zhì)定位效果。

表3 不同介質(zhì)模型中三測點定位坐標

在定位土壤中的次聲源，而接收傳感器位于空氣中時，由于次聲波在不同介質(zhì)之間傳播，其傳播路線將在界面產(chǎn)生偏折，導(dǎo)致定位誤差。由于折射角越大，定位誤差越小，故應(yīng)當將接收傳感器盡量放置于距離界面垂直距離較大處。另外，由于本文的多測點振速矢量直線匯聚聲源定位方法是一種數(shù)理統(tǒng)計優(yōu)化算法，可有效的減少由于傳播偏折產(chǎn)生的定位誤差，表明該方法應(yīng)用于工程快速定位具有很好的優(yōu)勢。

3 結(jié) 論

本文提出一種三維多測點振速矢量直線匯聚次聲源定位方法，并通過數(shù)值模擬驗證該算法有效性，得出如下結(jié)論：

(1) 本文基于聲場中質(zhì)點受聲波作用各向振速不同的原理，提出經(jīng)過測點位置且法線為振速方向矢量多條直線指向并匯聚于聲源發(fā)生處的假設(shè)，通過求多點質(zhì)心預(yù)測為聲源位置，該方法創(chuàng)新性強，原理簡單，應(yīng)用巧妙。

(2) 本文設(shè)計了數(shù)值計算實體滑坡模型和土壤空氣分層介質(zhì)模型中進行定位計算，通過三測點的實例驗證該方法的有效性，結(jié)果表明該方法預(yù)測出的聲源位置位于巖爆體內(nèi)，且定位誤差在4.26%以內(nèi)。

綜上，本文提出方法用于滑坡等泥石流監(jiān)測預(yù)警技術(shù)中，能夠提前準確預(yù)測滑災(zāi)害發(fā)生地點位置，在災(zāi)害預(yù)警技術(shù)中具有廣闊前景。