CN104820240A - 一种基于WSNs的地震层析成像方法 - Google Patents

Abstract

本发明公开了一种基于WSNs的地震层析成像方法，适用于地质勘探、地震监测及诊断领域。包括以下具体过程：用传感器捕获地震信号，通过工作于传感器节点上软件进行P波检测、P波识别和P波初至时间拾取，传感器节点把拾取后的P波初至时间转发给基站，所有捕获到P波的传感器节点联合在一起进行地震震源定位计算，基站根据地下地质体超立方体仿真辅助模型、P波穿越地下地质体的数学模型对其轨迹进行跟踪，通过求解上述数学模型获得等价描述地下地质体之物性参数的解，用最后的解来重构图像，从而实现层析成像。本发明能够满足实时成像、现场成像的应用需求，并能显示地下地质体在地震时的动态变化情况。

Description

一种基于WSNs的地震层析成像方法

技术领域

本发明属于地质勘探、地震监测与诊断领域，特别涉及一种基于WSNs结合地球物理专业知识领域的地震层析成像方法。

背景技术

根据地球物理学知识可知，地震发生时会产生很多种类型的波，其中P波会首先到达地球表面，因此如何快速地识别P波并捕获之，然后拾取其初至时间，根据P波的识别和初至时间来进行地震的检测和预警，在地球物理学界已经形成共识，并有比较长的研究历史。

传统的方法主要是通过地震仪来捕获地震信号，并对捕获的地震数据进行事后处理，存在以下几点不足：①满足不了实时处理、现场处理的要求；②不同的地震仪之间独立工作、相互之间协同性比较差，很难保证数据之间的关联性；③实现不了对地下地质体的实时层析成像，更实现不了实时地显示成像的结果，甚至地下地质体动态变化情况；④传统方法是采用数据统一存储、统一处理的模式，必将遇到大数据处理的瓶颈问题。

随着无线传感器网络(Wireless Sensor Networks，简称WSNs)技术的日趋成熟，大范围、大面积监控地质对象的需求完全可以实现；另外，由于传感器节点计算能力的显著提高，使基于无线传感器网络平台的地震层析成像技术的实现完全成为可能。

发明内容

鉴于现有地震监测、地震层析成像技术的不足，本发明提供一种效率更高、实现成本更低，能够满足实时、现场处理要求的基于WSNs的地震层析成像方法。

本发明的技术方案如下：

一种基于WSNs的地震层析成像方法，包括如下步骤：

①用传感器捕获地震信号，通过工作于传感器节点上的软件算法首先对地震信号进行滤波处理，然后进行P波检测、P波识别和P波初至时间拾取，传感器节点把拾取后的P波初至时间转发给基站；

②所有捕获到P波的传感器节点联合在一起进行地震震源定位计算；

③基站根据地下地质体超立方体仿真辅助模型、地震射线穿越地下地质体的数学模型对地震射线的传播轨迹进行跟踪；

④通过求解上述数学模型获得能够等价描述地下地质体之物性参数的解，用最终的解实现图像重构及层析成像。

在步骤①中，所述对地震信号进行滤波处理、P波检测、P波识别的具体方式为：首先对地震信号进行卡尔曼滤波处理，衰减环境噪声，提高地震信号的信噪比；在此P波检测、P波识别的具体方式提供三种方法：方法一，把地震信号看成沿时间轴展开的时间序列，用时间序列的异常检测来发现和识别P波；方法二，根据地震信号包括频谱、相移、信号短时窗能量、信号长时窗能量、P波发生的地理位置在内的属性特征，通过对这些属性进行训练，用人工智能、模式识别的方法进行P波的识别；方法三，对地震信号用低通滤波器进行滤波，根据信号处理结果进行P波识别。

在步骤①中，所述P波初至时间拾取的具体方式为：在P波检测、P波识别的同时，如果是用时间序列异常检测的方法进行P波识别，则检测到异常的时刻就为P波的初至时间；如果用人工智能、模式识别的方法进行P波识别，则识别P波的时刻视为P波初至时间；如果用信号处理的方法进行P波识别，假设捕获到的地震信号服从泊松分布，对其进行最大似然估计计算，根据地震波到达时信号发生变化的时刻确定P波初至时间。

在步骤②中，所述地震震源定位计算的具体方式为：计算所有捕获到P波的传感器节点所监测的圆形覆盖区域的重叠部分之几何中心，即为地震震源的位置。

在步骤③中，所述地下地质体超立方体仿真辅助模型为：包含多个小立方体，每一个独立的小立方体代表一个地质体，通过多个小立方体的组合仿真任意形状和大小的地下地质体的组成和结构；布局在地表的多个传感器节点组成WSNs，通过WSNs定位地震震源并估算其深度，并在模型中标记地震震源；用地震震源到传感器节点之间的连线表示地震射线的传播轨迹。

在步骤③中，所述地震射线穿越地下地质体的数学模型为：As＝t，其中，A是地震射线从上一个地质体穿越到下一个地质体时所经过的距离，s是地震射线在穿越该地质体时的速度的倒数，t是穿越该地质体所用的时间。

在步骤④中，所述图像重构及层析成像的具体方式为：求解每一个地震震源向外发散出的多根地震射线传播轨迹的数学模型，得到每一根地震射线所穿越每一个地下地质体时的“s”，在此“s”是一个与地质体密度有密切关系的参数，根据每一层的“s”来重构图像，达到层析成像的效果。

本发明的有益技术效果是：

与传统的通过事后处理地震数据的层析成像方法相比较，本发明基于WSNs的地震层析成像方法，在传感器节点之间的协同性、同一时间所采集的数据对同一事件的描述的可靠性、计算平台的计算能力等方面都有显著地提高，从而满足实时成像、现场成像的应用需求，并能显示地下地质体在地震时的动态变化情况。

本发明的优点将在下面具体实施方式部分的描述中给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

图1是用无线传感器网络(WSNs)进行地震层析成像的原理示意图。

图2是传感器节点的工作原理图。

图3是地下地质体超立方体仿真辅助模型。

图4是工作于基站的层析成像软件的计算过程。

具体实施方式

下面结合附图对本发明的具体实施方式做进一步说明。

本发明是一种利用WSNs进行地震层析成像的方法，其实现原理如图1所示。主要步骤包括：①通过传感器捕获地震信号，通过工作于传感器节点上的软件首先对地震信号进行滤波处理，然后进行P波检测、P波识别和P波初至时间拾取，传感器节点把拾取后的P波初至时间转发给基站；②所有捕获到P波的传感器节点联合在一起进行地震震源定位计算；③基站根据地下地质体超立方体仿真辅助模型、地震射线穿越地下地质体的数学模型对地震射线的传播轨迹进行跟踪；④通过求解上述数学模型获得能够等价描述地下地质体之物性参数的解，用最终的解实现图像重构及层析成像。

本发明的实现需要依赖一个完整的平台，平台从整体上而言，由感知层、无线传感器网络和工作于基站的层析成像软件三层体系结构组成，即平台包括硬件和软件两大部分。

布局在地表的传感器组成了平台的感知层，负责采集地震信号；无线传感器网络负责数据的处理和转发；基站负责存储各种模型，并由工作于基站的层析成像软件进行地震震源定位、地震射线的传播轨迹跟踪、图像重构相关的计算。

用“Beaglebone Black Board(BBB)”与“MSP430F6779”相结合的方式设计传感器节点，如图2所示。MSP430F6779的接口电路部分与GPS定位、Xbee通信、SD Micro存储卡以及Geophone传感器等连接，中间件层软件工作于Beaglebone Black Board(BBB)，主要完成的工作包括：P波识别、P波检测、P波初至时间拾取和参与地震震源定位计算，并把拾取后的P波初至时间转发给基站。

一、P波识别、P波检测和P波初至时间拾取。

传感器捕获到的地震信号，根据不同的传感器工作原理，可以是电压或磁场能量或加速度大小等不同的物性参数，但不论是哪一种性质的信号，首先用卡尔漫滤波法滤除信号中的环境噪声等干扰数据。另一方面，由地球物理学的常识可知，正常情况下所采集的地质环境背景数据是平坦的波形，一旦有地震发生，波形的振幅会由缓慢到剧烈的增加、然后又由剧烈的增加变化成缓慢最后恢复至平坦。在波形变化之初，即其振幅由平坦开始有增加阶段说明有P波到达。从数据处理的角度看，地震信号是一个沿时间轴展开的时间序列，其中P波的检测、识别和初至时间拾取可以看成是一个数据挖掘的过程，可通过时间序列的异常检测方法来检测和识别P波，检测到异常的时间点可认为是P波的初至时间；地震信号具有各种属性特征，如频谱、相移、信号短时窗能量、信号长时窗能量、P波发生的地理位置，通过对这些属性进行训练，用人工智能、模式识别等前沿理论知识研究对P波的智能识别方法，用P波被成功识别的时刻作为P波的初至时间；从信号处理的角度，可用低通滤波方法来检测和识别P波，检测到P波的那个时间点可认为是P波的初至时间；还可以通过在假设其服从一定的分布模型的前提下，对其进行最大似然估计计算，然后统计波形的变化，用最开始变化的时间点来表示P波的初至时间。

传感器节点把P波初至时间和当前传感器节点所在的地理位置(经度值、纬度值)封装成Xbee数据帧，并把该数据帧通过无线传感器网络转发给基站。

二、震源定位。

所有捕获到P波的传感器节点进行地震震源位置计算。每个传感器节点所监测的范围，即其覆盖区域是一个圆形面积。通过计算这些圆形覆盖区域的公共部分，即重叠部分，这些重叠部分的几何中心就是地震震源的位置。

三、地震射线轨迹跟踪。

图3示出了地下地质体超立方体仿真辅助模型，以8个子立方体为例，其中的每个子立方体还可以继续细分为8个更小的子立方体，如些类推、以至无穷。每一个独立的小立方体代表一个地质体，所以：①该超立方体模型可以包含无穷多个小立方体；②该辅助模型完全可以仿真任意形状和大小的地下地质体的组成和结构。图中“星号”代表地震震源；“三角形”代表布局在地表的传感器节点；立方体顶部的数字代表布局在地表的传感器节点的编号；从地震震源到传感器节点之间的连线表示地震射线传播轨迹。从图中可以看出，每条地震射线都穿越了多个不同的地质体后方能到达地表。布局在地表的多个传感器节点组成WSNs。

参照图3来具体说明本发明的实现方法。

假设超立方体仿真模型平均划分为C₁至C₈部分，其中C₂，C₄，C₆，C₈的顶部代表地表，每个小立方体上布局了4个传感器节点见图3。地震震源位于C₅，从左往右由地震震源向地球表面发射了四条地震射线，记为a₁，a₂，a₃，a₄，并且假设射线a₁穿越地质体(更细粒度的超立方体)C₅，C₁，C₂后被传感器节点1所捕获；射线a₂穿越C₅，C₁，C₃，C₄后被传感器节点6所捕获；射线a₃穿越C₅，C₆后被传感器节点12所捕获；射线a₄穿越C₅，C₇，C₈后被传感器节点16所捕获。并且假设{s₁，s₂，s₃，s₄，s₅，s₆，s₇，s₈}是地震射线分别穿越{C₁，C₂，C₃，C₄，C₅，C₆，C₇，C₈}的速度之倒数。根据地球物理学领域的地震射线“走时残差”理论和地震射线传播的速度模型，可得如下方程。表述如下：

$\left[\begin{array}{cccccccc}{a}_{\mathrm{1,1}}& a_{\mathrm{1,2}}& 0& 0& a_{\mathrm{1,5}}& 0& 0& 0\\ a_{\mathrm{2,1}}& 0& a_{\mathrm{2,3}}& a_{\mathrm{2,4}}& a_{\mathrm{2,5}}& 0& 0& 0\\ 0& 0& 0& 0& a_{\mathrm{3,5}}& a_{\mathrm{3,6}}& 0& 0\\ 0& 0& 0& 0& a_{\mathrm{4,5}}& 0& a_{\mathrm{4,7}}& a_{\mathrm{4,8}}\end{array}\right]$

$\begin{array}{c}\·{[s_1,s_2,s_3,s_4,s_5,s_6,s_7,s_8]}^T\\ =[t_1,t_2,t_3,t_4]\end{array}---\left(1\right)$

其中矩阵元素如“a_1，1”表示地震射线a₁在地质体C₁中所经过的距离，“a_4，7”表示地震射线a₄在地质体C₇中所经过的距离。t₁，t₂，t₃，t₄分别代表地震射线a₁，a₂，a₃，a₄在地质体中的传播时间。这个示例说清楚了以下几点：①方法不仅可以用于对地下地质体进行层析，还可以用于勘探地下地质体的物质密度，即可用于地质勘探领域；②虽然对于布局在地球表面的传感器而言，地震射线在地下地质体中传播的实际距离和所经过地质体的物质密度都是未知的，但是当只能在地表获得地震射线到达的时间的情况下，完全可以通过反演的方式逆向推算另外两个地下地质体的物性参数“a”和“s”；③对于真实环境而言，超立方体模型的规模相当大，即速度模型中的矩阵规模也相应得增加；④在本发明中，通过传感器网络捕获得到P波的初至时间，在此，初至时间并不是地震射线真正的传播时间。显然，在实际应用需求过程中，地震射线的轨迹跟踪要依赖于复杂的数学方法来进行求解。

对式(1)进行一般化处理后如式(2)所示，式(2)就是跟踪地震射线传播轨迹的速度模型

As＝t  (2)

其中A,s,t的含义同式(1)。

为了使方法更体现它的一般性和可用性，对上述方法进行进一步扩展。从而给出求解类似问题的通用数学方法。

震源所产生的地震射线穿越不同的地质体，由于不同地质体的物质“密度”不一定相同，所以穿越不同地质体时的时间也不相同，而每一条这样的地震射线从震源到被布局在地面的传感器节点所捕获，会经历一定的时间，这个时间表示成t_i＝[t_i1，t_i2，…，t_ij]^T，其中t_ij表示第“j”个地震震源往第“i”个传感器节点所传播地震射线所经历的时间。显然，这个时间从宏观上看是地面观测到的时间(即P波初至时间)与地震发生时的精确时间之差。问题的关键是在实际应用中很难精确地测定地震在地下地质体中精确的发生时间。

为了解决上述问题，首先假设两个时间：y^*表示观测到的地震射线传播时间，t⁰表示所预测的地震射线传播时间。则：表示传感器节点i所捕获到的第j个地震源的地震射线所经历的传播时间，其中表示第j个地震源的地震射线在往第i个传感器节点传播时所经历第m个地下地质体时的时间，A_i[j，m]表示地震射线在第m个地质体中所经历的距离,s^*[m]表示第m个地质体影响地震射线传播的物性参数，在速度模型中可认为是传播速度；另外，表示所预测的时间，显然在此满足根据速度模型式(2)，地震射线传播时间满足方程(3)和(4)。

A_is^*-A_is⁰＝A_is  其中A_is＝t_i  (3)

上式等价于

As＝t  (4)

其中A_i[j,m]∈A_i，A_i∈R^i×M， $A={[A_1^T,A_2^T,\·\·\·,A_i^T,\·\·\·,A_N^T]}^T,t={[t_1^T,t_2^T,\·\·\·,t_i^T,\·\·\·,t_N^T]}^T.$

式(4)是一个线性方程。这样问题变成如何利用诸如最小方差法求解上述线性方程。求解的过程简要描述成式(5)所示。

s＝arg  minimize||t-As||²  (5)

方程(4)、方程(5)的求解过程中另一个关键问题就是的确定，通常用P波的初至时间表示，需要估算。

以下给出一个有关的估算方法。结合图3，从上往下对超立方体沿纵向分割，即图3中的超立方体可分成四列，每一列包括上下两部分，即如图3所示，每一列包括两个更小的立方体。当然，当该立方体模型的粒度进一步变小时，列数和每一列包括的小立方体个数就会相应地增多。总之，通过这样分割后，假设每一列就是一个部分，记为“E_p”，那么地震射线a_l穿越E_p的估计时间用数学上的Kaczmarz方法计算如下，如式(6)所示。

$t_{l,p}=T_l\·\frac{{\stackrel{\→}{a}}_{l,p}\·{\hat{s}}_p\left(0\right)}{a_l,s^0}-{\stackrel{\→}{a}}_{l,p}\·{\hat{s}}_p\left(0\right)---\left(6\right)$

其中，表示地震射线a_l在E_p中穿越的部分，表示E_p中的s,T_l可以用传感器网络中相关节点捕获到的P波初至时间代替。

四、层析成像。

参照图3，当超立方体模型由无数个更细粒度的小立方体组成时。这样的超立方体沿“X-Y”平面或“Y-Z”平面或“X-Z”平面的任意一个横截面就对应地下地质体在某一经度或纬度上的某个层面。至此，层析成像问题就转换成了通过数学优化的方法求解方程(5)中的“s”，用“s”来实现图像重构。重构的图像就是层析成像的结果。

上述在无线传感器网络的基站进行层析成像计算的流程如图4所示。

以上所述的仅是本发明的优选实施方式，本发明不限于以上实施例。可以理解，本领域技术人员在不脱离本发明的基本构思的前提下直接导出或联想到的其他改进和变化，均应认为包含在本发明的保护范围之内。

Claims (7)

1.一种基于WSNs的地震层析成像方法，其特征在于，包括如下步骤：

①用传感器捕获地震信号，通过工作于传感器节点上的软件对所采集的地震信号进行去除噪声处理，然后进行P波检测、P波识别和P波初至时间拾取，传感器节点把拾取后的P波初至时间转发给基站；

②所有捕获到P波的传感器节点联合在一起进行地震震源定位计算；

③基站根据地下地质体超立方体仿真辅助模型、地震射线穿越地下地质体的数学模型对地震射线的传播轨迹进行跟踪；

④通过求解上述数学模型获得能够等价描述地下地质体之物性参数的解，用最终的解实现图像重构及层析成像。

2.根据权利要求1所述的基于WSNs的地震层析成像方法，其特征在于，在步骤①中，对所采集的地震信号进行去除噪声处理，然后进行P波检测、P波识别处理的方法为：首先对地震信号进行卡尔曼滤波处理，衰减环境噪声，提高地震信号的信噪比；P波检测、P波识别的具体方式包括：把地震信号看成沿时间轴展开的时间序列，用时间序列的异常检测来发现和识别P波；或者，根据地震信号包括频谱、相移、信号短时窗能量、信号长时窗能量、P波发生的地理位置在内的属性特征，通过对这些属性进行训练，用人工智能、模式识别的方法进行P波的识别；或者，对地震信号用低通滤波器进行滤波，根据信号处理结果进行P波识别。

3.根据权利要求1所述的基于WSNs的地震层析成像方法，其特征在于，在步骤①中，所述P波初至时间拾取的具体方式为：在P波检测、P波识别的同时，用时间序列异常检测的方法进行P波识别，则检测到异常的时刻就为P波的初至时间；或者，用人工智能、模式识别的方法进行P波识别，则P波被识别的时刻可视为P波初至时间；或者，用信号处理的方法进行P波识别，假设捕获到的地震信号服从泊松分布，对其进行最大似然估计计算，根据地震波到达时信号发生变化的时刻确定P波初至时间。

4.根据权利要求1所述的基于WSNs的地震层析成像方法，其特征在于，在步骤②中，所述地震震源定位计算的具体方式为：计算所有捕获到P波的传感器节点所监测的圆形覆盖区域的重叠部分之几何中心，即为地震震源的位置。

5.根据权利要求1所述的基于WSNs的地震层析成像方法，其特征在于，在步骤③中，所述地下地质体超立方体仿真辅助模型为：包含多个小立方体，每一个独立的小立方体代表一个地质体，通过多个小立方体的组合仿真任意形状和大小的地下地质体的组成和结构；布局在地表的多个传感器节点组成WSNs，通过WSNs定位地震震源并估算其深度，并在模型中标记地震震源；用地震震源到传感器节点之间的连线表示地震射线的传播轨迹。

6.根据权利要求1所述的基于WSNs的地震层析成像方法，其特征在于，在步骤③中，所述地震射线穿越地下地质体的数学模型为：As＝t，其中，A是地震射线从上一个地质体穿越到下一个地质体时所经过的距离，s是地震射线在穿越该地质体时的速度的倒数，t是穿越该地质体所用的时间。

7.根据权利要求1所述的基于WSNs的地震层析成像方法，其特征在于，在步骤④中，所述图像重构及层析成像的具体方式为：求解每一个地震震源向外发散出的多根地震射线传播轨迹的数学模型，得到每一根地震射线所穿越每一个地下地质体时的“s”，在此“s”是一个与地质体密度有密切关系的参数，根据每一层的“s”来重构图像，达到层析成像的效果。