煤礦井下微震震源高精度定位研究

程久龍， 宋廣東,,3*， 劉統(tǒng)玉， 胡賓鑫， 王紀(jì)強(qiáng)， 王金玉

1 中國(guó)礦業(yè)大學(xué)(北京)煤炭資源與安全開采國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 100083 2 山東省科學(xué)院激光研究所， 濟(jì)南 250014 3 山東微感光電子有限公司， 濟(jì)南 250103

?

煤礦井下微震震源高精度定位研究

程久龍1， 宋廣東1,2,3*， 劉統(tǒng)玉2， 胡賓鑫2， 王紀(jì)強(qiáng)2， 王金玉2

1 中國(guó)礦業(yè)大學(xué)(北京)煤炭資源與安全開采國(guó)家重點(diǎn)實(shí)驗(yàn)室， 北京 100083 2 山東省科學(xué)院激光研究所， 濟(jì)南 250014 3 山東微感光電子有限公司， 濟(jì)南 250103

煤礦井下微震震源準(zhǔn)確定位，對(duì)于動(dòng)力災(zāi)害監(jiān)測(cè)預(yù)警具有重要意義．由于微震震源需要通過井下傳感器接收信息反演確定，傳感器的安裝位置限制于煤礦井下巷道周圍，傳感器沿巷道近平面的不合理布置將大大降低震源定位精度．針對(duì)由傳感器信息反演震源位置引起的病態(tài)問題，本文提出了基于微震監(jiān)測(cè)測(cè)點(diǎn)優(yōu)化布置的震源高精度定位算法．首先通過計(jì)算系數(shù)矩陣條件數(shù)，判定病態(tài)問題；然后利用中心化法和行平衡法聯(lián)合進(jìn)行病態(tài)矩陣預(yù)處理．對(duì)預(yù)處理后的矩陣A、b利用L曲線法計(jì)算正則參數(shù)，結(jié)合Tikhonov正則化算法計(jì)算得到震源坐標(biāo)正則解．研究結(jié)果表明，中心化法有效降低了矩陣數(shù)量級(jí)，行平衡預(yù)處理降低了病態(tài)條件數(shù)，預(yù)處理后Tikhonov正則解的震源坐標(biāo)誤差最小可以達(dá)到3.09 m，與預(yù)處理前的高斯消去解相比誤差大大降低．通過上述優(yōu)化處理，實(shí)現(xiàn)了井下受限空間微震監(jiān)測(cè)震源高精度定位．

震源定位反演； 測(cè)點(diǎn)布置優(yōu)化； 微震； Tikhonov正則化； 高精度定位

1 引言

準(zhǔn)確定位震源位置是利用微震監(jiān)測(cè)系統(tǒng)對(duì)煤礦井下動(dòng)力災(zāi)害進(jìn)行監(jiān)測(cè)預(yù)警的關(guān)鍵問題之一．通過在煤礦井下布置傳感器，接收地下巖層破裂產(chǎn)生的震動(dòng)波信號(hào)，結(jié)合傳感器坐標(biāo)、巖層波速等數(shù)據(jù)可實(shí)現(xiàn)震源位置反演．而傳感器的布置受煤礦井下狹小巷道空間的制約，受此影響震源反演定位的誤差將會(huì)大大增加．

針對(duì)煤礦井下工作面的微震監(jiān)測(cè)，傳感器只能安裝在順槽周圍(葉根喜等，2008)．為使傳感器避免布置在同一水平面內(nèi)，需要鉆深孔安裝，增加測(cè)點(diǎn)間垂向差值．但是長(zhǎng)鉆孔卡鉆、塌孔等問題影響傳感器的順利安裝．井下施工條件惡劣，深孔安裝意味著施工難度的加大和成本費(fèi)用的提高，同時(shí)鉆孔深度是有限的．在地表進(jìn)行礦井采空區(qū)安全監(jiān)測(cè)和相鄰礦井防盜采微震監(jiān)測(cè)也面臨同樣問題．因而受限空間微震測(cè)點(diǎn)優(yōu)化布置和已有測(cè)點(diǎn)監(jiān)測(cè)數(shù)據(jù)非適定條件即病態(tài)條件下的震源高精度定位問題需要進(jìn)行深入研究．

目前震源定位方法有十多種，空間定位精度在5～50 m之間(林峰等，2010；董隴軍等，2011)，其中線性定位法是其中的一類重要分支．震源坐標(biāo)、傳感器坐標(biāo)與初始到時(shí)可以構(gòu)成二次非線性方程組，把非線性方程組通過降維轉(zhuǎn)化為線性方程組(李寅祺和祝永新，2009)，可以用高斯消去法(李慶揚(yáng)等，2008)求得精確解．但當(dāng)線性方程組系數(shù)矩陣條件數(shù)很大時(shí)，用高斯消去法求方程解誤差將會(huì)變大．當(dāng)系數(shù)矩陣A的元素?cái)?shù)值量級(jí)差別很大時(shí)，可采用行均衡或列均衡的方法，這時(shí)矩陣A的條件數(shù)可得到改善(劉建國(guó)，1999)．Geiger法、聯(lián)合定位法、相對(duì)定位法和雙重殘差法等線性定位方法都沒有對(duì)因測(cè)點(diǎn)布置引起的病態(tài)問題進(jìn)行研究(田玥和陳曉非，2002；金明培等，2014)，且當(dāng)只有有限測(cè)點(diǎn)監(jiān)測(cè)到震動(dòng)信號(hào)時(shí)，震源定位誤差將大大增加．

微震震源定位是根據(jù)監(jiān)測(cè)到時(shí)、波速等信息計(jì)算震源位置的反演問題，即數(shù)學(xué)中的反問題．Tikhonov 和 Arsenin (1977)和Phillips于20世紀(jì)60年代分別獨(dú)立提出了求數(shù)學(xué)反問題穩(wěn)定近似解的正則化方法．正則化方法包括截?cái)嗥娈愔捣纸?、改進(jìn)的截?cái)嗥娈愔捣纸狻V義截?cái)嗥娈愔捣纸?、衰減奇異值分解、最小熵原理、最大熵正則化、二次約束最小二乘法和截?cái)嗤耆钚《?Hansen，2007)．目前應(yīng)用較廣的正則化方法是Tikhonov法．正則化參數(shù)是正則化理論中的重要參數(shù)，正則化參數(shù)的選擇是正則化法的基礎(chǔ)(王彥飛，2007)．正則化參數(shù)λ控制賦予解范數(shù)和殘差范數(shù)最小化的權(quán)重，平衡不穩(wěn)定性和光滑性．正則化參數(shù)的選取分為兩類：(1)誤差估計(jì)，即偏差原則(馬成業(yè)等，2010)，這種簡(jiǎn)單方法的主要問題是當(dāng)選擇的參數(shù)過大時(shí)，得到的是具有巨大正則化誤差的過正則化解；(2)從右端項(xiàng)提取必要的信息，包括L曲線準(zhǔn)則，廣義交叉驗(yàn)證(Generalized Cross-Validation，GCV)和近似最優(yōu)準(zhǔn)則．GCV有時(shí)會(huì)把相關(guān)噪聲誤認(rèn)為是有效信號(hào)；而L曲線可以極好的分辨相關(guān)干擾(Hansen and O′Leary,1993；Hansen，2001)．逄煥東等(2004)針對(duì)微地震線性方程定位過程中的病態(tài)問題，采用Tikhonov正則化法進(jìn)行迭代計(jì)算，發(fā)現(xiàn)Tikhonov正則化法比最小二乘法和正規(guī)解法收斂速度更快．曾小牛等(2015)利用λ和截止波數(shù)ωc的關(guān)系，得到了改進(jìn)的正則低通濾波函數(shù)，成功應(yīng)用于位場(chǎng)各階垂向?qū)?shù)的換算．

最方便的圖形化分析離散非適定問題的工具是L曲線，包括了全部的有效正則參數(shù)．L曲線是使解范數(shù)和殘差范數(shù)同時(shí)最小化的折中方法，顯示了這些量對(duì)正則化參數(shù)的依賴程度(Hansen，1997)．Miller(1970)、Lawson and Hanson(1974)最早利用L曲線的方法處理非適定最小二乘問題．L曲線在聯(lián)系Tikhonov正則化和病態(tài)問題時(shí)起到了核心作用．只有有限測(cè)點(diǎn)監(jiān)測(cè)到震動(dòng)信號(hào)且測(cè)點(diǎn)布置近乎在同一平面時(shí)，由測(cè)點(diǎn)坐標(biāo)、到時(shí)等信息構(gòu)造的震源求解矩陣很可能是病態(tài)的，基于L曲線準(zhǔn)則的Tikhonov正則化法可提高此病態(tài)矩陣求解精度．

本文針對(duì)煤礦井下巷道受限空間微震監(jiān)測(cè)中存在的傳感器沿巷道、近乎平面布置問題，分析影響震源定位精度的因素，在對(duì)傳感器布置進(jìn)行優(yōu)化的基礎(chǔ)上，基于現(xiàn)場(chǎng)測(cè)試數(shù)據(jù)特點(diǎn)，對(duì)病態(tài)問題提出了先進(jìn)行預(yù)處理再求正則解的方法．提出了對(duì)傳感器安裝位置原始坐標(biāo)進(jìn)行中心化的算法，以降低系數(shù)矩陣A、b之間的數(shù)值量級(jí)差距；構(gòu)建了行平衡預(yù)處理的算子矩陣，以降低系數(shù)矩陣A的條件數(shù)；對(duì)經(jīng)上述中心化和行平衡預(yù)處理后的數(shù)據(jù)進(jìn)行奇異值分解(Singular Value Decomposition，SVD)，利用L曲線法計(jì)算正則條件數(shù)，再計(jì)算Tikhonov正則解，并驗(yàn)證了正則解的適用性．

2 理論基礎(chǔ)

設(shè)未知震源大地坐標(biāo)為(x，y，z)，第i個(gè)傳感器坐標(biāo)為xi，yi，zi．如圖1所示，震源與第i個(gè)傳感器Si之間的距離為

=Ki-2xix-2yiy-2ziz+x2+y2+z2,

(1)

圖1 震源與傳感器位置關(guān)系Fig.1 The position relation between seismic source and sensors

通過把二次非線性方程組降維變?yōu)橐淮尉€性方程組，可以求得震源坐標(biāo)的精確解．設(shè)ri,1為ri與r1差值，r1是距離震源最近傳感器S1與震源之間的距離(李寅祺和祝永新，2009)，則：

(2)

其中，xi,1=xi-x1，yi,1=yi-y1，zi,1=zi-z1，r1=(t1-t0)v，ri,1=(ti-t0)v-(t1-t0)v=(ti-t1)v，如圖1所示，t1是地震波到達(dá)S1的時(shí)間，ti是地震波到達(dá)Si的時(shí)間，t0為微震發(fā)生的時(shí)間是未知量，r1是t0的函數(shù)，v為波速．寫成矩陣的形式為

(3)

為論述方便上式可以簡(jiǎn)寫為

AX=b X∈Rn×n,b∈Rn,

(4)

對(duì)于線性方程組，設(shè)A為非奇異系數(shù)矩陣，則A的條件數(shù)為(李慶揚(yáng)，2008)：

cond(A)v=‖A-1‖v‖A‖v(v=1,2或∞).

(5)

當(dāng)A的條件數(shù)遠(yuǎn)大于1，即cond(A)?1時(shí)，則線性方程組是病態(tài)的，即A是病態(tài)矩陣．A的條件數(shù)越大，線性方程組越病態(tài)，用常規(guī)方法求得的解偏差越大(李慶揚(yáng)等，2008；肖庭延等，2003)．通常使用的條件數(shù)為cond(A)∞(劉建國(guó)，1999)，因而本文在計(jì)算中采用式cond(A)∞計(jì)算線性方程組系數(shù)矩陣條件數(shù)．

3 病態(tài)問題處理方法

測(cè)點(diǎn)布置是系數(shù)矩陣A病態(tài)與否的關(guān)鍵因素，為實(shí)現(xiàn)對(duì)目標(biāo)區(qū)域巖體活動(dòng)高精度監(jiān)測(cè)，首先對(duì)測(cè)點(diǎn)布置進(jìn)行優(yōu)化；對(duì)優(yōu)化后依然病態(tài)的問題，采用預(yù)處理加正則化處理的方法進(jìn)行求解．

3.1 線性無(wú)關(guān)優(yōu)化

為實(shí)現(xiàn)震源三維空間定位，測(cè)點(diǎn)布置時(shí)，四個(gè)以上的測(cè)點(diǎn)不要布置在同一水平面上，因測(cè)點(diǎn)間z坐標(biāo)相同，系數(shù)矩陣A的行列式detA=0，A為奇異矩陣，故A-1不存在，將導(dǎo)致方程組有無(wú)窮解或無(wú)解，即要保證detA>0．

系數(shù)矩陣A的行向量或列向量的近似相關(guān)性，即向量間夾角非常小時(shí)，是矩陣病態(tài)的直接原因．因而測(cè)點(diǎn)在空間三個(gè)方向都要不等間距設(shè)置，盡量無(wú)序排列，以便于構(gòu)建非奇異系數(shù)矩陣A，最終使基于直接法解線性方程組成為可能．即要滿足：ai≠q aj和ci≠q cj，其中ai、aj是行向量，ci，cj代表列向量，q是任意實(shí)數(shù)．系數(shù)矩陣A兩列或兩行成比例，趨于完全重合時(shí)，條件數(shù)趨于正無(wú)窮大．

為監(jiān)測(cè)足夠?qū)拸V的區(qū)域，在測(cè)點(diǎn)平面坐標(biāo)(xi，yi)基本確定的條件下，增大測(cè)點(diǎn)之間的安裝深度差，可以有效降低矩陣A的條件數(shù)．

3.2 初始坐標(biāo)中心化

如果常數(shù)項(xiàng)b的微小變化，引起方程組AX=b解的巨大變化，此方程組即為病態(tài)方程組．設(shè)b有擾動(dòng)δb，相應(yīng)解X的擾動(dòng)記為δx(劉建國(guó)，1999)，則：

A(X+δx)=b+δb,

(6)

δx=A-1δb.

(7)

對(duì)某一確定的微震事件，ri,1，xi,1，yi,1，zi,1是定值，大小決定于測(cè)點(diǎn)間的相對(duì)位置關(guān)系，不受測(cè)點(diǎn)原始坐標(biāo)取值大小的影響，但式(4)右端項(xiàng)可表示為

b= (x1xi,1+y1yi,1+z1zi,1)

(8)

從上式可以看出，x1，y1，z1的取值對(duì)b影響巨大．若初始坐標(biāo)采用8位整數(shù)大地坐標(biāo)表示x1、y1，將造成矩陣A和b之間數(shù)量級(jí)的巨大差距．因而需要把b轉(zhuǎn)化成與A一樣不受初始坐標(biāo)影響的常數(shù)項(xiàng)，保證在正則化求解中減小解的誤差．

對(duì)于給定的測(cè)點(diǎn)初始坐標(biāo)(xi，yi，zi)，本文提出了如下的中心化處理方法：

(9)

(10)

3.3 行平衡預(yù)處理

受制于監(jiān)測(cè)范圍寬廣和節(jié)省成本簡(jiǎn)單快速安裝的原因，xi,1或yi,1與zi,1數(shù)量級(jí)通常相差巨大．采用如下方法對(duì)式(4)進(jìn)行平衡化預(yù)處理，以降低系數(shù)矩陣的條件數(shù)．選擇非奇異矩陣P∈Rn×n，使式(4)化為等價(jià)方程組：

PAX=P b,

(11)

(12)

(13)

式(12)，是劉建國(guó)(1999)提出的，記為mmax．

本文提出了下面的基于2范數(shù)的對(duì)角矩陣構(gòu)造法，記為m2，具體如下：

mi=‖ai‖2,

(14)

ai為矩陣A的行向量，即A=[a1,a2,…,an]T，且存在‖PA‖2=1．在后面的數(shù)據(jù)處理部分，將比較兩種方法的實(shí)際應(yīng)用效果，選擇最優(yōu)的構(gòu)造方法．

3.4 正則化反演

式(4)中A是歐氏空間中的有界線性算子，是緊算子，當(dāng)detA≠0時(shí)方程組存在唯一解，條件數(shù)較大時(shí)，問題是不適定的，即病態(tài)的．求解病態(tài)問題，需要在b變化不大的情況下，尋求穩(wěn)定的近似解．正則化解xλ的定義(Tikhonov and Arsenin, 1977)為

(15)

λ為正則化參數(shù)，控制相對(duì)邊界約束和殘差范數(shù)最小化的權(quán)重，L為邊界限制矩陣，x*為估計(jì)初值．

Tikhonov 是直接化法，正則化模型即求解下面的最小二乘問題(ELDEN，1977)：

(16)

系數(shù)矩陣A的奇異值分解為

(17)

ui是U中的行向量，vi是V中的行向量，i是行標(biāo)．

標(biāo)準(zhǔn)形式的Tikhonov正則化模型，L=In，x*=0，正則解的計(jì)算方法(Hansen，2001)為

(18)

(19)

3.5 利用L曲線準(zhǔn)則計(jì)算λ

‖xλ‖2和‖Axλ-b‖2在判定λ時(shí)起到了重要的作用．過正則化作用于解時(shí)，將導(dǎo)致殘差‖Axλ-b‖2太大；然而當(dāng)欠正則化時(shí)，殘差‖Axλ-b‖2較小，但是正則解將受數(shù)據(jù)誤差的控制，導(dǎo)致‖xλ‖2將會(huì)很大．當(dāng)兩個(gè)量都需要控制時(shí)L曲線事實(shí)上是一條折中的曲線．

λ取為一組常量，代入式(18)和(19)，求得xλ和b-A xλ，并以‖Axλ-b‖2為橫坐標(biāo)，以‖xλ‖2為縱坐標(biāo)，得到L曲線．如圖2所示．

圖2 L曲線Fig.2 The L-curve

L曲線準(zhǔn)則，即在L拐角處最大曲率對(duì)應(yīng)的正則化參數(shù)即為最優(yōu)正則化參數(shù)．從圖中無(wú)法直接得到精確的最優(yōu)正則化參數(shù)，圖2中最優(yōu)正則化參數(shù)為0.039833，通過如下方法求解其精確數(shù)值解．

計(jì)算L曲線的曲率k(Hansen and O′Leary，2001)：

(20)

(21)

以k為縱坐標(biāo)，λ為橫坐標(biāo)做圖，曲率最大值對(duì)應(yīng)的λ即為最優(yōu)正則參數(shù)．

3.6 病態(tài)問題處理步驟

測(cè)點(diǎn)布局階段進(jìn)行線性無(wú)關(guān)優(yōu)化，從源頭上避免病態(tài)問題．隨著工作面推進(jìn)，可用測(cè)點(diǎn)減少，剩余測(cè)點(diǎn)或優(yōu)化后測(cè)點(diǎn)構(gòu)成的系數(shù)矩陣不可避免的出現(xiàn)病態(tài)問題時(shí)，需要進(jìn)行預(yù)處理和正則化求解．具體數(shù)據(jù)處理步驟如下：

(1) 選定初始到時(shí)最小的測(cè)點(diǎn)為基準(zhǔn)點(diǎn)，另外選定4個(gè)測(cè)點(diǎn)，構(gòu)建式(4)；

(2) 計(jì)算系數(shù)矩陣A，并求cond(A)∞，若cond(A)∞<10，則用高斯消去法直接求解，并結(jié)束計(jì)算；假如cond(A)∞≥10，則執(zhí)行步驟(3)；

(3) 進(jìn)行中心化預(yù)處理和行平衡；

(4) 進(jìn)行正則化處理，得到正則近似解．

4 現(xiàn)場(chǎng)試驗(yàn)

4.1 測(cè)點(diǎn)布置

某礦屬特厚煤層礦井，采礦方法為綜采放頂法，為適時(shí)監(jiān)測(cè)巖體破裂，預(yù)測(cè)礦井動(dòng)力災(zāi)害，在正在回采工作面安裝了一套微震監(jiān)測(cè)系統(tǒng)，共11個(gè)測(cè)點(diǎn)布置于軌道順槽頂?shù)装鍍?nèi)，測(cè)點(diǎn)布置如圖3所示，黑色填充圓點(diǎn)表示安裝在頂板的傳感器，正方形代表安裝在底板的傳感器，標(biāo)定炮孔在底板中．各測(cè)點(diǎn)坐標(biāo)及P波初始到時(shí)數(shù)據(jù)如表1所示．

表1 測(cè)點(diǎn)坐標(biāo)及到時(shí)數(shù)據(jù)

表1中S1、S2…S11為各測(cè)點(diǎn)編號(hào)，xi、yi、zi為各測(cè)點(diǎn)三維坐標(biāo)，ti為P波初至?xí)r刻．標(biāo)定炮的空間坐標(biāo)為(3899.89，5097.48，1009.30)，S1為距離標(biāo)定炮最近的基點(diǎn)，任取兩個(gè)傳感器Si和Sj，二者到標(biāo)定炮的距離分別為di和dj，震動(dòng)波到時(shí)分別為ti和tj，則基于上述兩個(gè)傳感器的波速計(jì)算公式為

(22)

分別取4組底板測(cè)點(diǎn)組合和4組頂板測(cè)點(diǎn)組合，選取間隔最大的測(cè)點(diǎn)、跨越軌順的測(cè)點(diǎn)、間隔最近的測(cè)點(diǎn)等分別計(jì)算對(duì)應(yīng)的波速vij，再對(duì)全部vij算術(shù)平均，得到波速為4.35 m·ms-1．具體如表2所示．

圖3 測(cè)點(diǎn)布置圖(a) 平面圖(x-y坐標(biāo))； (b) 剖面圖(y-z坐標(biāo)).Fig.3 Monitoring points arrangement graph(a) Planar graph (x-y coordinates); (b) Profile graph (y-z coordinates).

表2 波速計(jì)算

4.2 嚴(yán)重病態(tài)問題求解

通過設(shè)計(jì)傾角、深度變化的鉆孔實(shí)現(xiàn)測(cè)點(diǎn)三維方向差異性優(yōu)化布置．但仍然受限于巷道狹小空間，底板測(cè)點(diǎn)z坐標(biāo)最大間距僅為1.64 m，而頂板測(cè)點(diǎn)z坐標(biāo)最大間距為10.56 m．但y方向測(cè)點(diǎn)間距最大值為174 m．對(duì)因空間布置引起的病態(tài)問題，需要進(jìn)行分析處理．由于系統(tǒng)布置時(shí)未作線性無(wú)關(guān)優(yōu)化，因而S1、S3、S7、S10基本在同一水平面內(nèi)．在上述底板測(cè)點(diǎn)參與的震源定位反演中，下面四組數(shù)據(jù)出現(xiàn)了嚴(yán)重病態(tài)問題．各組合包含的測(cè)點(diǎn)編號(hào)為：組合一(S1、S3、S4、S7、S11)，組合二(S1、S3、S9、S10、S11)，組合三(S1、S3、S6、S7、S10)，組合四(S1、S4、S7、S9、S11)，其中S1為距離標(biāo)定炮最近的基點(diǎn)．具體方程解如表3所示．各組合分別采用m2、mmax構(gòu)造的對(duì)角矩陣進(jìn)行行平衡，計(jì)算的條件數(shù)如圖4所示．分別采用高斯消去算法(高斯消去)、正則化算法(正則化)、中心化處理后再求正則化解(中心化正則化)和

中心化加行平衡預(yù)處理再求正則化解(預(yù)處理正則化)得到的四組解的誤差如圖5所示．

從圖4可以看出，各組合的三類條件數(shù)基本是遞減的，行平衡處理前條件數(shù)最大，采用m2構(gòu)造的對(duì)角矩陣計(jì)算的條件數(shù)降低最明顯．從圖5可以看出，四種算法得到的解誤差也呈遞減變化規(guī)律，高斯消去解誤差最大，預(yù)處理正則化法震源定位反演效果最好．

表3 四組病態(tài)問題數(shù)據(jù)優(yōu)化處理

圖4 各組合條件數(shù)Fig.4 Condition number of four groups

圖5 不同算法解誤差Fig.5 The error of different algorithms

不做預(yù)處理，條件數(shù)較大，四組震源組合中，最小的條件數(shù)為6062.69，最大值為14834.2都屬于病態(tài)方程組．直接用高斯消去法進(jìn)行震源定位處理，得到的震源坐標(biāo)誤差最小為2411.34 m，最大為10558.34 m，結(jié)果都是沒有意義的．

不做預(yù)處理，直接用L曲線計(jì)算正則化參數(shù)、用Tikhonov正則化算法計(jì)算得到正則解，相比高斯消去法，對(duì)應(yīng)組合的解都得到了改善，誤差最小為2370.15 m，最大為2514.29 m，仍然與實(shí)際結(jié)果相差巨大，沒有應(yīng)用價(jià)值．

通過分析初始系數(shù)矩陣A和b的特點(diǎn)發(fā)現(xiàn)，A中元素絕對(duì)值最大量級(jí)為102而b中元素絕對(duì)值最大量級(jí)為106，初步分析是數(shù)量級(jí)差別大造成的解的巨大誤差．對(duì)四組數(shù)據(jù)進(jìn)行中心化預(yù)處理，處理后矩陣A沒有變化，因而條件數(shù)沒有降低，右端項(xiàng)b量級(jí)明顯降低，其元素絕對(duì)值最大量級(jí)為104，大多數(shù)為103．對(duì)中心化預(yù)處理后的數(shù)據(jù)，若用高斯消去法求解，定位誤差將得不到改善．利用Tikhonov正則化解法，大大提高了定位精度，組合一和組合四定位誤差縮小到10 m以內(nèi)，最小定位誤差僅為7.28 m，平均定位誤差為14.4 m，定位結(jié)果已經(jīng)達(dá)到可以實(shí)際應(yīng)用的水平．

對(duì)上述四組數(shù)據(jù)進(jìn)行中心化預(yù)處理后繼續(xù)進(jìn)行行平衡預(yù)處理，cond(A)∞得到明顯改善，最大值由14834.2降為10818．如圖4所示，組合一、二、四采用m2構(gòu)造的對(duì)角矩陣計(jì)算的條件數(shù)都比采用mmax構(gòu)造的對(duì)角矩陣計(jì)算的條件數(shù)小，只有組合三的條件數(shù)二者基本接近，因而本文提出的以2范數(shù)構(gòu)造對(duì)角矩陣進(jìn)行行平衡預(yù)處理效果更好．但是方程組仍然是病態(tài)的，需要進(jìn)行正則化處理．行平衡預(yù)處理后，組合一和組合二的Tikhonov正則解改善明顯，組合一誤差最小僅為3.09 m，組合三誤差減小了0.1 m，只有組合四的誤差增大了1.06 m．相比中心化正則化法定位精度得到了進(jìn)一步提升，平均定位誤差減小為11.6 m．測(cè)點(diǎn)布置、信號(hào)品質(zhì)和平均速度等因素，都會(huì)對(duì)解誤差產(chǎn)生影響．

5 結(jié)論

(1) 采用線性方法進(jìn)行震源定位反演，需要對(duì)測(cè)點(diǎn)布置做線性無(wú)關(guān)優(yōu)化，4個(gè)以上的測(cè)點(diǎn)不要布置在同一水平面上，測(cè)點(diǎn)布置三維方向非等間距，做到無(wú)序排列，增大測(cè)點(diǎn)之間的安裝深度差．

(2) 測(cè)點(diǎn)初始坐標(biāo)取值是影響震源定位結(jié)果誤差水平的最主要因素．基于微震數(shù)據(jù)的特點(diǎn)提出了中心化預(yù)處理法，在不影響系數(shù)矩陣A的情況下，極大的降低了常數(shù)項(xiàng)b的數(shù)量級(jí)，定位誤差大大降低．中心化前沒有任何意義的定位結(jié)果，中心化后得到的定位結(jié)果達(dá)到了可以實(shí)際應(yīng)用的水平．

(3) 提出了基于系數(shù)矩陣A行向量2范數(shù)構(gòu)建對(duì)角矩陣P對(duì)方程組進(jìn)行行平衡預(yù)處理，有效降低了系數(shù)矩陣A的條件數(shù)．但是條件數(shù)的改善是有限度的，對(duì)于達(dá)不到優(yōu)化預(yù)期的嚴(yán)重病態(tài)方程組，必須要進(jìn)行正則化處理．

(4) 若直接利用微震監(jiān)測(cè)數(shù)據(jù)求正則化解，則定位誤差依然很大．本文提出的對(duì)微震監(jiān)測(cè)數(shù)據(jù)進(jìn)行預(yù)處理后求正則化解的方法，可極大的提高解的精度．最小定位誤差可以達(dá)到3.09 m，四組數(shù)據(jù)中有三組的震源定位反演誤差小于10 m，平均定位誤差為11.6 m，達(dá)到了實(shí)際工程應(yīng)用要求．

總之，測(cè)點(diǎn)優(yōu)化布置、高質(zhì)量采樣信號(hào)、速度模型優(yōu)化等因素是獲得高精度震源位置的基礎(chǔ)．本文從常規(guī)方法認(rèn)為無(wú)效的無(wú)法處理的監(jiān)測(cè)數(shù)據(jù)中提取出了有用的信息，在有限測(cè)點(diǎn)監(jiān)測(cè)到信號(hào)的情況下，也能得到震源精確位置．本文為微震震源定位病態(tài)問題提供了一種全新的反演算法，促進(jìn)了微震監(jiān)測(cè)技術(shù)的發(fā)展．

DongL J, Li X B, Tang L Z, et al. 2011. Mathematical functions and parameters for microseismic source location without pre-measuring speed.ChineseJournalofRockMechanicsandEngineering(in Chinese), 30(10): 2057-2067.

Eldén L. 1977. Algorithms for the regularization of ill-conditioned least squares problems.BITNumericalMathematics, 17(2): 134-145.

Hansen P C, O'Leary D P. 1993. The use of the L-curve in the regularization of discrete ill-posed problems.SIAMJ.Sci.Comput., 14(6): 1487-1503.

Hansen P C. 1997. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. Philadelphia: SIAM.

Hansen P C. 2001. The L-Curve and its Use in the Numerical Treatment of Inverse Problems. Southampton: WIT Press.

Hansen P C. 2007. Regularization Tools Version 4.0 for Matlab 7.3.NumericalAlgorithms, 46(2): 189-194.

Jin M P, Wang R J, Dai S G. 2014. Improvement and case test of linear method for single event seismic location.ActaSeismologicaSinica(in Chinese), 36(5): 757-759. Lawson C L, Hanson R J. 1974. Solving Least Squares Problems. Englewood Cliffs: SIAM.

Li Q Y, Wang N C, Yi D Y. 2008. Numerical Analysis. 5th ed. (in Chinese). Beijing: Tsinghua University Press.

Li Y Q, Zhu Y X. 2009. Capsule in vivo environment's TDoA wireless location algorithm and simulation.MicrocomputerInformation(in Chinese), 25(22): 134-135.

Lin F, Li S L, Xue Y L, et al. 2010. Microseismic sources location methods based on different initial values.ChineseJournalofRockMechanicsandEngineering(in Chinese), 29(5): 996-1002. Liu J G. 1999. Matrix Condition Number and Ill-Posed Improvement [Master thesis] (in Chinese). Chongqing: Chongqing University. Ma C Y, Yang S L, Li S P. 2010. A regularization method of ill-conditioned linear equations.JournalofGansuScience(in Chinese), 22(4): 33-35.

Miller K. 1970. Least squares methods for ill-posed problems with a prescribed bound.SIAMJ.Math.Anal., 1(1): 52-74. Pang H D, Jiang F X, Zhang X M. 2004. Study on nonhomogeneous material′s AE by image processing method.RockandSoilMechanics(in Chinese), 25(Z1): 60-62. Tian Y, Chen X F. 2002. Review of seismic location study.ProgressinGeophysics(in Chinese), 17(1): 147-155. Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems. Washington: Winston. Wang Y F. 2007. Computational Methods for Inverse Problems and Their Applications (in Chinese). Beijing: Higher Education Press.

Xiao T Y, Yu S G, Wang Y F. 2003. Numerical Solution of Inverse Problem (in Chinese). Beijing: Science Press.

Ye G X, Jiang F X, Yang S H. 2008. Possibility of automatically picking first arrival of microseismic wave by energy eigenvalue method.ChineseJ.Geophys. (in Chinese), 51(5): 1574-1581. Zeng X N, Li X H, Jia W M, et al. 2015. A new regularization method for calculating the vertical derivatives of the potential

field.ChineseJ.Geophys. (in Chinese), 58(4): 1400-1410,doi: 10.6038/cjg20150426.

附中文參考文獻(xiàn)

董隴軍, 李夕兵, 唐禮忠等. 2011. 無(wú)需預(yù)先測(cè)速的微震震源定位的數(shù)學(xué)形式及震源參數(shù)確定. 巖石力學(xué)與工程學(xué)報(bào), 30(10): 2057-2067.

金明培, 汪榮江, 戴仕貴. 2014. 線性單事件地震定位方法的改進(jìn)及應(yīng)用嘗試. 地震學(xué)報(bào), 36(5): 757-759.

李慶揚(yáng), 王能超, 易大義. 2008. 數(shù)值分析(第5版). 北京: 清華大學(xué)出版社.

李寅祺, 祝永新. 2009. 體內(nèi)膠囊TDoA定位算法研究與仿真. 微計(jì)算機(jī)信息, 25(22): 134-135.

林峰, 李庶林, 薛云亮等. 2010. 基于不同初值的微震源定位方法. 巖石力學(xué)與工程學(xué)報(bào), 29(5): 996-1002.

劉建國(guó). 1999. 矩陣條件數(shù)及病態(tài)改善[碩士論文]. 重慶: 重慶大學(xué). 馬成業(yè), 楊勝良, 黎鎖平. 2010. 求解病態(tài)線性方程組的一個(gè)正則化方法. 甘肅科學(xué)學(xué)報(bào), 22(4): 33-35.

逄煥東, 姜福興, 張興民. 2004. 微地震的線性方程定位求解及其病態(tài)處理. 巖土力學(xué), 25(增刊): 60-62.

田玥, 陳曉非. 2002. 地震定位研究綜述. 地球物理學(xué)進(jìn)展, 17(1): 147-155.

王彥飛. 2007. 反演問題的計(jì)算方法及其應(yīng)用. 北京: 高等教育出版社.

肖庭延, 于慎根, 王彥飛. 2003. 反問題的數(shù)值解法. 北京: 科學(xué)出版社.

葉根喜, 姜福興, 楊淑華. 2008. 時(shí)窗能量特征法拾取微地震波初始到時(shí)的可行性研究. 地球物理學(xué)報(bào), 51(5): 1574-1581.

曾小牛, 李夕海, 賈維敏等. 2015. 位場(chǎng)各階垂向?qū)?shù)換算的新正則化方法. 地球物理學(xué)報(bào), 58(4): 1400-1410, doi: 10.6038/cjg20150426.

(本文編輯 汪海英)

High precision location of micro-seismic sources in underground coal mine

CHENG Jiu-Long1, SONG Guang-Dong1,2,3*, LIU Tong-Yu2, HU Bin-Xin2, WANG Ji-Qiang2, WANG Jin-Yu2

1StateKeyLaboratoryofCoalResourcesandSafeMining,ChinaUniversityofMining&Technology,Beijing100083,China2InstituteofLaser,ShandongAcademyofSciences,Ji′nan250014,China3ShandongMicro-sensorPhotonicsLimited，Ji′nan250103，China

Accurate location of micro-seismic sources in underground coal mine has great significance for monitoring and early warning of dynamic disaster. Micro-seismic source information is generally extracted through inversion of the data acquired by the underground sensors. The installation of the sensors is limited around underground roadway. The seismic location precision will be greatly decreased in the case of unreasonable arrangement of sensors along the roadway and approximately in a plane. As for the ill-posed problems caused by the inversion for source location with sensor-received information, high precision location algorithm of micro-seismic source based on monitoring points arrangement optimization is presented in this paper. Firstly, the ill-posed problems are judged by calculating the coefficient matrix condition numbers. The ill-posed matrix is then pretreated by centralization and row balance jointly. The regularization parameters of pretreated matrix A and b are calculated by the L-curve method. The seismic source coordinates regularization solution is obtained using Tikhonov regularization algorithm. The results show that the matrix magnitude is decreased effectively with centralization method, ill-posed condition numbers are reduced with row balance pretreatment. The minimum seismic source coordinate error of Tikhonov regularization solution after pretreatment is 3.09 m, and is greatly decreased compared to the error by Gaussian elimination solution before pretreatment. High precision seismic source location is realized in limited space of underground roadway by means of the above optimization processing.

Seismic source location inversion； Monitoring points arrangement optimization； Micro-seismic； Tikhonov regularization； High precision location

10.6038/cjg20161214.

國(guó)家重點(diǎn)研究發(fā)展計(jì)劃項(xiàng)目(2016YFC0801405)，國(guó)家自然科學(xué)基金項(xiàng)目(51574250)，山東省科技發(fā)展計(jì)劃項(xiàng)目(2014GSF120017)資助．作者簡(jiǎn)介 程久龍，男，1965年生，教授，博士生導(dǎo)師，主要從事應(yīng)用地球物理科研與教學(xué)工作.E-mail:jlcheng@126.com

10.6038/cjg20161214

P631

2016-08-04，2016-09-30收修定稿

程久龍， 宋廣東， 劉統(tǒng)玉等. 2016. 煤礦井下微震震源高精度定位研究. 地球物理學(xué)報(bào)，59(12):4513-4520,

Cheng J L, Song G D, Liu T Y, et al. 2016. High precision location of micro-seismic sources in underground coal mine.ChineseJ.Geophys. (in Chinese),59(12):4513-4520,doi:10.6038/cjg20161214.

*通訊作者 宋廣東，男，1981年生，中國(guó)礦業(yè)大學(xué)(北京)地球探測(cè)與信息技術(shù)專業(yè)在讀博士，研究方向?yàn)榈V山動(dòng)力災(zāi)害監(jiān)測(cè)與防治研究.

E-mail:guangdongsgd@163.com