CN103605151A - 基于相位测量的分布式群波浅层微震定位方法 - Google Patents

Abstract

本发明公开了一种基于相位测量的分布式群波浅层微震定位方法，采用分布式震源测试系统，经过布设震动传感器、提取震动信号的相位信息、校正相位信息、节点簇定位及节点群定位五个步骤，最终实现分布式群波浅层微震震源的定位。该方法在采用小区域节点簇布设和大范围节点群任意布设方式，且采用粒子群算法对节点簇内多种优势频率成分对应的相位信息进行自适应定位，采用主成分分析法对定位结果进行数据融合，具有传感器布设难度低、适用范围广、定位精度和可靠性高等优点。

Description

基于相位测量的分布式群波浅层微震定位方法

技术领域

本发明属于爆破震动测试技术、无源定位技术领域，具体涉及一种小区域浅层地下微震震源的定位方法。

背景技术

分布式地下震源定位技术是在监测浅层地下震动的小区域内，利用各节点的全球定位系统（GPS）秒脉冲实现节点之间的同步，并在起震后通过各探测节点拾取微震源产生的震动信号，利用地上控制系统的存储采集模块实现数据的存储。当地下微震过程结束后，地面数据转发系统将各节点的数据发送到中心处理平台，通过微震定位算法，实现震源定位。

与地震、煤矿深层开采、石油勘探等大区域、大当量、大深度震动监测相比较，小区域浅层分布式地下震源定位具有以下特点：（1）分布的范围相对较小，属于小区域监测；（2）地下震动的深度较浅，一般不超过地下100m，属于地下浅层震动；（3）传输介质的密度更加不均匀，往往包括浮土、岩石、沙石、密实土壤等多种类型；（4）对震动位置的定位精度要求更高，往往在1m以内。该技术可以有效地实现地下石油勘探、煤层采空区监测、隧道空洞检测、兵器试验场炮弹炸点定位等，并在工程爆破、国防军事、环境监测、煤矿安全生产等领域有着非常广扩的应用前景。

目前，在上述小区域浅层地下震源定位中，通常借鉴地震震源定位方法，利用拾震器获取震动信息，通过P波初至波的相关识别检测算法（如时间域的长短时平均比法、三分向识别算法、小波变换模极大值识别算法等），获取每个节点P波初至波的到达时间，采用到达时间差法进行震源定位。该方法虽然在地震、石油勘探等大区域深层地下震源定位中较为成熟，但将这种方法应用小区域浅层地下震源定位却存在以下问题：

1、在震动场远场，震动波优势频率稳定、群波速恒定、横纵波叠加现象不明显，因此采用P波初至波到达时间进行到达时间差（Time Difference ofArrival，TDOA）定位时精度较高，相比之下，震动场近场土壤本构特性呈弹塑性，震动波受地面反射、折射影响较大，频散现象严重，即不同频率成分的波形传播速度不同，群波速不是恒定值，因此无法采用现有的TDOA定位算法。

2、震动传感器布设在震源近场时，由于P波、S波、表面波速度不同，造成多种震动波形的混叠现象严重；同时传感器信号线上耦合有爆炸产生的复杂的电磁干扰信号，震动信号信噪比低，造成无法有效地提取精准的初至波到达时间，时差信息误差大。

因此，在小区域浅层地下震源定位中，无法采用现有的地震震源定位方法实现震源定位。

发明内容

本发明的目的是针对以上现有分布式地下震源定位技术所存在的问题，提出一种适应小区域浅层地下微震震源的基于相位测量的分布式群波浅层微震定位方法，该方法具有适用范围广、定位精度高和可靠性高的特点。

为了实现上述本发明的目的，本发明的技术方案是：

基于相位测量的分布式群波浅层微震定位方法，采用分布式震源测试系统，分布式震源测试系统包括控制主站和探测节点群两部分，其中，控制主站包括计算机、授时控制电路及无线传输模块；探测节点群包括至少两个独立的传感器节点簇，传感器节点簇包括五个传感器节点，传感器节点包括震动传感器、信号调理模块、大容量数据采集存储模块、时统模块、探测节点无线传输模块及电源模块，其特征在于：传感器节点簇的五个传感器节点中，四个传感器节点作为测量传感器节点，另一个传感器节点作为参考传感器节点，该方法经过布设震动传感器、提取震动信号的相位信息、校正相位信息、节点簇定位及节点群定位五个步骤，最终实现分布式群波浅层微震震源的定位：

1.布设震动传感器

震动传感器全部埋设于地下，测量传感器节点布设在水平面投影呈任意四边形的四个顶点，参考传感器节点位于四边形内部，在水平投影面上，测量传感器节点距参考传感器节点的最大距离小于震动波最大优势频率对应的波长；震动传感器在空间深度上任意布设，但任意两个传感器节点不能布设在同一深度；

2.提取震动信号的相位信息

对每个震动传感器所获得的微震时域信号，进行小波去噪，对去噪后的信号进行傅立叶变换，得到每个震动传感器所获得微震信号的频域功率谱信息，并计算其对应的优势频带；在优势频带中以频率步进或按照功率谱中所占能量的大小选择优势频率；将所得优势频率进行傅里叶反变换，得到其对应的时域信号；采用互相关算法，求取测量传感器节点与参考传感器节点在优势频率处的相位差；

3.校正相位信息

求出在上述优势频率点处每一个震动传感器的固有相频特性中相位偏移量，并求出各测量传感器节点与所属的节点簇中参考传感器节点的固有相位差；

4.节点簇定位

构建无预测速度条件下的定位方程，利用最小二乘法或chan算法，求得第N组节点簇在优势频率处的迭代初始值；采用粒子群算法求解出第N组节点簇在优势频率处的震源位置；

5.节点群定位

计算所有节点簇的定位结果，利用主成分分析算法进行加权融合，最终获得震源的位置。

本发明与现有技术相比，具有实质性的突出特点和显著效果是：

1、借鉴震动场近场明显的频散现象（即不同的频率成分，其传播速度不同），在小区域节点簇中，通过提取共有优势频率成分，利用相应的相位信息进行定位，解决了震源近场无法采用时差定位的难题。

2、通过提取多个优势频率成分下的相位信息，进行群波时差定位，极大地挖掘了在低信噪比条件下，有效时差信息的数据量，提高了时差信息的可靠度。

3、采用小区域节点簇布设和大范围节点群任意布设方式，使传感器布设难度降低，适用范围更广。

4、采用粒子群算法，通过利用节点簇内多种优势频率成分对应的相位信息进行自适应定位，保证了定位结果的多样性；同时利用主成分分析（PCA）法对定位结果进行数据融合，最大程度提取出了样本中包含的共有的震源特性参数，极大地提高了融合后定位结果的精度和可靠性。

附图说明

图1是本发明分布式震源测试系统结构图；

图2是本发明分布式群波浅层微震定位流程图；

图3是本发明实施例中震动场分布式传感器节点群布设图；

图4是本发明实施例中传感器节点簇布设图；

图5是本发明中震动信号相位提取流程图；

图6是传感器A1经小波处理后三轴合成后的时域图；

图7是传感器B1经小波处理后三轴合成后的时域图；

图8是传感器A1经小波处理后三轴合成后的频域图；

图9是传感器B1经小波处理后三轴合成后的频域图；

图10是传感器A1抽取20Hz后频谱图；

图11是传感器B1抽取20Hz后频谱图；

图12是传感器A1抽取20Hz后时域图；

图13是传感器B1抽取20Hz后时域图；

图14是本发明中相位特征的校正及时差测量流程图；

图15是本发明中节点群波定位流程图。

具体实施方式

以下结合附图详细说明本发明的技术方案。

如图1所示的分布式震源测试系统，包括控制主站和探测节点群两部分，其中，探测节点群包括至少两个独立的传感器节点簇，传感器节点簇包括五个传感器节点，传感器节点包括震动传感器、信号调理模块、大容量数据采集存储模块、时统模块、探测节点无线传输模块及电源模块。

基于相位测量的分布式群波浅层微震定位方法经过布设震动传感器、提取震动信号的相位信息、修正相位信息、节点簇定位及节点群定位五个步骤，最终实现分布式群波浅层微震震源的定位，如图2所示。

1.布设震动传感器

在震源近场，即在装药半径（7-120）倍的范围内布设传感器节点簇，以磁北方向和重力加速度方向组成的笛卡尔坐标系作为布设的统一坐标系，传感器节点簇的五个传感器节点中，四个传感器节点作为测量传感器节点，另一个传感器节点作为参考传感器节点，震动传感器全部埋设于地下。

为了简便，以两个传感器节点簇为例，如图3-4所示。通过2kgTNT在地下起爆产生地下震源，所述的TNT药包位于震动场的坐标原点（0，0，-3）。探测节点群由传感器节点簇1和传感器节点簇2组成，震动传感器A1、B1、C1、D1、E1组成传感器节点簇1，其中震动传感器A1作为参考传感器节点，在坐标系的水平投影面上，四个顶点距参考传感器节点A1的最大距离为3.7m；震动传感器A2、B2、C2、D2、E2组成传感器节点簇2，其中震动传感器A2作为参考传感器节点，在坐标系的投影面上，四个顶点距震动传感器A2的最大距离为3.6m，各震动传感器坐标如表1所示。

表1传感器节点布设位置表单位：m

震动传感器采用中北大学信息探测与处理技术研究所自主研发的地下全向震动传感器，存储设备的采样频率为100KHz，采用GPS秒脉冲保证各传感器节点之间的同步，无线控制触发。

2.提取震动信号的相位信息

分别求取传感器节点簇内四个测量传感器节点与参考传感器节点在优势频率处的相位信息，相位提取流程如图5所示。

①对每个震动传感器所获得的微震时域信号，进行小波去噪，对去噪后的信号进行傅立叶变换，得出每个震动传感器所获得微震信号的频域功率谱信息，即频谱图。在频谱图中计算出对应的优势频带，即每一个频谱图中幅值衰减为最大值的时所对应的频率范围。分别计算传感器节点簇内五个震动传感器的优势频带，并找出五个震动传感器的公共优势频带，设为[f_L,f_H]。

②在所述的优势频带中以一定的频率步进或按照功率谱中所占能量的大小选择优势频率。设f_k为优势频率，优势频率f_k处的单位脉冲为d（w-2pf_k），F_S为系统的采样率。以图4所示的传感器节点簇为例，在优势频率f_k处，对传感器节点簇内的五个震动传感器获取的震动信号进行频谱信息的抽样，抽样后的信号为F_i(w_k)。

$F_i\left(w_k\right)=\underset{-\∞}{\overset{\∞}{\mathrm{\Σ}}}(F_i\left(w\right).\mathrm{\δ}(w-2\mathrm{\π}{f}_k)+F_i\left(w\right).\mathrm{\δ}(w-2\mathrm{\π}(F_s-f_k)))---\left(1\right)$

其中：i=A、B、C、D、E。

③对抽样后的频域信号进行傅里叶反变换，得到其对应的时域信号S_ik(t)，即第i个传感器经第K个优势频率f_k抽取后的时域信号。

S_ik(t)=ifft(F_i(w_k))   （2）

其中，i=A、B、C、D、E。

④设第i个测量传感器节点所属传感器节点簇的参考传感器节点为j，采用互相关算法（式（3）-式（6）），求测量传感器节点B、C、D、E与参考传感器节点A在优势频率f_k处的相位差φ_ij：

$R_{s_i{s}_j}\left(0\right)=\frac{1}{T}{\∫}_0^T{s}_{\mathrm{ik}}\left(t\right)s_{\mathrm{jk}}\left(t\right)\mathrm{dt}---\left(3\right)$

$R_{s_{\mathrm{ik}}}\left(0\right)=\frac{1}{T}{\∫}_0^T{s}_{\mathrm{ik}}^2\left(t\right)\mathrm{dt}---\left(4\right)$

$R_{s_{\mathrm{jk}}}\left(0\right)=\frac{1}{T}{\∫}_0^T{s_{\mathrm{jk}}}^2\left(t\right)\mathrm{dt}---\left(5\right)$

${\mathrm{\φ}}_{\mathrm{ij}}=\arccos \left(\frac{R_{s_{\mathrm{jk}}{s}_{\mathrm{ik}}}\left(0\right)}{\sqrt{R_{s_{\mathrm{ik}}}\left(0\right)R_{s_{\mathrm{jk}}}\left(0\right)}}\right)---\left(6\right)$

其中，s_ik(t)为第i个测量传感器节点经第K个优势频率f_k抽取后的时域信号，s_jk(t)为第j个参考传感器节点经第K个优势频率f_k抽取后的时域信号，T指优势频率对应的周期，

(0)为所述测量传感器节点的时域信号与参考传感器节点的时域信号在时间差为0时的互相关函数值，(0)为所述测量传感器节点的自相关函数最大值，

(0)为所述参考传感器节点的自相关函数最大值，φ_ij为第i个测量传感器节点经频域抽取的时域信号与参考传感器节点j经频域抽取后的时域信号的相位差，i=B、C、D、E，j=A。

以传感器节点簇1中震动传感器A1和震动传感器B1获取的震动信号为例，将各个震动传感器获取的信号经‘dB6’小波分6层后去噪，并将各个传感器的三轴信号进行合成，图6为震动传感器A1经小波处理三轴合成后的时域震动信号，图7为震动传感器B1经小波处理三轴合成后的时域震动信号。

利用FFT变换将各节点经小波处理三轴合成后的信号转换为频域的双边谱图，图8为震动传感器A1经小波处理三轴合成后的局部频谱图，图9为震动传感器B1经小波处理三轴合成后的局部频谱图。结合图8和9可得，两者的共有优势频率在8Hz-130Hz，同时结合震动传感器C1、D1、E1的频谱图，可得传感器节点簇1中各震动传感器的共有优势频带在15Hz-110Hz，传感器节点簇2中各震动传感器的共有优势频带在20Hz-90Hz。

通过频谱分析，在传感器节点簇1中采用式（1）抽取20Hz、40Hz、60Hz和100Hz的频谱成分，如图10-11，采用公式（2）将所述频率下的频谱图转换为相应的时域图，如图12-13。同理抽取传感器节点簇2中30Hz、40Hz、60Hz和80Hz的频域谱信息，计算对应的时域信息。

3.校正相位信息

震动信号相位特征的校正及时差测量步骤如图14所示。

①以图3所示的节点簇为例，根据每一个震动传感器的固有相频特性，求出在上述的优势频率f_k处，每个震动传感器的相位偏移量，并求出各震动传感器与所属的传感器节点簇中参考传感器节点的固有相位差

：

Δφ_ij=φ_i-φ_j   （7）

φ_i为第i个震动传感器在优势频率f_k处的固有相位值,其中i=B、C、D、E，φ_j为所属传感器节点簇中参考传感器节点j在优势频率f_k处的固有相位值，其中j=A。

②修正后的相位为

Δφ_i=|φ_ij+Δφ_ij|   （8）

③第i个震动传感器与所属传感器节点簇中参考传感器在优势频率f_k处的时间差

其中f_k为所述的优势频率，Δφ_i为所述的相位差,其中i=B、C、D、E,j=A。

④利用上述方法求出每个传感器节点簇中四个测量传感器节点与所属传感器节点簇中参考传感器节点在优势频率f_k处的时差信息，以第N组传感器节点簇为例，在优势频率f_k处的时差信息为（τ_n1k,τ_n2k,τ_n3k,τ_n4k）。

在传感器节点簇1中将震动传感器A1设为参考传感器节点，采用式（3）-式（6），计算在20Hz、40Hz、60Hz、100Hz条件下传感器B1、C1、D1、E1与A1的相位差，结合各个传感器的相频特性，采用式（8）修正相位差，同时计算时间差信息，如表2所示。同理，在传感器节点簇2中将震动传感器A2设为参考传感器节点，采用（3）-式（6），计算在30Hz、40Hz、60Hz、80Hz条件下震动传感器B2、C2、D2、E2与A2的相位差，结合各个传感器的相频特性，采用式（8）修正相位差，同时计算时间差信息，如表3所示。

表2节点簇1中传感器B1、C1、D1、E1与A1的相位差和时间差

表3节点簇2中传感器B2、C2、D2、E2与A2的相位差和时间差

4.节点簇定位

基于同源相位信息的群波定位流程如图7所示，首先进行传感器节点簇定位，而后进行节点群定位。

①构建无预测速度条件下的定位方程（式9）：

$\begin{array}{c}\sqrt{{(X_k-x_1)}^2+{(Y_k-y_1)}^2+{(Z_k-z_1)}^2}-\sqrt{{(X_k-x_0)}^2+{(Y_k-y_0)}^2+{(Z_k-z_0)}^2}={\mathrm{v\τ}}_{n1k}\\ \sqrt{{(X_k-x_2)}^2+{(Y_k-y_2)}^2+{(Z_k-z_2)}^2}-\sqrt{{(X_k-x_0)}^2+{(Y_k-y_0)}^2+{(Z_k-z_0)}^2}={\mathrm{v\τ}}_{n2k}\\ \sqrt{{(X_k-x_3)}^2+{(Y_k-y_3)}^2+{(Z_k-z_3)}^2}-\sqrt{{(X_k-x_0)}^2+{(Y_k-y_0)}^2+{(Z_k-z_0)}^2}={\mathrm{v\τ}}_{n3k}\\ \sqrt{{(X_k-x_4)}^2+{(Y_k-y_4)}^2+{(Z_k-z_4)}^2}-\sqrt{{(X_k-x_0)}^2+{(Y_k-y_0)}^2+{(Z_k-z_0)}^2}={\mathrm{v\τ}}_{n4k}\end{array}---\left(9\right)$

式中，（X_k,Y_k,Z_k）为所求震源定位坐标，（x₀,y₀,z₀）为该传感器节点簇重心位置的参考传感器节点坐标，（x₁,y₁,z₁）、（x₂,y₂,z₂）、（x₃,y₃,z₃）、（x₄,y₄,z₄）分别为其余顶点处的4个震动传感器的坐标。

②采用粒子群算法求解出第N组传感器节点簇在优势频率f_k处的震源位置P_Nk=(X'_Nk,Y'_Nk,Z'_Nk)。

③重复步骤①-②，计算所有节点簇的定位结果P_ik=(X_ik,Y_ik,Z_ik)，其中下标i代表传感器节点簇序号，下标k代表对应传感器节点簇的第k个优势频率，将所有定位结果组成定位结果样本P_i=(X_i,Y_i,Z_i)，1≤i≤M,其中M为定位结果总数。

采用粒子群算法针对无预测速度条件下的定位模型利用表2、表3的数据进行定位，粒子群参数：Genration=20；算法重复次数max_iterm=300；迭代次数sizepop=40；种群数目scale1=5；X、Y搜索范围scale2=[-50]；Z搜索范围w_max=0.9、w_min=0.4；惯性权重。

得出传感器节点簇1、2的定位结果见表4。

表4各节点簇的定位结果

优化节点群样本

按照表4，构建节点群定位结果样本P=[P₁P₂P₃P₄P₅P₆P₇P₈]

第1组定位结果P₁=[-0.1247,0.1358,-3.0175]

第2组定位结果P₂=[-0.0202,0.0182,-2.9128]

第3组定位结果P₃=[0.0624,-0.0474,-3.8341]

第4组定位结果P₄=[-0.1621,-0.2902,-3.0891]

第5组定位结果P₅=[0.0086,-0.4245,-2.9200]

第6组定位结果P₆=[0.2271,0.6207,-2.8631]

第7组定位结果P₇=[-0.5576,-2.2282,-2.7198]

第8组定位结果P₈=[-0.0191,0.1183,-2.7303]

样本的期望E(P)=[-0.0732,-0.2621,-3.0108]

样本里的元素距期望的标准差L=[0.4013,0.3017,0.8615,0.1217,0.2032,0.9441,2.045,0.4757]

经过分析P₃,P₆,P₇误差较大，将其剔除，进行下一步数据融合。

5.节点群定位

①计算定位结果样本P_i的期值，同时计算上述定位结果样本P_i距期望的距离L_i，1≤i≤M。按照偏离期望的远近程度剔除样本中的粗大误差值，将其余值进行下一步的数据融合。

②采用PCA算法进行加权融合，最终计算出震源的定位结果P。

$P=\underset{i=1}{\overset{M-H}{\mathrm{\Σ}}}{\mathrm{\λ}}_i*P_i---\left(10\right)$

式中，λ_i为对应定位结果的权值，P_i为第i组定位结果。

通过对样本P构建协方差矩阵，计算得到最大特征值对的特征向量为[-0.4725,-0.4552,-0.4477,-0.4239,-0.4345]，通过归一化得到对应样本中每个元素的系数为[0.2114,0.2037,0.2003,0.1899,0.1944]。通过式（10）得到：

P=0.2144*P₁+0.2037*P₂+0.2003*P₃+0.1899*P₄+0.1944*P₅=[-0.0650,-0.0833,2.9361]

最后，本发明计算得到的震源结果为X=-0.0650m,Y=-0.0833m,Z=-2.9361m。实际的震源位置为X=0m,Y=0m,Z=-3m，均方根误差为12.3cm，满足近场的定位要求。

Claims (1)

1.基于相位测量的分布式群波浅层微震定位方法，采用分布式震源测试系统，分布式震源测试系统包括控制主站和探测节点群两部分，其中，控制主站包括计算机、授时控制电路及无线传输模块；探测节点群包括至少两个独立的传感器节点簇，传感器节点簇包括五个传感器节点，传感器节点包括震动传感器、信号调理模块、大容量数据采集存储模块、时统模块、探测节点无线传输模块及电源模块，其特征在于：传感器节点簇的五个传感器节点中，四个传感器节点作为测量传感器节点，另一个传感器节点作为参考传感器节点，该方法经过布设震动传感器、提取震动信号的相位信息、校正相位信息、节点簇定位及节点群定位五个步骤，最终实现分布式群波浅层微震震源的定位：

1）布设震动传感器

震动传感器全部埋设于地下，测量传感器节点布设在水平面投影呈任意四边形的四个顶点，参考传感器节点位于四边形内部，在水平投影面上，测量传感器节点距参考传感器节点的最大距离小于震动波最大优势频率对应的波长；震动传感器在空间深度上任意布设，但任意两个传感器节点不能布设在同一深度；

2）提取震动信号的相位信息

对每个震动传感器所获得的微震时域信号，进行小波去噪，对去噪后的信号进行傅立叶变换，得到每个震动传感器所获得微震信号的频域功率谱信息，并计算其对应的优势频带；在优势频带中以频率步进或按照功率谱中所占能量的大小选择优势频率；将所得优势频率进行傅里叶反变换，得到其对应的时域信号；采用互相关算法，求取测量传感器节点与参考传感器节点在优势频率处的相位差；

3）校正相位信息

求出在上述优势频率点处每一个震动传感器的固有相频特性中相位偏移量，并求出各测量传感器节点与所属的节点簇中参考传感器节点的固有相位差；

4）节点簇定位

构建无预测速度条件下的定位方程，利用最小二乘法或chan算法，求得第N组节点簇在优势频率处的迭代初始值；采用粒子群算法求解出第N组节点簇在优势频率处的震源位置；

5）节点群定位

计算所有节点簇的定位结果，利用主成分分析算法进行加权融合，最终获得震源的位置。

Images