薛 強(qiáng), 蔡承才, 王碩禾

(石家莊鐵道大學(xué) 電氣與電子工程學(xué)院，河北 石家莊 050043)

?

小波與相關(guān)函數(shù)算法在礦山震源定位中的應(yīng)用

薛 強(qiáng), 蔡承才, 王碩禾

(石家莊鐵道大學(xué) 電氣與電子工程學(xué)院，河北 石家莊 050043)

針對(duì)礦震震相提取分析中震相信號(hào)提取困難、噪聲影響大、精度低等問(wèn)題，提出從小波時(shí)頻分析和統(tǒng)計(jì)學(xué)的角度，根據(jù)小波域內(nèi)礦震信號(hào)與噪聲的不同特性以及不同采集通道間信號(hào)具有相關(guān)性特點(diǎn)，采用小波閾值濾波算法和相關(guān)函數(shù)算法對(duì)礦震信號(hào)進(jìn)行計(jì)算，得到了各個(gè)采集通道間較準(zhǔn)確的波到時(shí)差。試驗(yàn)結(jié)果表明，該算法對(duì)高斯噪聲和其他非相干噪聲具有抑制作用，能夠有效提高計(jì)算結(jié)果的精度。

數(shù)字信號(hào)處理；小波濾波；相關(guān)函數(shù)；震相提??；LabVIEW

0 引言

目前，我國(guó)許多礦體存在“一礦多開(kāi)”現(xiàn)象，越界開(kāi)采亦成為與之伴隨而來(lái)的一種違法行為，由于其隱蔽性，對(duì)其進(jìn)行地面監(jiān)測(cè)仍具困難。利用礦體炮采時(shí)所激發(fā)的地震波來(lái)計(jì)算作業(yè)點(diǎn)的時(shí)空參數(shù)是實(shí)現(xiàn)采掘面地面監(jiān)測(cè)和定位的有效手段之一[1]。在數(shù)字信號(hào)處理方面，針對(duì)地震波的時(shí)變性以及頻率成分復(fù)雜等特點(diǎn)，利用小波變換具有時(shí)頻分析的優(yōu)點(diǎn)可以有效消除隨機(jī)干擾因素的影響，觀察到信號(hào)的突變點(diǎn)，得到更好的P波到時(shí)識(shí)別效果[2-3]；并且由于在不同地點(diǎn)所接收到的同一震源所激發(fā)的地震信號(hào)具有波形相似性[4]，因此可以利用互相關(guān)函數(shù)求最大值的方法來(lái)精確測(cè)量?jī)蓚€(gè)地震波信號(hào)的到時(shí)差[5]。本文將小波閾值濾波算法和相關(guān)算法同時(shí)應(yīng)用于礦山震源定位，利用LabVIEW為應(yīng)用程序軟件實(shí)現(xiàn)平臺(tái)[6]，實(shí)現(xiàn)了礦山震源定位系統(tǒng)數(shù)據(jù)處理與計(jì)算，實(shí)驗(yàn)結(jié)果滿足國(guó)土資源廳關(guān)于越界開(kāi)采的10 m定位精度要求。

1 礦山微震定位系統(tǒng)架構(gòu)

當(dāng)?shù)V山地下炮采時(shí)激發(fā)地震波傳到地表，由布設(shè)在地表不同地點(diǎn)的多個(gè)地震檢波器采集信號(hào)，并將其發(fā)送到計(jì)算機(jī)[7]，進(jìn)而通過(guò)分析計(jì)算得到地下炮采精確位置坐標(biāo)。系統(tǒng)硬件結(jié)構(gòu)如圖1所示。其中，地震檢波器負(fù)責(zé)采集礦震信號(hào)；采集裝置通過(guò)A/D將模擬信號(hào)轉(zhuǎn)化為計(jì)算機(jī)數(shù)字信號(hào)；計(jì)算機(jī)完成各種計(jì)算。

圖1 硬件結(jié)構(gòu)圖

2 礦震信號(hào)的特點(diǎn)及處理算法

2.1 礦震信號(hào)的頻率特性及處理方法

礦震信號(hào)頻率成份復(fù)雜，包含各種干擾，為時(shí)變信號(hào)。通過(guò)對(duì)礦震信號(hào)的頻譜分析可知，其主頻成份集中在10～60 Hz，屬低頻信號(hào)，并同時(shí)具有兩個(gè)以15 Hz和45 Hz為中心的主頻帶，分別對(duì)應(yīng)S波和P波。在震源定位中，礦震波到時(shí)(或通道間的到時(shí)差)精確拾取最為重要，干擾的存在會(huì)影響震動(dòng)信號(hào)到時(shí)拾取的結(jié)果，因此將各種噪聲信號(hào)濾除掉并盡可能多的保留有用信號(hào)是精確定位的首要任務(wù)。但是礦震信號(hào)具有以下特點(diǎn)：(1)有效信號(hào)的頻率成分較多，在10～60 Hz之間連續(xù)存在有用信號(hào)；(2)礦山炮采時(shí)所產(chǎn)生的地震波主要由縱波、橫波和面波構(gòu)成，不同的傳播速度波速也不一樣，傳感器采集到的信號(hào)隨時(shí)間其頻率也在變化，為時(shí)變信號(hào)；(3)受周?chē)肼暤纫蛩氐挠绊?，其信?hào)變得更為復(fù)雜。傳統(tǒng)的傅立葉變換只能夠粗略得到信號(hào)的頻率成分，對(duì)于其時(shí)間特性則不能夠顯示。而小波變換是很好的時(shí)頻信號(hào)分析工具，克服了傅立葉變換不能夠自調(diào)節(jié)的同時(shí)分析時(shí)頻特性的缺點(diǎn)，并且小波波形與一些震動(dòng)波形非常相似，具有較好的匹配性。利用小波變換作為礦震信號(hào)的分析與濾波工具將得到較為滿意的結(jié)果。

2.2 礦震信號(hào)處理流程

因不同地點(diǎn)檢波器采集的同一震源發(fā)出的震動(dòng)信號(hào)具有相關(guān)性，故通道間信號(hào)到時(shí)差可采用互相關(guān)運(yùn)算來(lái)獲取。結(jié)合礦震信號(hào)噪聲影響大、頻率成分復(fù)雜等特點(diǎn)，作互相關(guān)運(yùn)算前首先要對(duì)原始信號(hào)進(jìn)行濾波，信號(hào)處理流程如圖2所示。首先，對(duì)原始信號(hào)采用小波閾值濾波，盡可能精確的還原有效信號(hào)。其次，由于高、低頻信號(hào)非主頻成分且相關(guān)性較主頻成分差，因此在互相關(guān)運(yùn)算前進(jìn)行帶通二次濾波；最后，對(duì)以上兩個(gè)步驟處理后的結(jié)果進(jìn)行互相關(guān)運(yùn)算，得到通道間信號(hào)的準(zhǔn)確到時(shí)差。

圖2 礦震信號(hào)處理流程

2.3 礦震信號(hào)的小波分析與濾波

小波變換是時(shí)變信號(hào)的很好分析工具，小波帶通濾波和小波閾值濾波是兩種現(xiàn)代數(shù)字信號(hào)處理方式。本文就礦震信號(hào)，分別應(yīng)用這兩種濾波方式進(jìn)行了討論，對(duì)其結(jié)果進(jìn)行分析并比較其優(yōu)略。

2.3.1 小波基函數(shù)的選取

選取小波基應(yīng)具有一般原則，首先小波應(yīng)具有緊支撐性，對(duì)稱(chēng)性和正則性；其次，選擇的小波基函數(shù)盡量與被分析的波形相近，從而在小波分解的過(guò)程中更容易逼近。事實(shí)上，有很多小波基函數(shù)與地震波形很相近，可以直接把它們作為小波基函數(shù)進(jìn)行數(shù)據(jù)的分析和處理。經(jīng)實(shí)驗(yàn)比對(duì)，本系統(tǒng)采用的小波基函數(shù)為db4小波。

2.3.2 經(jīng)典小波帶通濾波討論

以db4小波為母函數(shù)，對(duì)礦震信號(hào)進(jìn)行6層小波分解，然后對(duì)各層分別進(jìn)行單層重構(gòu)。從單層重構(gòu)中可以得到信號(hào)的能量主要集中在d2、d3、d4和d5層，在低頻段和高頻段能量都很小，認(rèn)為礦震波信號(hào)主要由這4層構(gòu)成。將能量較小的a1、d1和d6層小波系數(shù)置零，對(duì)原始信號(hào)進(jìn)行重構(gòu)。

帶通濾波結(jié)果如圖3所示，其中圖3(a)為原始信號(hào)，圖3(b)為采用帶通濾波后的信號(hào)。從圖3中濾波結(jié)果來(lái)看，帶通濾波后信號(hào)具有一定程度失真，濾波效果不理想，原因有以下兩點(diǎn)：(1)被剔除的a1、d1和d6層存在信號(hào)的有效成分，對(duì)其簡(jiǎn)單的切除，會(huì)引起信號(hào)的失真；(2)雖然d2、d3、d4和d5層的信噪比較大，但是這些層的信號(hào)仍然存在一定的噪聲，而帶通濾波未對(duì)這些噪聲做任何處理。為了盡可能精確的還原有效信號(hào)，需要在各層均保留有效信號(hào)并剔除噪聲，利用小波閾值濾波可達(dá)到該目的。

圖3 帶通濾波結(jié)果

2.3.3 小波閾值濾波及信號(hào)分析

(2)小波閾值濾波。本文采用軟閾值[10]的方法在小波域?qū)Ω鲗有〔ㄏ禂?shù)分別處理；然后對(duì)濾波后的各層系數(shù)進(jìn)行小波重構(gòu)，恢復(fù)原始信號(hào)，從而得到濾波后的信號(hào)。通過(guò)仿真分析可知，高頻部分的d1層和d2層濾波效果最為明顯，這是因?yàn)樵肼曄鄬?duì)于有效信號(hào)來(lái)說(shuō)，屬于高頻干擾，其成份多聚集在這兩層內(nèi)，因此這兩層的濾波效果也最好。圖4為小波閾值濾波結(jié)果前后對(duì)比，其中圖4(a)為原始信號(hào)，圖4(b)為濾波后的信號(hào)。從圖4中可以看到，信號(hào)上附帶的噪聲明顯得到了抑制，且避免了失真。

這種閾值濾波的優(yōu)點(diǎn)是在濾除各個(gè)尺度上噪聲的同時(shí)保留了有效信號(hào)。克服了前文所述帶通濾波存在的兩個(gè)缺點(diǎn)，使濾波后波形更接近于實(shí)際信號(hào)，提高了系統(tǒng)的抗干擾能力，這也正是閾值濾波優(yōu)于帶通濾波之處，因而可以減少由于礦震波到時(shí)誤判而引起的定位不準(zhǔn)確，為后續(xù)工作做好準(zhǔn)備。

圖4 小波閾值濾波結(jié)果

2.4 利用相關(guān)函數(shù)算法提取礦震信號(hào)的到時(shí)差

當(dāng)?shù)V體不太大，認(rèn)為礦震波傳播的介質(zhì)是均勻的且各向同性的，同一震源激發(fā)的礦震波在各個(gè)方向上具有相似性，而對(duì)于各通道噪聲而言，其相似程度很小，利用礦震信號(hào)的這一特性，求取兩通道間相關(guān)函數(shù)[9]，進(jìn)而得到兩通道間信號(hào)的到時(shí)差。因噪聲不具有相關(guān)性，因此這種算法也起到抑制噪聲的目的。

2.4.1 相關(guān)函數(shù)

設(shè)x(n)和y(n)為復(fù)數(shù)序列，定義

(1)

為x(n)和y(n)的互相關(guān)函數(shù)。其中，x*(n)為x(n)的共軛。rxy(m)是關(guān)于m的函數(shù)。 互相關(guān)和卷積在時(shí)域的關(guān)系為

(2)

相關(guān)技術(shù)是一種利用波形相似性來(lái)提取信號(hào)的一種技術(shù)，其特點(diǎn)是在信噪比較小的情況下，利用信號(hào)間具有很好的相關(guān)性而噪聲沒(méi)有相關(guān)性的特點(diǎn)，從強(qiáng)噪聲背景中有效的提取有效信號(hào)。這種方法對(duì)采集系統(tǒng)要求比較低，具有很強(qiáng)的抗干擾能力。

2.4.2 利用相關(guān)函數(shù)提取各通道信號(hào)間的到時(shí)差

礦震信號(hào)的低頻分量和高頻分量能量較小，且高頻信號(hào)各通道間的相關(guān)性變差。針對(duì)互相關(guān)運(yùn)算算法，切除該頻段信號(hào)相當(dāng)于做了二次濾波，參與運(yùn)算的量少了相關(guān)性較差信號(hào)的作用，結(jié)果精度反而會(huì)有所提高，操作分兩步：(1)帶通二次濾波。對(duì)軟閾值濾波后信號(hào)再剔除低頻和高頻成份，即去除能量較小、相關(guān)性較差的a1，d1和d6層，對(duì)剩余的主頻成份進(jìn)行重構(gòu)，恢復(fù)信號(hào)。此步驟可以結(jié)合軟閾值濾波一次性完成，為簡(jiǎn)少運(yùn)算量可先作帶通濾波，再作閾值濾波。而本文為了論述清晰，故分開(kāi)來(lái)討論。(2)互相關(guān)運(yùn)算。將恢復(fù)完成的各通道間信號(hào)兩兩進(jìn)行互相關(guān)運(yùn)算，得到相關(guān)函數(shù)并求取其最大值點(diǎn)坐標(biāo)，該最大值點(diǎn)坐標(biāo)即為兩通道間信號(hào)的到時(shí)差(數(shù)據(jù)以采樣點(diǎn)為單位，可以結(jié)合采樣頻率換算為時(shí)間)。

3 結(jié)果分析與定位系統(tǒng)實(shí)現(xiàn)

3.1 定位算法

本系統(tǒng)采用幾何定位算法，即對(duì)于每一個(gè)震源監(jiān)測(cè)點(diǎn)列寫(xiě)地震波在介質(zhì)中的傳播方程

(3)

式中，地面A(x1,y1,z1)坐標(biāo)為檢波器坐標(biāo)是已知量，震源位置P(x,y,z)是待求變量；t1為傳感器接收信號(hào)時(shí)刻為已知量；t為發(fā)震時(shí)刻亦為變量；v為信號(hào)傳播速度，其值由傳播介質(zhì)決定，亦為變量，故在地表至少布置5個(gè)以上傳感器，才能夠列寫(xiě)出5個(gè)以上的方程以求解5個(gè)變量。本系統(tǒng)首先應(yīng)用幾何定位方法進(jìn)行初步定位，求得足夠接近真實(shí)解的震源位置以及速度v。并以該解作為Geiger定位[10]的迭代初值進(jìn)行修訂。這樣既可以避免幾何定位的不準(zhǔn)確，又可以為Geiger定位提供足夠接近真實(shí)值的迭代初值。

3.2 定位試驗(yàn)

表1～表3為一組試驗(yàn)數(shù)據(jù)，震源的獲取采用人為地面放炮方式，地震波的傳播介質(zhì)為沙地，采用5通道檢波器測(cè)量信號(hào)，傳感器1～5布設(shè)按照近似正方形平面布設(shè)，其中心亦設(shè)一個(gè)傳感器，炮點(diǎn)和傳感器采用GPS進(jìn)行絕對(duì)坐標(biāo)標(biāo)定，傳感器包圍面積約為3 000 m2。采集卡采樣頻率設(shè)置為500 Hz。本次試驗(yàn)共收集到6組有效數(shù)據(jù)。其中炮(1～4)為同一放炮點(diǎn)，炮(5～6)為同一放炮點(diǎn)，傳感器坐標(biāo)固定不變。

將GPS絕對(duì)坐標(biāo)轉(zhuǎn)換為適合于式(3)的相對(duì)坐標(biāo)，傳感器和震源實(shí)際坐標(biāo)如表1所示，表中坐標(biāo)為換算后的相對(duì)坐標(biāo)，并以傳感器5坐標(biāo)點(diǎn)為參考坐標(biāo)點(diǎn)，定義該點(diǎn)坐標(biāo)為(0,0)。

表1 檢波器坐標(biāo)和震源坐標(biāo) m

以點(diǎn)5為參考點(diǎn)，定義其接收到信號(hào)時(shí)刻為0，以本文所述方法計(jì)算信號(hào)到時(shí)差，數(shù)據(jù)如表2所示。

表3為根據(jù)相關(guān)算法求到時(shí)差，然后根據(jù)幾何定位算法和Geiger算法計(jì)算出的結(jié)果對(duì)比，從表3數(shù)據(jù)可以看出，修正前計(jì)算值與真實(shí)值平均誤差為1.36 m，經(jīng)過(guò)修正以后其誤差平均值減小為0.98 m。定位結(jié)果滿足項(xiàng)目的定位精度要求。

本系統(tǒng)算法程序主要是在LabVIEW平臺(tái)下完成的，少數(shù)算法是通過(guò)LabVIEW的MATLAB Script節(jié)點(diǎn)，由MATLAB來(lái)完成復(fù)雜算法計(jì)算，該部分非本文所述重點(diǎn)，不再贅述。

表2 信號(hào)到時(shí)差 s

表3 相關(guān)算法定位結(jié)果 m

4 結(jié)論

各通道有效信號(hào)到時(shí)差的準(zhǔn)確拾取對(duì)于震源準(zhǔn)確定位具有重要意義。聯(lián)合小波濾波與相關(guān)函數(shù)算法能夠?qū)υ肼暭捌洚a(chǎn)生的影響進(jìn)行有效濾除或抑制。小波濾波能夠?yàn)V除各個(gè)尺度上的噪聲，同時(shí)在各個(gè)尺度上保留了有效信號(hào)；互相關(guān)運(yùn)算算法從宏觀或整體的角度可以找到兩個(gè)信號(hào)的到時(shí)差，并從統(tǒng)計(jì)學(xué)上的角度上對(duì)隨機(jī)高斯噪聲起到了進(jìn)一步抑制的作用，相當(dāng)于對(duì)信號(hào)的二次濾波，提高了算法的健壯性和抗干擾能力。實(shí)驗(yàn)表明，將兩種算法相結(jié)合應(yīng)用于礦山微震定位系統(tǒng)，其定位結(jié)果可以滿足系統(tǒng)的精度要求。

[1]馮國(guó)磊.煤礦盜采監(jiān)測(cè)系統(tǒng)的示范研究[D].長(zhǎng)春:吉林大學(xué),2006.

[2]劉曉紅.地震P波震相識(shí)別[D].北京:北京師范大學(xué),2008.

[3]徐宏斌,李庶林,陳際經(jīng).基于小波變換的大尺度巖體結(jié)構(gòu)微震監(jiān)測(cè)信號(hào)去噪方法研究[J].地震學(xué)報(bào),2012,34(1):85-96,127.

[4]魏富勝,許忠淮,郭祥云.基于波形相關(guān)性的地震事件分類(lèi)[J].地震地磁觀測(cè)與研究,2010,31(3):52-57.

[5]羅桂純.利用相關(guān)檢測(cè)法進(jìn)行地震波速及其變化的精確測(cè)量[D].北京:中國(guó)地震局地震預(yù)測(cè)研究所,2006.

[6]王斌.基于虛擬儀器的微震監(jiān)測(cè)系統(tǒng)研究[D].天津:天津大學(xué),2007.

[7]薛強(qiáng),葉明,范君菲.無(wú)線傳感器網(wǎng)絡(luò)在礦山微震定位系統(tǒng)中的應(yīng)用[J] .石家莊鐵道大學(xué)學(xué)報(bào)：自然科學(xué)版,2015,28(2):102-105,110.

[8]David L,Donoho.De-Noising by Soft-Thresholding[J]. IEEE Transactions on Information Theory, 1995, 41(3): 613-627

[9]胡廣書(shū).數(shù)字信號(hào)處理 [M].北京:清華大學(xué)出版社,2012.

[10]薛強(qiáng).震源定位系統(tǒng)在礦山越界開(kāi)采監(jiān)控中的應(yīng)用研究[D].石家莊:石家莊鐵道大學(xué),2011.

Application of Wavelet and Correlation Function Algorithm in Mining Source Localization

Xue Qiang， Cai Chengcai， Wang Shuohe

(School of Electrical and Electronic Engineering, Shijiazhuang Tiedao University, Shijiazhuang 050043, China)

According to the analysis of seismic phase extraction, the phase signal extraction is difficult, with big noise and low precision. This paper presents a new algorithm based on Wavelet time-frequency analysis and statistics. According to the different characteristics of wavelet domain of seismic signal and noise, and the correlation betweent different signal acquisition channels, using the wavelet threshold filtering algorithm and correlation algorithm to calculate the seismic signal, more accurate the time difference between the wave is obtained for each acquisition channel. The results show that this method provides an inhibition to gauss noise and other non coherent noise, and improves the accuracy of calculation effectively.

digital signal processing;wavelet filtering;correlation function;seismic phase extraction;LabVIEW

2015-09-11 責(zé)任編輯:劉憲福

10.13319/j.cnki.sjztddxxbzrb.2016.04.15

河北省國(guó)土資源廳項(xiàng)目(200903)

薛強(qiáng)(1977-)，男，講師，研究方向?yàn)閿?shù)字信號(hào)處理與模式識(shí)別，電能質(zhì)量測(cè)量。E-mail：xueq77@163.com

TH762

A

2095-0373(2016)04-0090-06

薛強(qiáng), 蔡承才, 王碩禾.小波與相關(guān)函數(shù)算法在礦山震源定位中的應(yīng)用[J].石家莊鐵道大學(xué)學(xué)報(bào):自然科學(xué)版,2016,29(4):90-94.