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

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

?

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

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

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

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

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

1 引言

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

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

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

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

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

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

2 理論基礎(chǔ)

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

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

(1)

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

通過把二次非線性方程組降維變?yōu)橐淮尉€性方程組，可以求得震源坐標的精確解．設(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是地震波到達S1的時間，ti是地震波到達Si的時間，t0為微震發(fā)生的時間是未知量，r1是t0的函數(shù)，v為波速．寫成矩陣的形式為

(3)

為論述方便上式可以簡寫為

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

(4)

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

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

(5)

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

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

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

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

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

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

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

3.2 初始坐標中心化

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

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

(6)

δx=A-1δb.

(7)

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

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

(8)

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

對于給定的測點初始坐標(xi，yi，zi)，本文提出了如下的中心化處理方法：

(9)

(10)

3.3 行平衡預處理

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

PAX=P b,

(11)

(12)

(13)

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

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

mi=‖ai‖2,

(14)

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

3.4 正則化反演

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

(15)

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

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

(16)

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

(17)

ui是U中的行向量，vi是V中的行向量，i是行標．

標準形式的Tikhonov正則化模型，L=In，x*=0，正則解的計算方法(Hansen，2001)為

(18)

(19)

3.5 利用L曲線準則計算λ

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

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

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

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

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

(20)

(21)

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

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

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

(1) 選定初始到時最小的測點為基準點，另外選定4個測點，構(gòu)建式(4)；

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

(3) 進行中心化預處理和行平衡；

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

4 現(xiàn)場試驗

4.1 測點布置

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

表1 測點坐標及到時數(shù)據(jù)

表1中S1、S2…S11為各測點編號，xi、yi、zi為各測點三維坐標，ti為P波初至時刻．標定炮的空間坐標為(3899.89，5097.48，1009.30)，S1為距離標定炮最近的基點，任取兩個傳感器Si和Sj，二者到標定炮的距離分別為di和dj，震動波到時分別為ti和tj，則基于上述兩個傳感器的波速計算公式為

(22)

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

圖3 測點布置圖(a) 平面圖(x-y坐標)； (b) 剖面圖(y-z坐標).Fig.3 Monitoring points arrangement graph(a) Planar graph (x-y coordinates); (b) Profile graph (y-z coordinates).

表2 波速計算

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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

附中文參考文獻

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

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

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

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

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

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

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

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

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

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

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

曾小牛, 李夕海, 賈維敏等. 2015. 位場各階垂向?qū)?shù)換算的新正則化方法. 地球物理學報, 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.

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

10.6038/cjg20161214

P631

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

程久龍， 宋廣東， 劉統(tǒng)玉等. 2016. 煤礦井下微震震源高精度定位研究. 地球物理學報，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年生，中國礦業(yè)大學(北京)地球探測與信息技術(shù)專業(yè)在讀博士，研究方向為礦山動力災(zāi)害監(jiān)測與防治研究.

E-mail:guangdongsgd@163.com