0 引言

在水力压裂微震监测方法中, 微地震“事件”反演定位方法的科学性、算法的稳定性、计算精度和效率等, 始终是微震监测方法研究及其工程技术研究所需面对的核心技术问题.

目前流行的单事件线性定位算法大都源于Geiger(1912) 提出的经典定位方法, 如HYPO系列定位程序.在此基础上又发展了多种定位方法, 如多事件定位方法(震源与台站校正的联合反演, 震源与速度结构的联合反演)(Spence, 1980)、相对定位法(主事件定位法), 空间域定位方法、非线性定位方法(牛顿法、震源搜索算法)、双差定位算法(Waldhauser and Ellsworth, 2000)等等.从数学原理上讲, 地震定位的实质是求目标函数的极小值.上述大多数方法在求解过程中都要计算偏导数、逆矩阵、解方程组等, 若台阵布局不合理、观测资料不理想时, 会出现求解失败、迭代不收敛等问题.

震源扫描算法作为“穷举法”的一种全局最优算法, 由于其适应性较强, 计算无需求偏导及逆矩阵、不用解方程等优势, 仍被广泛用于地震定位中.

传统的震源扫描算法是在整个时空内安照一定的尺度将解空间网格化, 在反演过程中计算所有网格点处的微地震记录在时间域的能量“亮度”, 将“亮度”最大的点作为定位解(Kao and Shan, 2004).震源扫描算法无需拾取微震事件, 属于直接波形反演, 从计算方法来看, 属于约束条件下的全局最优化问题.震源扫描算法的反演精度和可靠性严重地依赖于划分网格的尺度大小, 尺度越小, 精度越高, 反之亦然.但是该算法由于要遍历解空间的所有网格点, 所以计算量随着搜索范围的增大及网格的细化, 计算量越来越大, 特别是对于多参数优化问题, 计算量随着网格尺寸的减小而呈几何级数量增加.因此, 传统震源扫描算法不适合大批量事件的多参数反演.Jones等人(1993) 首次提出了DIRECT算法, 该算法用来寻找已知解空间内的多变量函数的全局最优解, 既不需要明确的目标函数表达式也不需要计算导数, 而是依据每次迭代中采样点的函数值和超矩形体的细分情况决定下一步搜索的区域, 非常适合求解具有确定变量空间的“黑盒”最优问题.DIRECT算法在求解问题时, 无需具体划分网格, 因此解的精度不依赖于网格尺寸.搜索过程中, 该算法在真实解附近会加大搜索密度, 而在远离真实解附近会减小搜索密度.因此, 搜索次数会远少于传统网格搜索算法.

本文基于传统网格搜索算法的思想, 在求解震源最优解时采用DIRECT搜索算法, 提出了基于DIRECT算法的快速网格搜索震源扫描定位方法.通过理论的微震地面监测和井中监测模型测试对比, 在相同条件下, 基于DIRECT算法的震源扫描定位方法比传统震源扫描算法定位精度高、搜索次数少、收敛速度快, 是一种高效、准确的定位算法, 能更好地适用于微震监测大数据实时处理技术需求.

1 震源扫描定位算法

震源扫描算法(Source Scanning Algorithm, SSA)是Kao和Shan(2004) 提出的.提出该方法的理论基础是：微地震有效事件的信号在传播特点上或形式上类似于常规地震勘探中的绕射波, 借鉴偏移的思想由扫描叠加偏移发展而来.震源扫描算法是一种可以对震源分布进行成像的方法, 这种方法不需要预先知道断层产状和规模, 它可以方便地利用相对振幅、到时等波形数据检测给定的时间和位置上是否有微地震震源.通过系统地扫描一个试探性的位置范围和起震时刻, 可以发现整个震源序列分布而不需要拾取到时, 又无需计算理论地震图.

震源扫描算法属于直接网格搜索定位方法, 虽然该方法工作量大, 但它对于求解具有相对少量参数反演问题, 是一种行之有效方法, 可以准确地收敛到最优解, 且在均匀介质或层状介质情况下, 搜索速度也可以接受.而当介质模型比较复杂时, 单纯使用直接网格搜索法, 则需要耗费巨大的工作量和计算时间, 使得此法不能满足实际使用的要求.下面简单介绍震源扫描算法的基本原理.

震源扫描算法就是在整个可能的时间空间内搜寻可能的震源位置的方法.基本原理可用图 1描述.

假设某一微震事件由N个(图 1中N=3, A、B、C点)台站记录到, 首先对每个台站接收到的地震记录u_n进行归一化处理, 然后按式(1) 计算某个空间点(η)在某个时刻(τ)的“亮度”函数为

(1)

式中, N是接收台站个数, u_n是第n个台站接收微震信号的归一化记录, τ是微震事件的发震时刻, t_ητ为从点η处到台站n的走时.如果所有的最大振幅都是由点η和时间τ产生, 则b_η(G, S)=1.同样, b_η(G, S)=0.1意味着此处只有10%的可能性为震源.通常情况下, b_η(G, S)=1的情况很少见, 可以认为整个时空中“亮度”函数取最大时的空间及时间参数, 分别是其震源位置和事件发生的时刻.

在实测数据反演过程中, 由于速度误差等原因, 会造成计算走时与理论走时有一定的偏差, 一些学者对“亮度”函数进行了改进, 如：李文军和陈棋福(2006) 在“亮度”函数中加入光滑函数; Kao和Shan(2007) 在“亮度”函数中加入时窗长度因子等.

利用网格搜索算法求解上述问题时, 需要首先给定解的存在空间, 然后按照一定的尺度将解空间进行网格划分, 最后通过遍历所有网格点来寻找整个空间中的全局最优解, 即“亮度”值最大对应的解.由传统网格搜索法原理可知, 其反演精度和可靠性严重依赖于划分网格时所用尺度的大小; 尺度越小, 反演精度越高; 尺度越大, 反演精度越低.但是尺度越小, 网格数越大, 计算量随之巨增.特别是对于多参数优化问题, 计算量会随着尺度的减小而呈几何级数增加.微地震方法中待优化的参数包括速度(V)、震源空间位置(X、Y、Z)以及震源激发时刻(t₀)等; 求解时, 每个参数都需给出解空间, 并按一定尺度进行网格划分, 因此计算量非常大.虽然一些学者提出了一些减少计算量的策略, 但总体来说, 随着网格尺度的减小, 计算量依然很大.

上述震源扫描法是基于能量扫描的, 在实际应用中也可以用拾取走时与计算走时差作为目标函数, 利用震源扫描算法来求解目标函数的最优解.公式为

(2)

式中：t₀表示震源发震时刻, x、y、z表示震源坐标, t_k^cal表示第k个台站的模拟波至时间, t_k^obs表示第k个台站的观测波至时间, N为台站个数, p为范数, 通常p取1或2.

2 基于DIRECT算法的震源扫描定位算法 2.1 DIRECT算法

DIRECT算法是一种新的确定性全局最优化算法, 由Jones等人(1993) 首次提出, 该算法的名字来自于“Dividing Rectangles”(分割矩形)的首字母缩写词, 同时也指出它是一种直接(DIRECT)搜索算法.DIRECT算法用于搜索简单约束条件下的多维(多变量)函数的全局最小值, 非常适合于具有确定变量空间的“黑盒”问题的最优化求解.该算法是对标准的Lipschitz优化算法的一种改造, 使其不依赖于Lipschitz常数, 在每次迭代中利用多个Lipschitz常数同时进行搜索, 实现了问题的最优化, 并且DIRECT算法比传统的Lipschitz算法收敛速度更快.

设有如下优化问题：

(3)

式中, x∈Rⁿ为向量x的取值范围.

DIRECT算法解上述最优化问题的基本思想是：首先将变量空间的取值范围归一化, 变成一个n维超立方, 对于一维情况(单变量), 搜索域就变成了[0, 1]区间; 对于二维情况(两个变量), 搜索域就变成了边长为1的正方形; 对于三维情况(三个变量), 搜索域就是单位立方体; 对于多维情况, 搜索域就是一个超立方体(hyper-cube); 变换完成后, 计算变量空间中心处的最优化目标函数值, 随后将超立方体沿着每一维的方向分割成三部分, 形成许多超矩形体(hyper-rectangles), 对所有的超矩形的中心采样, 同时记录最优化目标函数值; 其次, 在此基础上, 依据超矩形体的中心函数值确定出一系列的所谓“潜在最优矩形体”(potentially optimal rectangles), 意即这些“潜在最优矩形体”所包含的区间内可能存在最优化目标函数的全局最优解.这些“潜在最优矩形体”被进一步分割成更小的超矩形体, 同样对它们的中心采样, 重新评估所有函数值得到新的“潜在最优化矩形体”, 依此类推, 随着搜索域的缩小, 函数逐渐趋近于目标函数的全局最小值.

潜在最优矩形体定义如下：

假设已经把单位超立方体分割成了m个超矩形体, c_i为第i个超矩形体的中心, d_i为第i个超矩形体的中心到其顶点的距离.令e为一正常数, 通常为1.0e^-4, f_min为当前目标函数的最优解.若存在某个变化率常数K, 使得超矩形体j满足：

(4)

我们称超矩形体j为潜在最优超矩形体, 该矩形体在下次迭代中将进一步被细分.DIRECT算法先对初始单位超立方中心采样, 在每次迭代中识别潜在最优超矩形集合并将其细分, 这个过程将一直持续, 直至达到预设迭代次数限制.

确定了待分割空间(潜在最优超矩形体), 需要确定空间的分割方法.DIRECT算法定义的分割方法为三等分法, 对于一维问题, 直接将分割空间三等分即可; 对于n维问题, 需将n个分割方向三等分.其基本步骤描述如下：

先沿着最大边长的维度方向将超矩形体三等分, 若最大边长的维度方向的个数大于1, 则：

1) 确定具有最大边长的维度集合I, 令δ为这个最大边长的三分之一;

2) 对点c±δe_i(i∈I)取样并对其进行函数估值, 其中c为矩形中心, e_i为第i个方向的单位矢量;

3) 沿着I中每个维度方向, 将超矩形体分割成三等分, 分割方向根据该维度方向的w_i值, 从小到大依次进行, 其中：

(5)

下面给出多维DIRECT算法的具体步骤：

步骤1：将搜索空间正规化为单位超立方体, 令c_i为超立方体的中心, 并对函数f(c_i)估值.令f_min=f(c₁), m=1, t=0(t为迭代次数);

步骤2：识别潜在最优超矩形体, 并放入集合S;

步骤3：选择任意超矩形体j∈S;

步骤4：按前面所述方法细分超矩形体, 决定从超矩形体j何处位置采样以及如何将超矩形体分割为更小的子超矩形体.更新f_min并令m=m+Δm, 其中Δm为新采样点的个数;

步骤5：令S=S-{j}, 如果S≠φ则转到步骤3;

步骤6：令t=t+1, 若t=T, 即迭代次数达到限制次数, 停止迭代, 否则转到步骤2.

为了说明Direct算法迭代的详细过程.图 2描述了DIRECT算法优化二维测试函数Brain的前三次迭代网格搜索情况.图中第一列表示一次迭代细分矩形的初始状态, 第二列表示识别潜在最优矩形集合(图中带阴影矩形), 第三列表示如何对潜在最优矩形采样及进一步将其细分.

图 3给出了三维情况下的超矩形体分割的基本示意图.对于多维的情况, 按上述空间分割方法在各维划分即可.

2.2 基于DIRECT算法的微震定位方法

基于DIRECT算法的微震定位方法, 就是利用传统震源扫描算法的思想, 将(1) 式或者(2) 式作为目标函数, 在求解最优解时, 利用DIRECT搜索算法, 搜索整个解空间的全局最优解.该方法无需划分网格, 在计算过程中能根据解的位置, 在靠近最优解的地方搜索密度大, 在远离最优解的地方搜索密度小, 减少了搜索次数, 加快了计算效率, 且精度较高.其基本步骤如下：

(1) 构造目标函数, 可用(1) 式或(2) 式作为目标函数;

(2) 确定解空间, 根据实际资料, 确定震源x、y、z和发震时刻t₀的大致区间范围;

(3) 利用DIRECT算法求解目标函数全局最优解, 得到最终的定位结果.

3 数值算例

为测试基于DIRECT算法的微震定位方法的应用效果, 本节设计了地面微震监测模型(模型1) 和井中微震监测模型(模型2) , 基于射线追踪方法合成了微震记录, 用于反演测试.本文主要讨论定位算法准确性及运算效率, 因速度模型复杂与否与该测试相关性不大, 故仅设计了单一速度的均匀模型, 未考虑复杂速度模型下的微震定位问题.

3.1 地面微震监测定位模型

图 4是地面微震监测三维示意图, 地面布置环网状单分量检波器25个, 设计震源10个.图 5是利用射线追踪正演模拟方法合成的一个微震记录(Z分量).利用正演模拟得到的微震记录, 将(1) 式作为最优化目标函数.在速度模型已知的情况下, 微震定位需要求解震源的空间坐标(X、Y、Z)和震源的发震时刻(t₀).将(1) 式作为目标函数时, 若t₀未知, 传统的震源扫描算法将会非常费时, 本次基于模型的对比未将t₀作为反演变量, 仅反演X、Y、Z三个变量.搜索范围限制为X[-50.0, 50.0], Y[-50, 50], Z[300, 400], t₀=0, v=2500.0.

图 6是其中一个震源点A(-2.186, -21.310, 369.740) 的Dierect算法所有网格搜索点的三维示意图, 颜色代表搜索顺序, 一共20007次搜索, 图 7-图 9是所有搜索位置的XY、YZ、XZ平面投影图, 从图中可以看出, 远离震源点的位置搜索密度小, 在震源位置点附近的搜索密度大.表 1是对10个震源位置定位点的定位结果, 表中给出了计算次数及耗时.表 2是对10个震源点利用传统的网格搜索法的定位结果,选用的网格大小是1 m×1 m×1 m; 从表 1和表 2可以看出, DIRECT方法的定位结果在X、Y、Z方向上的误差非常小, 均在10^-2 m级, 而传统网格搜索法的定位结果依赖于网格的尺寸, 由于网格大小是1 m×1 m×1 m, 所以其定位结果均落在网格点上, 要想提高精度, 只能减小网格尺寸.DIRECT方法每个震源点定位计算次数均在20000次左右, 时间在0.3 s左右, 而传统网格搜索定位方法每个震源要循环计算1000000次, 耗时4 s左右.相比传统网格搜索定位方法, DIRECT方法搜索次数提高了约50倍, 时间提高了13～15倍左右.DIRECT算法是一种直接网格搜索算法, 无需划分网格, 定位精度与网格无关, 传统网格搜索方法定位精度取决于网格的尺寸, 本例中如果将网格定义为0.5 m×0.5 m×0.5 m, 则每个震源点定位要循环8000000次, 耗时将更多, 且定位结果只能落在网格点上, 如果震源不在网格点上, 误差更大; 对于DIRECT算法则不存在这个问题.

3.2 井中微震监测定位

本例测试DIRECT搜索定位算法在井中微震监测的定位效果.表 3是井中传感器的坐标, 包括两个观测井, 图 10是对应的三维观测系统图.图 11是利用射线追踪合成的一个震源点的微震记录.仅反演X、Y、Z三个变量, 搜索范围限制为X[-100.0, 100.0], Y[-100, 100], Z[300, 400], t₀=0, v=2500.0.

图 12是其中一个震源点A(-2.186, -21.310, 369.740) 的DIRECT算法所有网格搜索点的三维示意图, 颜色代表搜索顺序, 一共20007次搜索, 虽然搜索范围扩大了, 但是搜索次数并未增加, 图 13～图 15是所有搜索位置的XY、YZ、XZ平面投影图, 从图中可以看出, 远离震源点的位置搜索密度小, 在震源位置点附近的搜索密度大.表 4是对10个震源位置的定位结果, 表中给出了计算次数及耗时.表 5是对10个震源点利用传统的网格搜索法的定位结果, 选用的网格大小是0.5 m×0.5 m ×0.5 m; DIRECT方法每个震源点定位搜索次数均在20000次左右, 时间在0.3 s左右, 而传统网格搜索定位法由于扩大了搜索范围, 减小了网格尺寸, 每个震源点定位要循环8000000次, 耗时6 s左右.相比传统网格搜索定位方法, 本文方法搜索效率提高了近400倍, 时间提高了20倍左右.从表 4和表 5的定位误差来看, 二者基本相当.

3.3 讨论

在实际微震监测中, 反演定位需要解的变量至少包括震源坐标(X、Y、Z)及发震时刻t₀, 但本节所测试的地面监测及井中监测定位方法时, 未将发震时刻作为反演变量传统震源扫描算法与DIRECT震源扫描算法进行测试对比.这是因为若将发震时刻也当做反演变量, 仍然将(1) 作为最优化目标函数, 传统震源扫描算法将耗费大量的时间, 因为在搜索的发震时刻范围内, 每给定一个t₀就要重复所有网格点, 计算目标函数的最优值.若给定的t₀扫描范围为[t₀₁, t₀₂], 扫描步长为dt₀, 则需要重复n=(t₀₂-t₀₁)/dt₀次网格点运算, 时间将会是t₀已知情况的n倍.但是对于DIRECT震源扫描算法, 仅仅是变量空间多了一维, 并不会成倍增加其搜索次数.

4 实例

本文算法已应用在煤层气井水平分段压裂微震监测系统当中, 图 16是对淮北某煤矿煤层气井水平多段压裂微震监测定位效果图.

实例应用效果表明, 本文算法运行稳定, 运算效率高, 定位精度高, 可应用于微震监测实时定位系统中.

5 结论

震源扫描算法是一种经典的地震定位算法, 被广泛应用于地震定位中, 故该方法也是微地震震源定位中的重要理论基础.但是传统的震源扫描算法, 是对解空间按照一定的网格尺寸网格化之后, 遍历所有的网格点, 求解目标函数的全局最优化解, 其缺点是：定位精度依赖于网格划分的尺寸, 计算效率随着网格的细化而降低, 特别是对于多参数优化问题, 计算量会随着尺度的减小而呈几何级数增加.

本文基于传统震源扫描算法, 在搜索过程中利用DIRECT矩形分割搜索算法, 实现了基于DIRECT算法的微震定位方法, 该方法是一种直接搜索算法, 无需划分网格, 定位精度与网格尺寸无关, 其搜索点随着最优解的靠近而加密, 按照DIRECT算法的搜索路径, 实现了全局最优解.

分别对地面与井中微震监测模型进行定位测试, 试验结果表明DIRECT算法是一种快速、准确的定位算法, 与传统的网格搜索算法对比, DIRECT算法搜索次数少, 计算耗时短, 计算精度高, 尤其适合大空间内搜索微震震源位置, 为水力压裂微震监测过程中大量数据的实时、快速、准确定位提供了可能.

致谢 感谢审稿专家提出的宝贵修改意见和编辑部的大力支持！