堤防隐患探测的地震层析成像应用研究

Application of Seismic Tomography in Detecting Hidden Hazards of Embankments

摘要: 简要介绍了地球物理反演中图像重建方法的基本原理、SIRT算法的原理和实现步骤，在现有地震层析成像的方法上研究堤防层析成像。本研究利用MATLAB软件进行了堤防地震层析成像的数值模拟，并分析程序的反演效果，并做了实际工程验证。结果表明SIRT算法取得了良好的反演效果，验证了地震层析成像探测堤防隐患的可行性。

Abstract: This paper briefly introduces the basic principle of image reconstruction method in geophysical inversion, the principle and implementation steps of SIRT algorithm, and studies levee tomography based on the existing seismic tomography methods. In this paper, the numerical simulation of seismic tomography of dike is carried out with MATLAB software, and the inversion effect of the program is analyzed, and the practical engineering verification is done. The results show that the SIRT algorithm has obtained a good inversion effect and verified the feasibility of seismic tomography to detect hidden dangers of levees.

1. 引言

在防洪减灾领域，堤防作为关键的基础设施，对于保障河流及沿海区域居民的生命和财产安全具有至关重要的作用。堤防的构建与维护质量直接决定了其防洪能力及灾害减轻的效果。鉴于全球气候变化所引发的极端天气事件频发，洪水灾害的发生率也随之上升，这对堤防的安全性提出了更为严格的要求。因此，研究和发展高效的堤防隐患探测技术，及时发现和处置隐患，对于提高防洪减灾能力和保障社会可持续发展具有重要意义[1]。

20世纪80年代以来，我国开始采用物探方法[2]开展堤防隐患探测的研究工作，如高密度电法[3]、地质雷达法[4]、瞬变电磁法等[5]-[7]。物探方法在堤防质量检测中起到了重要的作用。

2. 基本理论

2.1. 离散图像重建原理

在层析成像技术中，逆问题求解是关键步骤，它涉及从一系列测量的投影数据中恢复物体内部特定物理量的空间分布。

该过程正是Radon逆变换理论的实践应用。1917年，奥地利数学家约翰·雷当提出了一个重要的理论，即图像函数可以通过其Radon变换的无限集合$\left(Rf\right)$ 被唯一且精确地重建，为层析成像技术提供了坚实的数学基础。在二维连续图像的数学描述中，图像函数是关键。该图像函数的Radon变换，表示为$Rf$ ，$Rf$ 是指$f$ 沿特定直线方向的线积分，这个积分值也被称作投影值。

$\left[Rf\right]\left(l,\theta \right)={\displaystyle \underset{L}{\int }f\left(r,\phi \right)}ds\left(\left|l\right|\le E,0\le \theta \le \pi \right)$ (1)

在此，$l$ 和$\theta$ 代表直线的空间参数，其中$l$ 表示从原点到直线的垂直距离，即直线的法线长度，$\theta$ 表示法线与参考轴之间的夹角。

在极坐标系统$\left(r,\phi \right)$ 中，假定成像区域为$r<a,\phi <2\pi$ ，该区域被划分成$J$ 个不相交的像元。对于每一个像元，可以定义其对应的基函数为：

$b_j=\{\begin{array}{l}1\left(r,\varphi \right)在第j像元\\ 0其他\end{array}$ (2)

则离散图像重建问题可以转化为在给定的投影数据$b$ 条件下，对图像向量$f$ 进行估计的问题。

$b=Af+e$ (3)

这里提到的$A$ 被称为投影矩阵。解决这一问题的一些常见方法包括代数重建技术(ART)和联合迭代重建技术(SIRT)等。

2.2. SIRT算法

迭代重建技术(SIRT)是在代数重建技术(ART)的基础上发展而来的，最初由Gilbert提出。SIRT算法的核心思想与ART的主要区别在于，SIRT不是直接使用修改值$f_j^{\left(q,i+1\right)}$ 来调整投影数据，而是将每一轮迭代中所有射线产生的修改值累积起来，并存储在一个特定的空间内。当一轮迭代完成后，SIRT算法会计算这些修改值的某种平均值，然后应用这些平均值$\Delta \overline{f}$ 来进行图像的更新。这种方法使得SIRT在处理射线数据时更为精确和有效。

$f^{\left(q+1\right)}=f^{\left(q\right)}+\Delta \overline{f}$ (4)

对每个像素点的衰减系数进行调整，并保留这些调整值以供后续迭代利用。在最基本的SIRT迭代修正模型中，通过合理选择平滑矩阵以及设定一个恒定的松弛参数µ，可以得到如下公式：

$f_j^{k+1}=f_j^k+\frac{{\displaystyle \sum _{i=1}^I\left[a_{ij}\left({\tau }_i-\left(a_i,f^k\right)\right)/{\displaystyle \sum _{j=1}^J{a}_{ij}}\right]}}{\mu +{\displaystyle \sum _{i=1}^I{a}_{ij}}}$ (5)

对于编号为$j=1,2,\cdot \cdot \cdot ,J$ 的各个像素单元，以及编号为$I=1,2,\cdot \cdot \cdot ,I$ 的射线，系数矩阵A表示为$A=\left(a_1,a_2,\cdot \cdot \cdot a_I\right)$ ，代表迭代次数。根据公式(5)，可以通过计算通过像素单元j的所有射线的修改增量的加权平均值，来确定衰减系数的调整增量。

尽管SIRT算法存在对计算机内存要求较高的不足，但其良好的收敛特性和较高的图像重建精度是其显著优势。此外，SIRT方法在穿透孔径CT (cross-hole CT)领域也得到了广泛应用。本研究采纳了SIRT算法来构建堤防CT程序。

3. 程序介绍

该程序使用MATLAB软件进行编写，在本研究中采用了SIRT (Simultaneous Iterative Reconstruction Technique)算法，程序主要分为两个部分，分别为堤防的数值剖切面(梯形界面)和水平剖切面(矩形界面)。

本研究提出的梯形模型软件，旨在通过构建堤防的梯形几何结构模型，实现对堤防结构的地震波传播特性的全面分析。程序的主要功能涵盖了模型的建立、正演模拟、反演分析等关键步骤。在模型建立阶段，程序能够精确地模拟堤防的几何形态，为后续的地震波传播模拟提供基础。正演模拟过程涉及地震波在堤防模型中的传播，通过模拟波的传播路径和速度，为反演分析提供理论基础。反演分析则利用观测数据，推断堤防内部结构的特性，进而生成地震波在堤防中的速度分布图。该速度分布图是识别堤防隐患位置的关键，它不仅揭示了地震波在堤防中的传播速度变化，而且通过可视化手段，直观地展现了潜在的隐患区域。最终，程序能够准确地定位堤防中的隐患，显著提升了堤防隐患探测的效率，并增强了数据解释的可视化程度，为堤防的安全评估和管理提供了强有力的技术支持。

4. 数值模拟

在对堤防工程进行地震层析成像探测其内部隐患时，研究通常从水平和垂直两个维度进行。尽管已有学者对垂直方向的层析成像进行了深入探究，本研究专注于水平方向的层析成像技术。在堤防结构中，异常区域与堤防填筑料之间的平均速度差异对成像质量具有显著影响，速度差异越大，成像质量越佳。地震波在穿越堤防时，若遇到裂缝、空洞等异常结构，其传播速度会降低，这一现象对于识别堤防内部的异常结构至关重要。

第一步，建立一个网格化的堤防矩形几何模型，并在此模型中确定激发点和采集点的位置；第二步，设定假定的背景探测速度、实际探测速度，通过计算射线穿过每个像素单元的长度，进而推算出每条射线的实际传播时间与理论传播时间；第三步，利用SIRT算法对这些数据进行反演处理，直至达到预定的迭代停止标准；第四步，将反演结果进行可视化处理，绘制出探测区域的速度分布图。通过分析速度分布图，可以识别出隐患异常体的具体位置。这一流程不仅提高了堤防隐患探测的效率，而且通过数据的可视化，增强了对隐患位置的识别准确性。

基于MATLAB软件建立模型，并设置低速异常体，模型参数见表1：

Table 1. Model parameters table

表1. 模型参数表

模型尺寸(m)	网格划分	正常区域波速(m/s)	异常区域波速(m/s)	异常体数量(个)	震源数量 (个)	检波器数量(个)
10 × 12	1 × 1	300	200	2	12	12

堤防模型的左侧界定为迎水侧，而右侧则为背水侧。激发点均匀布置在模型的左侧，采集点则位于模型的右侧，两者之间的间隔设定为1米。重建区域的网格划分为120个。在坐标位置(2, 6)和(5, 5)处，特别设置了两个低波速异常区域，其波速被设定为200米/秒，而模型中其余区域的波速则被设定为300米/秒。模型详情见图1(a)。

本文采用SIRT算法进行计算，反演结果见图1(b)。

Figure 1. Includes low-velocity volume model and inversion diagram

图1. 含低速体模型及反演图

地震层析成像技术通过对观测到的地震波的运动学(如走时、射线路径)和动力学(如波形、振幅)特性进行分析，利用反演方法来推断地下介质的速度结构及其他物性参数。这种方法依赖于对地震射线的分布和精度的分析，成像结果的分辨率与研究区射线分布情况密切相关。在图像重建领域，精确评估重建结果的准确性是极为关键的。本研究引入了三种图像差异度量指标：归一化均方根误差(d)、归一化平均绝对误差(r)和最大误差(e)。这些指标为量化重建图像与真实图像之间的差异提供了方法。具体来说，归一化均方根误差(d)反映了重建图像与真实图像之间的均方根误差，d值接近0表示重建图像能够精确复制真实图像，而较高的d值则暗示两者之间存在较大差异。归一化平均绝对误差(r)衡量了重建图像与真实图像之间的平均绝对误差，r值接近0意味着误差较小，而r值的增加则显示误差的增加。最大误差(e)则代表了原始图像与重建图像之间的最大平均衰减差异，它是评估重建质量的一个严格标准。

Table 2. Image measurement distance values of inversion results

表2. 反演结果的图像测量距离值

迭代次数	d	r	e
30	0.5545	0.0078	0.00082
60	0.5008	0.0065	0.00071
90	0.5007	0.0063	0.00068

迭代次数越多，经过每一个像元的射线条数越多，反演结果就越准确，从表2中可以看出，迭代次数为90次时，SIRT算法的反演结果趋于稳定，可以看作已收敛。由图1可以看出，反演程序准确地反演出了堤防隐患的位置。

5. 工程应用

选择某地黄河大堤裂缝处进行实地探测，仪器选用WZG-24A工程地震仪，配套检波器12组，使用锤击震源。在迎水面布设震源，两相邻震源点水平距离为1 m，一次测量共布设12个震源点。背水面布设检波器，两相邻检波器水平距离为1 m，一次测量共布设12个检波器。从大堤顶部开始布设测线，测线1探测完成后，将震源和检波器沿坡面下移0.5 m，继续测量测线2，以此类推。反演结果见图2，裂缝发育深度小于0.5 m，未深入堤防内部，后续应尽快采取措施修复。

Figure 2. Defect inversion map of Yellow River levee

图2. 黄河大堤缺陷反演图

6. 结语

本研究基于SIRT算法对堤防隐患的地震层析成像做了应用研究，主要结论包括：

1) SIRT算法的有效性：通过SIRT算法进行的地震层析成像在堤防隐患探测中取得了良好的反演效果，验证了该算法在探测堤防隐患中的可行性。

2) 地震层析成像技术的应用价值：地震层析成像技术在堤防隐患探测领域的应用显著提升了探测工作的效率和准确性，并有效降低了钻探成本。

3) 复杂路径下的探测能力：即使地震波在传播过程中的实际路径并非直线，通过分析速度场内相邻区域的差异，仍能辨识出堤防内部异常体的空间位置和形态特征。

这些结论表明，地震层析成像技术结合SIRT算法在堤防隐患探测中具有重要的应用潜力和实际价值。

地震层析成像技术在堤防隐患探测领域的应用，极大地提升了探测工作的效率和准确性，并且有效降低了钻探成本。在实际操作过程中，若能在观测系统的设计、数据的采集以及后续的数据处理这三个关键环节中选用恰当的技术方法，那么堤防层析成像技术就能够实现更为出色的探测成效。这种技术的有效运用不仅增强了隐患探测的科学性和精确性，而且为堤防的安全监测和管理提供了重要的技术支撑。

基金项目

黄河水利科学研究院科技发展基金(202214)。

NOTES

^*第一作者。

^#通讯作者。