使用多步比较R方法优化二维电阻率测量线的数组 姚 曲（Yao Qu） 林 凯德（Lin Caide） 刘 海（Liu Hai） 刘 向泰（Liu Xiangtai） 孟 徐（Xu Meng） 张 上阳（Zhang Shangyang） 尹 心欣（Yin Zixin） 胡 和松（Hu Hesong）

《Geosciences》：Optimized Arrays for 2-D Resistivity Survey Lines Using a Multi-Step Compare R Method Yao Qu, Caide Lin, Hai Liu, Xiangtai Liu, Xu Meng, Shangyang Zhang, Zixin Yin and Hesong Hu

　　摘要 电阻率层析成像（ERT）的成像质量在很大程度上取决于电极阵列的配置。尽管对称优化的“Compare R”（CR）方法提高了计算效率，但将搜索范围限制在对称数据集内本质上会限制成像精度。为了解决这一限制，本文提出了一种多步骤优化的CR方法，该方法逐步探索对称和非对称阵列，以扩展搜索空间并进一步提高成像精度。数值实验表明，多步骤优化的阵列获得了最高的平均相对模型分辨率（0.646）和结构相似性指数（0.668），超过了对称优化阵列（分别为0.615和0.630）。在管道检测的实地实验中证实，所提出的阵列能够准确识别地下异常的位置和几何形状，并实现了更高的成像精度。在喀斯特洞穴探索中的应用进一步验证了所提出的阵列能够有效检测到深部喀斯特洞穴和基岩界面，这一点通过钻孔得到了验证。此外，还评估了两种优化阵列在不同深度下的检测性能。结果表明，多步骤优化的阵列在更深的区域更可靠地保留了异常的几何形状和电阻率，这归因于深部区域非对称数据点的积累，从而显著提高了数据密度。

1. 引言
电阻率层析成像（ERT）是一种广泛用于成像地下结构的地球物理技术[1,2]。然而，地下结构的成像质量受到复杂结构条件[3,4]和所选电极配置的强烈影响。标准阵列（如Wenner阵列和偶极-偶极阵列）受到固定配置的限制，无法充分利用仪器的多通道能力，特别是在成像复杂地下结构时[5,6,7,8]。全面的数据集包含丰富的地下信息，并通过大幅增加数据密度来改善反演结果[9,10]。然而，在大多数实地调查中获取如此庞大的全面数据集既耗时又不切实际。因此，开发优化算法至关重要，这些算法能够在使用相同数量的数据点的情况下选择最具信息量的电极配置，从而在保持调查效率的同时增强深部成像能力。
大多数现有的阵列优化算法要么专注于最大化灵敏度矩阵元素的总和[11,12,13]，要么专注于最大化模型分辨率矩阵对角元素的总和[14,15,16]。Stummer等人[9]提出了Goodness Function（GF）优化算法，该算法使用排名函数根据灵敏度矩阵和模型分辨率矩阵来评估阵列性能。数值和实地测试表明，优化后的阵列比标准阵列提供了更多的地下信息。为了进一步提高成像质量，Wilkinson等人[17]引入了“Compare R”（CR）方法。尽管CR算法产生的阵列具有最高的平均模型分辨率，但它在整个优化过程中搜索整个配置空间，因此需要比GF算法更多的计算时间。
为了提高CR方法的计算效率，以往的研究集中在四个关键方向上：利用2-D测量线和网格模型的对称性、利用GPU加速、优化Sherman-Morrison更新以及修改全面数据集。Loke等人[16]通过引入模型对称性约束并基于模型分辨率增量矩阵推导出修改后的Sherman-Morrison Rank-1更新公式，从而减少了全面数据集的大小。他们还利用了GPU并行架构，显著减少了计算时间和内存消耗。随后，Loke等人[18]将Sherman-Morrison Rank-1更新从矩阵-向量形式扩展到矩阵-矩阵形式，实现了多个阵列的批量处理。通过结合GPU加速策略，阵列优化算法的计算效率提高了大约两个数量级。最近，Loke等人[19]进一步改进了CR更新过程，将Sherman-Morrison Rank-1更新公式的四电极配置形式替换为两电极对的线性组合形式。通过将优化限制在对称数据集内，对称优化的CR方法在计算成本上又降低了一个数量级。
然而，对称优化的CR方法在可达到的模型分辨率上仍然存在上限。对称数据集仅占全面数据集的约1%到3%，其可达到的模型分辨率大约是全面数据集提供的75%到80%，从而阻碍了优化达到全局最优[19]。此外，在存在非对称地形或具有横向不均匀电阻率分布的复杂地质结构的情况下，从对称优化阵列获得的图像容易发生非对称失真[19,20]。此外，相对较少的实地研究考察了优化阵列在复杂环境中的实际性能[21]。例如，Abdullah等人[22]使用噪声加权的“Compare R”（NWCR）优化阵列，在划分弱区和基岩界面方面表现出优于标准阵列的性能。
为了解决这些限制，本文提出了一种多步骤优化的CR方法，该方法通过引入多步骤策略扩展了对称优化的CR方法。在优化的早期阶段，该方法使用对称数据集以保持高计算效率并建立稳定的数据集；而在后期阶段，引入有限数量的非对称数据集以克服搜索空间的限制，从而在更深区域进一步提高模型分辨率和成像精度。在相同数量的数据点下，所提出的方法将选定的数据点重新分配到更深区域，从而实现比标准阵列更好的地下分辨率和更深的调查深度。通过将其应用于喀斯特洞穴探索，进一步验证了所提出方法的有效性和实际适用性。
本文的其余部分组织如下。第2节详细介绍了所提出的方法。第3节基于数值和实地实验进行了比较分析。第4节展示了其在喀斯特洞穴探索中的应用。第5节讨论了参数敏感性和深度成像性能。最后，第6节总结了结论。

2. 方法论
2.1 多步骤优化CR方法
Menke[23]引入了模型分辨率矩阵R，其主要对角元素提供了模型单元分辨率的估计。模型分辨率矩阵R表示为：
(1)
其中J是雅可比矩阵，λ是阻尼因子，C是粗糙度过滤约束矩阵。J和C定义在相同的模型离散化上，水平网格间距等于一个电极间距。
使用Sherman-Morrison Rank-1更新方法[24]来计算当新测试配置添加到基础数据集时模型分辨率矩阵ΔRb的变化：
(2)
其中g是新阵列的模型灵敏度向量，z = Bbg，y = Abz，μ = g.z。
排名函数FCR用于评估和排名由于附加配置而带来的模型分辨率改进：
(3)
其中m是模型单元的数量，Rc是全面数据集的模型分辨率矩阵。
在计算了所有候选配置的排名函数后，采用Wilkinson等人[17]描述的选择程序。在每次迭代中，将排名最高的配置（由灵敏度向量g1表示）添加到基础数据集中。随后，只有当第二个配置与第一个配置有适当的正交度时，才添加下一个排名最高的配置（g2）：
(4)
其中θ是一个指定的阈值（0.97）。计算出的点积必须低于此值，以确保配置之间的足够独立性，同时避免选定的配置之间出现过度冗余[16,17]。
对称优化的CR方法将搜索限制在对称候选集内以降低计算成本。在每次迭代中，添加的新配置数量设置为0.05 nb，其中nb表示当前基础数据集的大小，0.05表示步长。一旦基础数据集达到预定的总大小，优化过程就终止。添加新配置后，使用方程(1)重新计算扩展后的基础集的模型分辨率矩阵。
在对称优化CR方法的基础上，本文提出了一种多步骤优化的CR方法。工作流程如图1所示。首先，根据几何对称性将全面数据集分为两个子集：对称数据集和非对称数据集。优化从由偶极-偶极阵列组成的初始数据集开始，其中偶极长度a固定为一个电极间距，分离因子n的范围从1到10。
图1. 多步骤优化CR方法的流程图。所提出的算法引入了两阶段选择过程。在第一阶段，使用对称候选集进行优化。迭代过程持续进行，直到迭代次数超过预定义的阈值。然后，算法切换到非对称候选集，重复该过程，直到基础数据集达到预定的总大小。每次迭代中添加的新配置数量与对称优化CR方法保持一致。
平均相对模型分辨率Sr[25,26]定义为：
(5)
其中Sr通常用于量化优化数据集相对于全面数据集的性能。Sr值越接近1，表示优化阵列的分辨率性能越接近全面阵列的性能。

2.2 反演方法
在反演之前，检查视电阻率数据并移除异常值。反演使用L1范数反演算法进行，该算法适用于恢复具有明显对比度的电阻率分布的模型[27]。L1范数反演算法的目标函数表示为：
(6)
其中Δd表示观测值和计算值之间的视电阻率差异，Δm是模型更新向量，m表示模型参数向量，λ是正则化参数。
这个逆问题通过线性化迭代方案来解决。通过对目标函数关于Δm求导并将导数设置为零，得到相应的线性方程组：
(7)
其中 和 分别表示数据权重和模型权重对角矩阵。
采用由平均视电阻率定义的均匀半空间模型作为起始模型。L1范数反演用于迭代更新地下电阻率分布，并减少观测值和计算值之间的均方根误差（RMSE）。当RMSE降低到5%以下或最大迭代次数达到7时，反演停止。

2.3 评估指标
为了定量评估反演图像与真实电阻率模型之间的结构相似性，采用了结构相似性指数（SSIM）[28]。对于SSIM，只关注结构项s(x,y)，而忽略亮度项l(x,y)和对比度项c(x,y)。因此，当α = β = 0且γ = 1时，方程简化为：
(8)
其中x和y分别表示标准化的反演图像和标准化的真实模型。σx和σy是标准差，σxy是电阻率值的协方差。引入常数C1是为了防止分母接近零时的数值不稳定。

3. 对比研究
3.1 数值实验
构建了一个合成的双层地球模型（图2a）。上层的电阻率为500 Ω·m，厚度为40 m。上层嵌入了六个低电阻率异常（50 Ω·m）。下层的电阻率为100 Ω·m，厚度为120 m，并包含三个电阻率为1000 Ω·m的电阻异常。使用85个电极进行正向建模，电极间距为5 m，总测量长度为420 m。随后，使用有限元方法（FEM）模拟合成视电阻率数据，并添加5%的高斯噪声以模拟真实的现场条件。图2. 电阻率模型和反演结果：(a) 真实模型；(b) Wenner–Schlumberger阵列；(c) 对称优化阵列；(d) 多步骤优化阵列。黑色矩形表示异常的位置。比较研究包括Wenner–Schlumberger阵列、对称优化阵列和所提出的多步骤优化阵列。为了保持比较的一致性，优化后的数组的数据点数量被限制为与Wenner–Schlumberger数组（4992个数据点）相匹配。优化结果显示，多步优化数组的平均相对模型分辨率为0.646，超过了对称优化数组（0.615）。这些结果证明了其在解析地下结构方面的优越能力。如图2所示，所有三个数组都清晰地勾勒出了双层结构，并成功识别出了浅层低电阻率异常。多步优化数组受边界效应的影响较小，并且减轻了模型边缘附近浅层低电阻率异常的几何畸变。Wenner–Schlumberger数组在解析深层高电阻率异常的几何形状和电阻率特性方面存在困难，其SSIM值仅为0.469。相比之下，两个优化后的数组显著提高了深层成像性能，有效恢复了深层异常的几何形状、位置和电阻率值，其SSIM值分别为0.630和0.668，RMSE分别为1.4%和1.3%。值得注意的是，多步优化数组产生的反演结果更接近真实电阻率值，重建的异常边界与真实模型几何形状更为一致（图2d）。总之，与对称优化数组相比，多步优化数组在成像精度上具有明显优势。这一优势主要归因于多步优化策略。与仅限于对称数据集的对称优化CR方法不同，所提出的方法通过结合选定的非对称数据集扩展了搜索空间，这些额外的非对称数据集提供了关键的互补敏感度信息。

3.2. 管道现场实验
现场实验研究了一条直径为0.65米的排水管道。管道中心位于调查线起点的水平距离5.6米处，埋藏深度为0.7米。ERT调查线垂直于管道轴线布置，包含48个电极，电极间距为0.3米，如图3所示。在采集现场数据之前，电极阵列配置已经过优化。选择了偶极-偶极阵列作为基线配置，偶极长度a=1米，分离因子n分别为2、4、6。在优化过程中，每次迭代都会向基线阵列中添加固定数量的新数组（增量为0.02nb）。当两个优化数组的总数据点数与标准偶极-偶极数组（900个数据点）相等时，迭代过程终止。为了适应现场条件，原始电极间距为1米的优化数组配置被重新缩放以匹配实际的0.3米电极间距。优化结果显示，多步优化数组的平均相对模型分辨率为0.660，超过了对称优化数组（0.624）。

图4中的反演结果一致表明，在4.8米至6.0米的水平距离范围内存在一个明显的高电阻率异常，这与已知空心管道的位置（黑色圆圈）相符。然而，这些数组在成像精度上存在显著差异。如图4a所示，偶极-偶极数组的反演结果显示管道的响应较弱，异常区域显得模糊，内部最大电阻率仅为600.3 Ω·m，RMSE为0.9%。相比之下，对称优化数组的反演结果（图4b）识别出一个封闭的高电阻率异常，内部最大电阻率为735.4 Ω·m，RMSE降低到0.5%。多步优化数组的反演结果（图4c）进一步描绘出了一个更大且更明显的异常区域，内部最大电阻率上升至767.8 Ω·m，同时保持了较低的RMSE（0.6%）。重建的异常形状与空心管道的位置和几何形状非常吻合，这与充满空气的空腔的预期地球电响应一致。这些现场结果证实，与偶极-偶极数组和对称优化数组相比，多步优化数组提供了更高的成像精度。

4. 应用于喀斯特洞穴探索
研究地点位于中国广州市白云区（图5a）。该地区属于湿润的亚热带季风气候区，降水量丰富。根据气象记录（http://data.cma.cn），年平均降雨量和最大年降雨量分别为1800毫米和2200毫米。研究区域附近的主要河流是珠江和石井河。水文过程在喀斯特环境的发展中起着关键作用[29]。丰富的地表水和活跃的地下水循环为喀斯特的发展创造了有利条件。研究区域位于广 Cong 断层带和 Guangsan 断层带的交汇处，伴有多个次级断层（图5b）。现有地质数据显示，该地区下伏地层为石炭纪 Hutian 组（C2ht），主要由石灰岩、白云质石灰岩和白云石组成。这些活跃的构造结构与高溶解性的岩石类型共同促进了喀斯特特征的发展[30,31]。ERT调查线的布局如图5c所示，剖面从西向东方向延伸。共部署了48个电极，电极间距为3.5米，总剖面长度约为164.5米。为了改善电极与地面的耦合并确保电流注入的稳定性，在电极位置注入了盐水。在本部分研究中，使用了多步优化数组进行地下勘探。

图6显示了一个从地表延伸至大约15米深度的浅层低电阻率层（<90 Ω·m），这被解释为第四纪覆盖层。下伏层显示出相对较高的电阻率，表明其由断裂或完整的石灰岩组成。在L3和L4调查线中，覆盖层与石灰岩层之间的边界相对平滑。相比之下，L1和L2调查线在该界面显示出向下的U形变化。值得注意的是，L1 ERT剖面勾勒出了两个垂直延伸的低电阻率区域（A1和A2）：一个位于15.0–23.8米深度（距离60.0–83.0米），另一个位于16.0–22.0米深度（距离92.0–120.0米）。在L2调查线中，还描绘出了另一个向下的U形低电阻率异常（A3），位于60.0至85.0米之间，并延伸至15.0至23.8米深度。这些U形低电阻率异常被解释为含水的砂填充喀斯特洞穴的存在。

5. 讨论
5.1. 参数敏感性分析
本节研究了优化性能对关键参数的敏感性，包括步长、数组类型和网格大小。所有计算都在配备2.30 GHz Intel i7 14核Intel CPU和NVIDIA GeForce RTX 3060 GPU的工作站上进行。测试使用了具有4092个选定数据点的48电极模型。对于对称优化、多步优化和综合数组类型，分别测试了五种步长，即0.02、0.03、0.05、0.07和0.09。迭代次数、计算时间和平均相对模型分辨率总结在表1中。表1显示了使用48电极模型和4092个选定数据点时不同数组类型和步长的比较结果。对于对称优化数组，所有测试步长的平均相对模型分辨率几乎保持不变，为0.702；而计算时间从13.40秒减少到5.81秒，步长从0.02增加到0.09。多步优化数组和综合数组也表现出类似的趋势。较小的步长产生了略高的平均相对模型分辨率，但需要更多的迭代次数和更长的计算时间。这一趋势与之前的优化数组研究结果一致[18]。因此，步长主要影响迭代次数和计算时间。

不同数组类型之间的比较表明，在步长为0.05时，综合数组需要最长的计算时间（33.65秒），并实现了最高的平均相对模型分辨率（0.785），因为它搜索了所有可能的对称和非对称配置。相比之下，多步优化数组需要16.54秒，而对称优化数组仅需7.18秒。与对称优化数组相比，多步优化数组将平均相对模型分辨率从0.702提高到0.768。尽管多步优化数组在后期结合了非对称配置，从而增加了计算时间，但它在对称优化数组的计算效率和综合数组的更高分辨率之间提供了实际的折中。

除了步长和数组类型外，还在保持物理调查几何形状、步长、数组类型和选定数据点数量不变的情况下，进一步研究了网格大小的影响。如表2所示，细化水平网格使计算时间从16.54秒增加到24.33秒和46.53秒，而平均相对模型分辨率分别从0.768降低到0.753和0.735。这些变化与网格细化导致的模型单元数量增加有关。具体来说，细化网格增加了雅可比矩阵J和约束矩阵C的尺寸，从而增加了每次迭代的计算成本。同时，固定数量的数据点必须限制更多的模型单元，这增加了参数相关性，从而降低了单元模型的分辨率。因此，在本研究中考虑的多步优化数组中，1.0电极间距的网格提供了最合适的离散化方案。

5.2. 深度成像性能
本节研究了优化数组在不同深度的合成地下模型上的成像效果。所有模型都由两层结构组成：上层电阻率为100 Ω·m，厚度为10米；下层电阻率为900 Ω·m，厚度为35米。一个直径为15米、电阻率为100 Ω·m的圆形异常用黑色圆圈表示。该异常位于水平距离115米处，分别埋藏在15米、20米和25米的深度（图8）。正向建模使用了48个电极，电极间距为5米，总调查长度为235米。随后，使用FEM模拟了合成视电阻率数据。为了模拟真实的测量条件，在视电阻率数据中添加了5%的高斯噪声。反演结果清楚地表明，随着深度的增加，两种数组的成像精度都会下降。这一趋势与理论上的深度敏感度降低一致[32,33]。然而，两个阵列之间存在显著差异，这突显了阵列配置的重要性。当异常位于15米深度时（图9a,d），多步骤优化阵列能够产生一个更加明显的U形低电阻率异常，其内部最小电阻率为371.2 Ω·m，显示出其在分辨地下目标方面的优越能力。相比之下，对称优化阵列的分辨率较低，尤其是在界定异常的下边界时。当异常位于25米深度时（图9c,f），多步骤优化阵列仍然能够清晰地分辨出异常，而对称优化阵列的结果中的目标与背景难以区分。此外，多步骤优化阵列的均方根误差（RMSE）为3.8%，而对称优化阵列为3.9%。因此，在更大的深度下，多步骤优化阵列能够更可靠地保持异常的几何形状和电阻率。图9显示了使用对称优化阵列（左侧面板）和多步骤优化阵列（右侧面板）对埋藏在（a,d）15米、（b,e）20米和（c,f）25米深度的异常进行反演的结果。图10展示了两种阵列的数据点分布情况。规定的总数据点数为4092个。对称优化阵列完全由4092个对称配置组成，而多步骤优化阵列包括2288个（55.9%）对称配置和1804个（44.1%）非对称配置。每个数据点的水平位置是根据参与测量的四个电极的平均水平位置计算的，而垂直位置则是根据调查特征曲线的中值估算的[34,35]。对数据分布的统计分析表明，在大约25米以下的伪深度范围内，对称优化阵列占用了总数据量的31.0%（1270个点）。相比之下，多步骤优化阵列将这一数据比例提高到了33.0%（1349个点）。这一定量证据表明，多步骤优化阵列在深层包含更多的数据点，有效捕获了更多的深层电学信息，从而对深层结构提供了更强的数据约束。因此，综合考虑检测深度和成像精度，多步骤优化阵列的表现优于对称优化阵列。图10显示了（a）对称优化阵列和（b）多步骤优化阵列的数据点分布情况。

6. 结论与未来工作
6.1 结论
本文提出了一种多步骤优化的CR方法，以解决对称优化CR方法在成像精度方面的局限性。通过对Wenner–Schlumberger阵列、对称优化阵列和多步骤优化阵列进行数值模拟，结果表明多步骤优化阵列的平均相对模型分辨率（0.646）和SSIM（0.668）最高，超过了对称优化阵列（分别为0.615和0.630）以及Wenner–Schlumberger阵列（SSIM为0.469）。管道实验证实，多步骤优化阵列在成像精度上优于偶极-偶极阵列和对称优化阵列。应用于喀斯特洞穴勘探时，该阵列能够可靠地识别基岩界面和深层喀斯特洞穴。最后，通过引入选定的非对称数据集来扩展搜索空间，多步骤优化阵列在更深的伪深度下包含了更多的非对称数据点，从而在调查深度和成像精度方面表现出更优的性能。

6.2 未来工作
未来的工作将致力于将多步骤优化CR方法扩展到3D ERT（电阻率成像）调查中，利用其效率来降低计算成本。此外，还将开发一种基于平均相对模型分辨率增量或信息增益的自适应切换机制，以自动确定对称数据集和非对称数据集之间的最佳转换点，从而在计算效率和成像精度之间取得平衡。