非结构网格上的非梯度正则化：让视走时与重力反演"即插即用"

论文信息：Zhang, Y., Xu, Y., Yang, B., Mooney, W. D., & Wang, F. (2022). Geophysical inversions on unstructured meshes using non-gradient based regularization. Geophysical Journal International, 230, 1864–1877. https://doi.org/10.1093/gji/ggac155

地球物理反演十有八九是病态问题，正则化几乎是"标配"。然而当模型网格从规整的长方体变成三角形、四面体构成的非结构网格时，传统基于空间梯度（x/y/z 方向偏导）的正则化算子会遇到一连串麻烦。本文由张壹担任第一作者，与合作者在三角形、四面体非结构网格上提出了一套不依赖空间梯度（non-gradient based）的正则化新算法，涵盖光滑度约束与结构相似性约束，并用二维地震走时和三维重力、三维地震走时的合成算例做了系统验证。下面围绕论文梳理其动机、方法与主要结论。

为什么非结构网格需要"另起炉灶"的正则化

在地球物理反演中，地下介质通常被离散为许多小网格单元，反演就是优化每个单元物理属性、使其既能拟合观测又具有地质意义的过程。传统的结构化网格（如规整的矩形网格、规则网格）数学处理方便，但局限也很明显：难以刻画起伏地形，难以与以三角面片组成的线框地质模型衔接，也无法根据数据覆盖或建模需求自适应地调整单元大小。

非结构网格（即三角形、四面体网格）能天然地贴合地质与地形模型，因此在电磁、大地电磁、地震走时、波形、位场等各类反演中越来越受青睐。麻烦在于，正则化公式往往需要模型参数的空间梯度（x/y/z 方向偏导），而结构化网格上单元沿轴向整齐排列，梯度很容易算；非结构网格的相邻单元界面（公共边/公共面）朝向各不相同、通常不垂直于坐标轴，直接构造轴向梯度并不顺利。若采用逐界面差分（如 Günther 等人的做法）又无法给出轴向梯度值，也就难以应付"某些方向更光滑、另一些方向更锐利"的各向异性约束；且把所有邻居都纳入计算的算法还容易遇到所谓的拼图问题（patchwork problem）——靠近等边三角形的中心单元有可能在梯度计算中被"漏掉"，完全得不到光滑度约束（图1）。

方法：非梯度正则化怎么在两个算子上实现

本文的核心思想是绕开轴向梯度，直接用非结构网格上单元的"邻居关系"来构造正则化项。

• 光滑度约束：粗糙度不再由梯度衡量，而是在相邻单元界面（三角形网格的公共边、四面体网格的公共面）上计算两侧单元的模型参数之差，并以两单元质心距离的倒数、公共边/公共面的长度（面积）加权。由于每一条公共边在粗糙度矩阵中单独成行，即使出现棋盘状分布也能被正确度量，因而天然免疫拼图问题，边界单元也无需特殊处理。

• 方向性光滑（向任意方向收紧/放松）：通过引入一个权重向量，把公共边/公共面的单位法向量与权重向量的点积（或叉积）作为权重（由指数 β 控制变化快慢）。当需要模型沿某个方向更平滑时，把权重向量转动到该方向即可——只需"旋转"一个向量，非常直观。这正弥补了传统逐界面差分无法施加各向异性约束的缺陷。

• 结构相似性约束（把参考模型的"形状"搬进来）：先获取每个单元与其邻居的若干参数差，将它们映射到参数空间。对三角形网格，三个差值构成三维参数空间中的向量，可构造准跨梯度（quasi-cross-gradient）公式（类比经典的跨梯度构造，但求叉积的坐标方向随每个三角形单元自身的三条边而定）；对四面体网格，四个差值构成四维向量，由于四维空间没有叉积，改用点积公式度量两个模型的相似性。这样即使算法已引入其他地球物理/地质数据的复杂约束，也能有效度量反演模型的结构信息。

整个反演仍是标准的 Tikhonov 框架，目标函数为数据拟合项与模型目标项之和、以正则化因子平衡；只是光滑度项与参考模型项都用上述非梯度算子构造。由于全程不出现轴向空间梯度，该方法对不同坐标系、甚至更复杂的混合（非一致）网格都易于推广。

算例与结果：二维走时、三维重力与三维走时

论文用三个合成算例检验算法。

• 二维地震首波走时（三角形网格）：地下埋有两个主要异常体与一条连接它们的窄通道。仅用基本光滑度约束时，两个主体恢复得不错，但窄连接通道丢失（高慢度小结构在反演中权重小、本就难恢复）；采用方向性光滑约束后，中间连接部分大致恢复，且反演的高慢度异常呈现与权重方向一致的优势延长方向；再叠加参考模型约束（准跨梯度、点积或经典跨梯度均可），异常体（含窄连接通道）被非常完整地恢复。不同方法得到的结构几乎相同，验证了新算子能有效度量模型的结构信息；在参考模型无结构的区域，等效于仅施加光滑度约束。

• 三维重力（四面体网格）：一个具明显延长方向的密度多面体埋于起伏地形之下。仅用最小结构+光滑度约束时基本看不出延长方向（三维重力反演常见现象，深度加权也补不了走向/倾角信息）；给定一个固定权重向量后，反演密度体的优势延长方向被"拨正"到权重方向，整体结构更贴近真实；再叠加参考模型约束，异常体轮廓基本恢复、延伸得到良好界定。

• 三维地震首波走时（四面体网格）：地下埋有三个异常体。仅用光滑度约束能恢复位置与大致取向，但形态（如表面起伏）不够准；Key (2016) 公式的结果相对发散，其中扁平异常体几乎无法从背景中区分；加入参考模型约束后，异常体边界基本跟随真实模型轮廓、显著改善形状分辨率。

讨论与结论

与已有方法相比，本文算法的主要优势在于：

• 粗糙度在单元界面上度量，无拼图问题，相邻单元之差总能被正确评估，边界单元无需特殊照顾；
• 方向性光滑约束只需提供一个权重向量即可实现，直观且易于按需给不同方向赋予不同约束（可叠加多个权重方向以表达多个优势延伸方向）；
• 结构相似性用邻居关系与参数空间的点积/叉积度量，不依赖网格几何架构，理论上对不同质量的网格都较稳健。

论文也坦承了局限：点积形式的结构相似性是通过辅助变量间接度量的，联合反演中每个迭代都要重算这些辅助变量，代价较高，因此建议点积公式仅用于构造参考模型约束；而准跨梯度算子在三角形网格上可方便地推广到采用跨梯度约束的联合反演中。作者还展望：该方法有望推广到其他坐标系，以及含多种单元类型的非一致混合网格；并计划在作者主页（person.zju.edu.cn）公布数据与代码。

延伸阅读与引用

• Günther, T., Rücker, C., & Spitzer, K. (2006). Three-dimensional modelling and inversion of dc resistivity data incorporating topography—II. Inversion. GJI, 166(2), 506–517.
• Key, K. (2016). MARE2DEM: a 2-D inversion code for controlled-source electromagnetic and magnetotelluric data. GJI, 207(1), 571–588.
• Lelièvre, P. G., & Farquharson, C. G. (2013). Gradient and smoothness regularization operators for geophysical inversion on unstructured meshes. GJI, 195(1), 330–341.
• Lelièvre, P. G., & Oldenburg, D. W. (2009). A comprehensive study of including structural orientation information in geophysical inversions. GJI, 178(2), 623–637.
• Zhang, Y., Xu, Y., & Yang, B. (2021). Lévy gradient descent: augmented random search for geophysical inverse problems. Surv. Geophys., 42(4), 899–921.

本文为对该论文的中文解读，图片均引自论文原文。