CN112257241B - 一种三角网菲涅尔带时差层析反演方法 - Google Patents

Abstract

本发明公开了一种三角网菲涅尔带时差层析反演方法，本发明使用的层析反演方法不易受环境因素影响，工作难度小，施工成本较低；并且不同于传统层析反演使用矩形网格来剖分模型，本发明能够较好地贴合边缘不规则的桥梁模型，减少了正演的误差，使得后续反演过程也能取得良好的效果。本发明使用离炮点最近的道作为参考道，计算其理论初至时间与实际测量的初至时间之差，并以此为基准推算出所有检波点的真实初至时间。利用该方法可以较准确地拾取初至走时，有利于后续反演。本发明使用的菲涅尔带技术避免了传统层析反演需要大量的检波设备来确保精确度的问题，使得在桥梁检测这类相对小型的项目中使用层析反演方法依旧能取得较好的桥梁检测效果。

Description

一种三角网菲涅尔带时差层析反演方法

技术领域

本发明涉及桥梁检测领域，尤其涉及一种三角网菲涅尔带时差层析反演方法。

背景技术

传统的桥梁检测方法有大地测量法、物理学传感器法或光纤传感器法等方法，这些方法或易受环境影响，或工作难度大。层析反演作为石油地震勘探中的一种反演地表结构的常用方法，使用层析反演技术进行桥梁检测不仅现场操作简单，而且精确度较高。然而石油勘探中使用的常规层析反演技术是基于矩形网格的大尺度反演，而桥梁模型具有复杂的外形轮廓且尺度较小，使用规则的矩形网格来建立反演模型，由于不能精确贴合桥梁模型，因而会造成不可忽略的误差。

在桥梁检测领域，回弹检测法，动力试验法以及荷载试验法都是常用的桥梁检测技术，然而，以上方法都不能很好反映桥梁的内部结构与状态信息。现有技术中袁浩,唐英等人采用激振点与受信点之间两两对穿形成弹性波射线，将弹性波射线编织成的射线网作为检测剖面，反演剖面波速图来检测故障区域，其使用的层析反演是通过观测数据获取桥梁内部结构信息的核心手段，但采用的是以数学射线来处理地震波的传播，这种以抽象的高频近似不能反映真实的物理过程，且它的射线分布不均匀。本发明引入菲涅尔带，一段菲涅尔带内，不同部分对地震波能量的作用是不均等的，一般离菲涅尔带轴心越远，能流密度越小，更能真实反映波在桥梁内部传播速度的衰减。此外，由于桥梁结构复杂，造型多样，边缘不规则，地学上常采用的矩形网划分不能很好的贴合桥梁模型边缘，且在射线分布稀疏情况下的反演精度较低，计算中产生不可忽略的离散化误差。此外，由于检测设备无法精确测量振动的触发时间，每次记录的振动触发时间都有不同程度的延时误差，使得同一次振动到达各检波点的初至时间都有相同的误差。本发明采用时差法，利用理论走时与实际测量走时之差计算出真实初至时间，并通过迭代反演，获得理想的桥梁速度模型。

发明内容

本发明一种三角网菲涅尔带时差层析反演方法，包括以下步骤：

S1：建立桥梁模型，对桥梁模型进行三角剖分，在模型中布设炮点和检波点；

S2：给定初始模型慢度S0；

S3：设定总迭代次数N，设定反演终止的平均误差门槛值ξ₀；

S4：将三角网格的累计慢度修改量设为0；

S5：遍历每个炮点开始反演计算；

S6：通过正演计算得到第i个炮点到每个检波点的旅行时表和射线轨迹,使用时差法对拾取的初至波的旅行时进行计算，得到合成后的旅行时；

S7：找出每根射线菲涅尔带范围内的三角网格，加入表L中；

S8：计算表L中每个三角网格的慢度修改量，并累加慢度修改量；

S9：重复步骤S5-S8，完成炮点循环，使用累计慢度修改量修改模型；

S10：利用新的慢度模型，计算所有炮点到检波点的旅行时t′_ij；

S11：使用评价函数

其中n_S为炮点个数，n_R为每炮影响到的检波点个数，t_ij为计算合成后各道的真实初至旅行时，t′_ij为正演的初至旅行时；

S12：重复步骤S4-S11，直到满足ξ＜ξ₀或迭代次数大于N，退出迭代，反演计算完成。

进一步，所述步骤S6中的射线为地震波从炮点出发到达检波点的轨迹。所述步骤S7中的菲涅尔带是由地震波沿中心射线一定范围内呈条带状传播的射线构成。所述菲涅尔带范围判定参数为α＝1.0/(2.0*f)，其中f为振动波的频率，与菲涅尔带范围的大小成反比，由于桥梁体积相对较小，使用高频2000HZ进行测量效果较为理想，通过将三角网格的旅行时与射线的旅行时的差与α进行比较来判断三角网格是否在当前射线的菲涅尔带范围内。所述步骤S6中的时差法，是由于设备对振动起始时间记录不准确，导致记录的振动到达各检波点的初至时间整体上都有偏差。选出距离起振点最近的道作为参考道，记初至时间为t₀，理论初至时间为T₀，计算各道的合成初至时t_j，t_j＝T₀+Δt_j，其中Δt_j为拾取的各道初至时与参考道初至时间t₀的差值，j为检波道，j∈[1,n]，n为起振点的检波道数。

进一步，所述步骤S7-S8具体包括以下步骤：

S7-1：要找出射线的菲涅尔带范围内的所有三角网格，先要找出当前射线通过的三角网格，每条射线由多个节点连接而成，为此先建立一条队列Q；

S7-2：遍历射线的每个节点，得到节点所在的三角网格，将其加入队列Q，并标记为菲涅尔带范围内的三角网格，下述简称为加上flag标记，每条射线的flag值不同，并加入表L当中，这些射线直接通过的三角网格的权值Ω＝1；

S7-3：计算射线的理论初至时间t′_ij与真实初至时间t_ij的差Δt＝t′_ij-t_ij；

S7-4：取出队列Q中的第一个三角网格，找出其邻接三角网格，若果该邻接网格没有flag标记，则计算该邻接网格的顶点到当前射线的炮点和检波点的时间之和t，将三个顶点的t的平均值at记为该邻接三角网格的初至时间，计算at与射线理论初至时间t′_ij的差的绝对值T_f；如果T_f<α，说明邻接三角网格在菲涅尔带的范围内，则将邻接三角网格加入队列Q，并加入表L中，在这一步中计算的三角网格均射线菲涅尔带范围内，但未被射线直接通过，它们的权值Ω＝1.0-2.0*T_f*f；

S7-5：判断队列Q是否为空，如果不为空，重复步骤S8-4，如果队列Q为空，说明当前射线菲涅尔带范围内的三角网格已经被全部找出并加入表L中，开始下一步计算；

S7-6：遍历表L中的每个三角网格，将每个网格的权值Ω除以表L中所有三角网格的总权值，得到当前三角网格的权重百分比Ω，则当前三角网格的慢度修改量ΔS＝Ω*Δt；

S7-7：表L的所有三角网格计算完毕后，当前射线便计算完成。

进一步，所述拾取的初至波由数字工程地震仪测量获取，所述参考道为最先接收到初至波的检波道，其初至时记为t₀。所述各道拾取初至时t_j由数字工程地震仪测量并提取得到，其中j为炮点发出的振动影响到的检波点数，即每炮的检波道数。

本发明的有益效果是：

本发明是一种，与现有技术相比，本发明使用的层析反演方法不易受环境因素影响，工作难度小，施工成本较低；并且不同于传统层析反演使用矩形网格来剖分模型，本发明使用了Delaunay三角网剖分技术对桥梁模型进行剖分，能够较好地贴合边缘不规则的桥梁模型，减少了正演的误差，使得后续反演过程也能取得良好的效果，此外，本发明使用的菲涅尔带技术避免了传统层析反演需要大量的检波设备来确保精确度的问题，使得在桥梁检测这类相对小型的项目中使用层析反演方法依旧能取得较好的桥梁检测效果。

附图说明

图1为对桥梁模型进行三角剖分后的效果图；

图2为在桥梁模型上布设测点的效果图，其中空心点为炮点，实心点为检波点；

图3为一根射线的菲涅尔带样式图，其中椭圆内线条为射线，环绕射线的椭圆曲线为菲涅尔带；

图4为三角形网格旅行时at的计算示意图，at＝t_a+t_b；

图5为射线的菲涅尔带范围内的三角网格慢度修改量计算流程图；

图6为正演效果图，红色线为射线追踪路径；

图7为初至拾取示意图，即计算时差Δt并得到真实初至时间t_ij。

具体实施方式

下面结合附图以及具体实施例对本发明作进一步描述，在此发明的示意性实施例以及说明用来解释本发明，但并不作为对本发明的限定。

如图1-7所示：本发明包括以下步骤：

S1：建立桥梁模型，对桥梁模型进行三角剖分，在模型中布设炮点和检波点；

S2：给定初始模型慢度S0；

S3：设定总迭代次数N，设定反演终止的平均误差门槛值ξ₀；

S4：将三角网格的累计慢度修改量设为0；

S5：遍历每个炮点开始反演计算；

S6：通过正演计算得到第i个炮点到每个检波点的旅行时表和射线轨迹,使用时差法对拾取的初至波的旅行时进行计算，得到合成后的旅行时；

S7：找出每根射线菲涅尔带范围内的三角网格，加入表L中；

S8：计算表L中每个三角网格的慢度修改量，并累加慢度修改量；

S9：重复步骤S5-S8，完成炮点循环，使用累计慢度修改量修改模型；

S10：利用新的慢度模型，计算所有炮点到检波点的旅行时t′_ij；

S11：使用评价函数

其中n_S为炮点个数，n_R为每炮影响到的检波点个数，t_ij为计算合成后各道的真实初至旅行时，t′_ij为正演的初至旅行时；

S12：重复步骤S4-S11，直到满足ξ＜ξ₀或迭代次数大于N，退出迭代，反演计算完成。

进一步，所述步骤S6中的射线为地震波从炮点出发到达检波点的轨迹。所述步骤S7中的菲涅尔带是由地震波沿中心射线一定范围内呈条带状传播的射线构成。所述菲涅尔带范围判定参数为α＝1.0/(2.0*f)，其中f为振动波的频率，与菲涅尔带范围的大小成反比，由于桥梁体积相对较小，使用高频2000HZ进行测量效果较为理想，通过将三角网格的旅行时与射线的旅行时的差与α进行比较来判断三角网格是否在当前射线的菲涅尔带范围内。所述步骤S6中的时差法，是由于设备对振动起始时间记录不准确，导致记录的振动到达各检波点的初至时间整体上都有偏差。选出距离起振点最近的道作为参考道，记初至时间为t₀，理论初至时间为T₀，计算各道的合成初至时t_j，t_j＝T₀+Δt_j，其中Δt_j为拾取的各道初至时与参考道初至时间t₀的差值，j为检波道，j∈[1,n]，n为起振点的检波道数。

进一步，所述步骤S7-S8具体包括以下步骤：

S7-1：要找出射线的菲涅尔带范围内的所有三角网格，先要找出当前射线通过的三角网格，每条射线由多个节点连接而成，为此先建立一条队列Q；

S7-2：遍历射线的每个节点，得到节点所在的三角网格，将其加入队列Q，并标记为菲涅尔带范围内的三角网格，下述简称为加上flag标记，每条射线的flag值不同，并加入表L当中，这些射线直接通过的三角网格的权值Ω＝1；

S7-3：计算射线的理论初至时间t′_ij与真实初至时间t_ij的差Δt＝t′_ij-t_ij；

S7-4：取出队列Q中的第一个三角网格，找出其邻接三角网格，若果该邻接网格没有flag标记，则计算该邻接网格的顶点到当前射线的炮点和检波点的时间之和t，将三个顶点的t的平均值at记为该邻接三角网格的初至时间，计算at与射线理论初至时间t′_ij的差的绝对值T_f；如果T_f<α，说明邻接三角网格在菲涅尔带的范围内，则将邻接三角网格加入队列Q，并加入表L中，在这一步中计算的三角网格均射线菲涅尔带范围内，但未被射线直接通过，它们的权值Ω＝1.0-2.0*T_f*f；

S7-5：判断队列Q是否为空，如果不为空，重复步骤S8-4，如果队列Q为空，说明当前射线菲涅尔带范围内的三角网格已经被全部找出并加入表L中，开始下一步计算；

S7-6：遍历表L中的每个三角网格，将每个网格的权值Ω除以表L中所有三角网格的总权值，得到当前三角网格的权重百分比Ω，则当前三角网格的慢度修改量ΔS＝Ω*Δt；

S7-7：表L的所有三角网格计算完毕后，当前射线便计算完成。

进一步，所述拾取的初至波由数字工程地震仪测量获取，所述参考道为最先接收到初至波的检波道，其初至时记为t₀。所述各道拾取初至时t_j由数字工程地震仪测量并提取得到，其中j为炮点发出的振动影响到的检波点数，即每炮的检波道数。

石油勘探中使用的常规层析反演技术是基于矩形网格的大尺度反演，而桥梁模型具有复杂的外形轮廓且尺度较小，使用规则的矩形网格来建立反演模型，由于不能精确贴合桥梁模型，因而会造成不可忽略的误差。为此，本发明使用了自适应三角网剖分技术对桥梁模型进行剖分，能较好地贴合桥梁形状，从而有效地减少反演计算产生的离散化误差。

此外，由于检测设备存在不可避免的延时，导致每次记录的振动到达各检波点的初至时间都有整体偏差。本发明计算离炮点最近的几个检波点的理论初至时间与实际测量的初至时间之差，并以此为基准推算出所有检波点的真实初至时间。利用该方法可以较准确地拾取初至走时，有利于后续反演。

最后，由于桥梁检测受成本或其他因素的制约，接收设备有限，因此相对于地震勘探，桥梁检测时的射线密度较为稀疏，为此本发明使用了菲涅尔带层析反演方法，使得位于射线中心及其菲涅尔带范围内的三角网格均参与反演计算，因此在射线相对稀疏的情况下也可以取得较好的反演效果，从而降低施工成本。

本发明为了克服桥梁检测领域现有技术不能在真实反映桥梁内部波传播状态、使用规则网络对复杂结构桥梁剖分存在较大误差的缺陷以及传感器无法精确记录振动的触发起始时间的缺陷，本发明提供了一种三角网菲涅尔带时差层析反演方法，采用Delaunay三角网对桥梁模型进行剖分，减少了规则网格模型正演计算时在模型边界产生的较大误差，引入菲涅尔带反演技术可增大射线覆盖范围，改善射线覆盖不均匀的程度，降低反演矩阵的稀疏性，最后根据射线理论走时与实际测量的时间差来修正射线初至时间，并根据该时间采用迭代求解来修正桥梁速度模型，尽可能真实地反映桥梁内部波速度情况，在不用求解射线精确旅行时、射线相对稀疏的情况和低施工成本下也可以取得较好的反演效果。

本发明的技术方案不限于上述具体实施例的限制，凡是根据本发明的技术方案做出的技术变形，均落入本发明的保护范围之内。

Claims (8)

1.一种三角网菲涅尔带时差层析反演方法，其特征在于：包括以下步骤：

S1：建立桥梁模型，对桥梁模型进行三角剖分，在模型中布设炮点和检波点；

S2：给定初始模型慢度S0；

S3：设定总迭代次数N，设定反演终止的平均误差门槛值ξ₀；

S4：将三角网格的累计慢度修改量设为0；

S5：遍历每个炮点开始反演计算；

S6：通过正演计算得到第i个炮点到每个检波点的旅行时表和射线轨迹,使用时差法对拾取的初至波的旅行时进行计算，得到合成后的旅行时；

S7：找出每根射线菲涅尔带范围内的三角网格，加入表L中；

S8：计算表L中每个三角网格的慢度修改量，并累加慢度修改量；

S9：重复步骤S5-S8，完成炮点循环，使用累计慢度修改量修改模型；

S10：利用新的慢度模型，计算所有炮点到检波点的旅行时t′_ij；

S11：使用评价函数

其中n_S为炮点个数，n_R为每炮影响到的检波点个数，t_ij为计算合成后各道的真实初至旅行时，t′_ij为正演的初至旅行时；

S12：重复步骤S4-S11，直到满足ξ＜ξ₀或迭代次数大于N，退出迭代，反演计算完成。

2.根据权利要求1所述的三角网菲涅尔带时差层析反演方法，其特征在于：所述步骤S6中的射线为地震波从炮点出发到达检波点的轨迹。

3.根据权利要求1所述的三角网菲涅尔带时差层析反演方法，其特征在于：所述步骤S7中的菲涅尔带是由地震波沿中心射线一定范围内呈条带状传播的射线构成。

4.根据权利要求3所述的三角网菲涅尔带时差层析反演方法，其特征在于：所述菲涅尔带范围判定参数为α＝1.0/(2.0*f)，其中f为振动波的频率，与菲涅尔带范围的大小成反比，由于桥梁体积相对较小，使用高频2000HZ进行测量效果较为理想，通过将三角网格的旅行时与射线的旅行时的差与α进行比较来判断三角网格是否在当前射线的菲涅尔带范围内。

5.根据权利要求1所述的三角网菲涅尔带时差层析反演方法，其特征在于：所述步骤S6中的时差法，是由于设备对振动起始时间记录不准确，导致记录的振动到达各检波点的初至时间整体上都有偏差；选出距离起振点最近的道作为参考道，记初至时间为t₀，理论初至时间为T₀，计算各道的合成初至时t_j，t_j＝T₀+Δt_j，其中Δt_j为拾取的各道初至时与参考道初至时间t₀的差值，j为检波道，j∈[1,n]，n为起振点的检波道数。

6.根据权利要求4所述的三角网菲涅尔带时差层析反演方法，其特征在于：所述步骤S7-S8具体包括以下步骤：

S7-1：要找出射线的菲涅尔带范围内的所有三角网格，先要找出当前射线通过的三角网格，每条射线由多个节点连接而成，为此先建立一条队列Q；

S7-2：遍历射线的每个节点，得到节点所在的三角网格，将其加入队列Q，并标记为菲涅尔带范围内的三角网格，下述简称为加上flag标记，每条射线的flag值不同，并加入表L当中，这些射线直接通过的三角网格的权值Ω＝1；

S7-3：计算射线的理论初至时间t′_ij与真实初至时间t_ij的差Δt＝t′_ij-t_ij；

S7-4：取出队列Q中的第一个三角网格，找出其邻接三角网格，若果该邻接三角网格没有flag标记，则计算该邻接三角网格的顶点到当前射线的炮点和检波点的时间之和t，将三个顶点的t的平均值at记为该邻接三角网格的初至时间，计算at与射线理论初至时间t′_ij的差的绝对值T_f；如果T_f<α，说明邻接三角网格在菲涅尔带的范围内，则将邻接三角网格加入队列Q，并加入表L中，在这一步中计算的三角网格均射线菲涅尔带范围内，但未被射线直接通过，它们的权值Ω＝1.0-2.0*T_f*f；

S7-5：判断队列Q是否为空，如果不为空，重复步骤S8-4，如果队列Q为空，说明当前射线菲涅尔带范围内的三角网格已经被全部找出并加入表L中，开始下一步计算；

S7-6：遍历表L中的每个三角网格，将每个网格的权值Ω除以表L中所有三角网格的总权值，得到当前三角网格的权重百分比Ω，则当前三角网格的慢度修改量ΔS＝Ω*Δt；

S7-7：表L的所有三角网格计算完毕后，当前射线便计算完成。

7.根据权利要求5所述的三角网菲涅尔带时差层析反演方法，其特征在于：所述拾取的初至波由数字工程地震仪测量获取，所述参考道为最先接收到初至波的检波道，其初至时记为t₀。

8.根据权利要求5所述的三角网菲涅尔带时差层析反演方法，其特征在于：所述各道拾取初至时t_j由数字工程地震仪测量并提取得到，其中j为炮点发出的振动影响到的检波点数，即每炮的检波道数。

Images