Research

论文信息：Peng, Z., Yang, B., Xu, Y., Wang, F., Liu, L., & Zhang, Y. (2022). Rapid Surrogate Modeling of Electromagnetic Data in Frequency Domain Using Neural Operator. IEEE Transactions on Geoscience and Remote Sensing, 60, 2007912. https://doi.org/10.1109/TGRS.2022.3222507 电磁数据的正演与反演堪称一对"既生瑜何生亮"的矛盾：反演要想高效快速，必须依赖同样快速的正演；而频域电磁正演恰恰又是出了名的"重活"。本文我们课题组（浙江大学地球科学学院，张壹为本文作者之一）与第一作者彭中等合作，提出了一种用神经算子构建的快速正演"代理模型"，把原本需要反复求解线性方程组的频域电磁正演，变成了一次次近乎实时的神经网络推理。下面围绕这篇 IEEE TGRS 论文，梳理它的动机、方法要点与主要结果。 研究背景：为什么频域电磁正演是反演的"卡脖子"环节 在准静态假设下，频域电磁（EM）地球物理观测的控制方程是一组偏微分方程——curl-curl 方程 ∇×∇×E + iωμσE = S。给定空间变化的电导率结构 σ(x)，正演就是要在地表各测点 s、各频率 f 上算出电场与磁场。传统的有限差分法（FDM）等方法必须先对计算区域做网格离散，再对每个频率求解一个大型线性方程组。当网格加密时，计算代价急剧上升。 更关键的是，电磁反演的本质就是反复调用正演：每一次迭代都要为当时的地下模型重新计算正演响应，因此反演的总耗时几乎完全被正演主导。这就构成了一个严酷的现实——反演迭代过程中，我们其实并不需要每次都把正演算到极高精度，只要是"够用"的响应即可。这给了"代理模型"一个巨大的施展空间：能否训练一个网络，让它学习从电导率结构到观测响应的映射，从而用一次快速推理替代一次昂贵的数值求解？ 方法：神经算子与"扩展"出的任意测点/频率外推 本文采用的核心是傅里叶神经算子（Fourier Neural Operator, FNO）。它把正演视为函数空间之间的映射（算子）来学习，而不是在单一网格上拟合。FNO 由三部分组成：一个提升层（lifting）把输入电导率 σ 映射到高维表征，若干傅里叶层，以及一个投影层。傅里叶层借鉴卷积定理，把卷积核在傅里叶域中做逐点相乘（用 FFT 加速），从而以很小的代价捕获远距离的非局部特征——非常契合电导率随空间高度异质这一特点，也是其计算高效的关键。 图1　本文提出的扩展傅里叶神经算子（EFNO）架构示意（引自论文 Figure 2）在标准 FNO 之上增加一个“位置网络”（location network），把测点与频率 (s, f) 作为输入，与 FNO 展平后的输出相乘，使网络能直接预测任意测点/频率上的电场磁场（或视电阻率与相位）。 ...

论文信息：Zhang, Y., Xu, Y., Mooney, W. D., & Chen, C. (2021). Local separation of potential field anomalies using equivalent sources: application for the 3-D structure of mantle uplift beneath Von Kármán crater, the Moon. Geophysical Journal International, 227(3), 1612–1623. https://doi.org/10.1093/gji/ggab307 位场（重力与磁力）观测是研究地下密度与磁性结构的重要途径，但在实际资料中，浅部小尺度目标的异常总是叠加在深部大尺度场的背景之上，如何把它们可靠地分离开来一直是位场资料处理与解释中的经典难题。本文以一作张壹为代表的研究，提出了一种基于等效源技术 + 迭代反演的位场异常局部分离新方法，并将其应用于月球 Von Kármán 撞击坑（VKC，嫦娥四号着陆区）下方地幔隆起的 3D 结构重建。下面围绕这篇 Geophysical Journal International 论文，梳理其动机、方法要点与主要结论。 为什么要做"局部分离"：区域—剩余异常分离的局限 观测位场数据是多个互相叠加的地下源体所产生的重力或磁效应之和，各源体的异常往往难以单独提取。传统上，我们把反映深部大尺度背景的异常称为区域异常，把浅部小尺度目标所对应的异常称为剩余异常，二者的分离（regional–residual separation）是位场解释中极有价值的一步。几十年来发展出了大量方法，大体可分两类： 基于数据的方法：利用区域与剩余异常在图形特征或频率上的差异来分离，例如多项式拟合（Beltrao et al., 1991）、最小曲率法（Mickus et al., 1991）、有限元（Mallick & Sharma, 1999）、三维主成分与纹理分析（Zhang et al., 2009）、双维经验模态分解（Hou et al., 2012）、维纳滤波（Pawlowski & Hansen, 1990）、独立成分分析（Forootan & Kusche, 2012）以及各类小波分析（Fedi & Quarta, 1998；Xu et al., 2009）等。 基于模型的方法：利用位场的固有特性来分离，例如向上延拓作为标准分离滤波（Jacobsen, 1987；Zeng et al., 2007）、优化滤波参数（Pilkington & Cowan, 2006）、以及基于**等效层（equivalent layer）**概念的带通滤波等（Pawlowski, 1994；Guo et al., 2013）。 这些方法虽多，却普遍存在几个不足：其一，对被分离异常的质量缺乏可量化的评价指标，选择何种分离结果仍带有较强的主观性；其二，无论是小波分析还是球谐分析本质上都是信号处理手段，所生成异常的空间物理合理性（physical plausibility）并不总能得到保证。更重要的是，正如 Li & Oldenburg (1998) 指出的，以往绝大多数分离方法针对的是埋深明显不同的源体所产生、波长差异较大的异常；而当我们想要从观测中分离出波长相近的局部异常时（例如几处相邻源体的异常叠加），基于频率差异的方法便难以奏效，这正构成了本文要解决的"局部分离"问题。 ...

论文信息：Liu, L., Yang, B., Zhang, Y., Xu, Y., Peng, Z., & Yang, D. (2024). Calculating sensitivity or gradient for geophysical inverse problems using automatic and implicit differentiation. Computers & Geosciences, 193, 105736. https://doi.org/10.1016/j.cageo.2024.105736 在地球物理反演里，灵敏度矩阵与目标函数梯度是几乎所有确定性优化算法的"燃料"，它们的准确与否直接决定反演能否收敛、收敛得多快。然而，对每个算子手工推导、编码并调试其导数公式，往往是耗时且极易出错的一环。本文（第一作者刘连，来自南方科技大学；浙江大学团队共同完成，张壹参与其中）提出并系统检验了一套自动微分 + 隐式微分（Automatic Differentiation with Implicit Differentiation，简称 ADID）的方案，能够自动而精确地算出灵敏度与梯度，并将其接入 3D 电磁法与直流电阻率反演。下面围绕这篇 Computers & Geosciences 论文，梳理其动机、方法与数值结果。 灵敏度/梯度为什么是反演绕不开的"硬骨头" 地球物理反演是非线性的、且因观测覆盖不足与模型空间存在零空间而严重欠定。求解这类问题要么走确定性路线（如 Occam、非线性共轭梯度、L-BFGS），要么走随机路线（如贝叶斯 / 转维抽样）。无论哪种，确定性框架下都要反复计算两类量：正演响应与灵敏度/梯度——后者正是目标函数对模型参数的偏导数。 对灵敏度/梯度的计算，传统上有三条路径（McGillivray & Oldenburg, 1990）： 扰动法（perturbation）：对每个模型参数额外求解一个全新的正演问题，代价随参数个数线性暴涨。 灵敏度方程法（sensitivity equation）：计算量与扰动法相当。 伴随方程法（adjoint equation，AE）：只需解少量额外的子正演来拼出灵敏度矩阵，或多次解一个伴随正演来得到梯度向量，是过去数十年里电磁反演（MT、DC、CSEM）的首选。 伴随法与灵敏度方程法都利用了线性方程组系数矩阵可复用的特性，但它们的痛点在于：需要手工推导、编码并调试相关的导数公式。尤其是当代码"看起来能跑、结果不离谱"时，隐藏在导数实现里的错误格外难揪出来。 方法：让自动微分"接管"难啃的线性求解部分 自动微分（AD）本质上是链式法则 + 编程技术：复杂的函数被拆成一系列初等代数运算（指数、对数、三角等），它们的导数都是良定义的，再用链式法则拼起来即可得到任意复合函数的导数，而且保证精确、无需手工推导。自 Sambridge 等（2007）把 AD 引入地球物理以来，它已广泛用于地震全波形反演等优化问题。 但 AD 在地球物理中一直"卡壳"在一点：地球物理正演（如 MT/DC 的电磁场）通常要求解线性方程组 Ku = s，其中 K 是系数矩阵、u 是物性场。线性求解器本身是迭代过程，不属于初等代数函数，AD 无法显式地对它的解求导。常见的两个补丁都不理想： ...

论文信息：Liu, L., Yang, B., Zhang, Y., Xu, Y., Peng, Z., & Wang, F. (2023). Solving Electromagnetic Inverse Problem Using Adaptive Gradient Descent Algorithm. IEEE Transactions on Geoscience and Remote Sensing, 61, 5902415. https://doi.org/10.1109/TGRS.2023.3239106 大地电磁（MT）等频率域电磁（EM）反演属于典型的非线性、非唯一地球物理反问题，工程上多采用线性化迭代求解。其中非线性共轭梯度（NLCG）算法因实现简单、内存需求低而广受欢迎，但它每一步都要做费力费时的线搜索来确定步长——这正是本文想要"动刀"的地方。本文作者之一张壹参与的研究把深度学习里早已成熟的自适应梯度下降（AGD）算法引入 EM 反演，用简单的代数累积代替线搜索，在保持结果可比的同时，把计算成本砍掉三分之一以上。下面围绕这篇 IEEE TGRS 论文，梳理其动机、方法要点与主要结论。 为什么 NLCG 反演"贵"在线搜索 在 NLCG 的每次迭代里，模型梯度决定了搜索方向，而线搜索负责确定步长。线搜索每次通常要多次求解相应的正演问题；即便采用高效的伴随方程法，梯度估计本身也大致要花一次正演两倍的计算量，而线搜索至少还要再额外付出一到两次正演。粗略估算下来，线搜索约占 NLCG 总成本的三分之一。对于需要反复迭代、每个样本都要正演的大型 EM 反演，这一块开销相当可观。 一个自然的想法是：能不能不做线搜索？深度学习社区在过去十年里给出了答案——自适应梯度下降（AGD）。它的核心是把历次迭代累积的梯度与模型更新用简单的代数运算来估计当前迭代的模型参数，从而彻底躲开耗时的线搜索。从 Adagrad 到 Adadelta、RMSProp，再到广泛使用的 Adam 及其变体（Adamax、NAdam、AMSGrad），这些算法在深度学习里已被验证得非常成熟。尽管在地震全波形反演中已有研究指出 AGD 相对 NLCG / L-BFGS 在时间与分辨率上的优势，但不同研究也存在分歧，且此前尚无将 AGD 用于 EM 数据反演的系统工作。这正是本研究要填补的空白。 方法：用累积梯度替代线搜索的 AGD 反演框架 AGD 家族共享一个朴素而漂亮的思想：不再费力去"试"一个标量步长，而是把每次迭代的梯度（及其平方）累积起来，用它们自适应地调节每个模型参数的学习率，从而把步长从标量变成向量（逐参数不同）。以最简洁的 Adagrad 为例，模型更新形如 ...

论文信息：Liu, L., Yang, B., & Zhang, Y. (2024). Inverting magnetotelluric data using a physics-guided auto-encoder with scaling laws extension. Frontiers in Earth Science, 12, 1510962. https://doi.org/10.3389/feart.2024.1510962 大地电磁（MT）反演本质上是强非线性、多解的逆问题，传统上依赖奥卡姆（Occam）等确定性与贝叶斯等随机方法反复迭代，计算代价高昂。神经网络（ANN）号称能在训练后"瞬时"出结果，却长期困于两大顽疾：拟合不准（恢复的模型无法贴合实测数据的反演误差（RMS））与不可复用（换个观察系统就得从头重训）。本文作者之一是本站博主张壹，这篇 Frontiers in Earth Science 论文正是针对这两个痛点而作——用一条**物理引导自编码器（PGAE）解决"物理一致性"，再用尺度定律（scaling laws）**解决"跨系统复用"，把 ANN 型 MT 反演从"实验室玩具"推向复杂真实环境。 为什么要动 ANN 的"物理一致性"与"泛化性" MT 逆问题在地球物理中出了名的难缠：观测数据有限、噪声不可避免，加上模型固有的零空间，反演结果并不唯一（Backus & Gilbert, 1967; Parker, 1983）。传统做法要么走确定性路线（Tikhonov、Occam、NLCG、ModEM 等），要么走随机采样路线（Bayesian / trans-dimensional），无一例外需要频繁调用正演、反复迭代，代价高昂。 深度学习的兴起为 MT 反演带来了"训练一次、秒出结果"的希望。常规 ANN 反演大多走监督学习路线：先构造大量地电模型作为标签，正演算出对应的 MT 响应作为特征，再训练网络拟合数据→模型的映射。其问题也随之而来：其一，标签（地电模型）本身很难构造——要么融入测井、近地表地质等先验，要么依赖"模型应当平滑"这类人为约束；其二，即便训练好了，网络拟合的是统计关系而非法则，恢复出的模型常常无法真正匹配实测数据的差异；其三，这样的网络对观察系统极其敏感，频率范围一换，原来的网络就无法复用，只能推倒重训（Ling et al., 2023; Pan et al., 2024 等）。 正是这三点，让 ANN 反演在真实与复杂场景中举步维艰——本文要逐一破解。 方法：把正演物理"焊死"进网络，再借尺度定律"搬家" **物理引导自编码器（PGAE）**的构思非常直白（图1）：自编码器本是编码器把输入映射到一个隐参数空间、解码器再从隐空间重建输入的"无标签"结构。本文做的关键改动是——把解码器直接替换成 MT 正演算子。这样一来： ...

论文信息：朱明涛，张壹，马险，王林松（2025）。考虑起伏地形的区域水文重力效应模拟与校验：以三峡库首区为例。测绘学报，54(5)，819–830。DOI: 10.11947/j.AGCS.2025.20240449 三峡大坝的截流与周期性蓄水调度，是改变长江中上游陆地水循环最直接的人为过程之一。巨大水体质量在库区聚集、涨落，其最直观的地球物理场响应当属区域重力场的变化——它不仅反映了各圈层物质的不均匀分布与运动信息，也是监测与评价库区地质构造环境稳定性（如滑坡、地震等蓄水后效）的重要参考资料。本文以三峡库首区（秭归至巴东段）为例，针对以往蓄水负荷模型"过于粗糙"、正演方法"近似误差大"两个痛点，建立了一套考虑起伏地形的区域水文重力效应高精度模拟框架，并结合绝对重力与台站连续重力观测加以校验。作者之一**张壹（浙江大学地球科学学院）**参与了本项工作，下面梳理其建模思路与主要结论。 背景动机：被"简化"掉的地形与边坡 三峡库首区长期布设了多条流动绝对重力测线（A10 型绝对重力仪，每年夏枯、冬蓄两次观测）与台站连续重力观测（gPhone 潮汐重力仪），为库水与地壳相互作用研究提供了丰富数据。然而，解读这些观测需要精确的蓄水负荷模型，而前人工作主要有两方面不足： 一是水体负荷模型过于粗糙。多数研究基于分辨率较低的数字化地形图或高程模型，即便有些借助 Landsat 影像修正了库区不连续性，仍将蓄水负荷离散化为直立长方体或棱柱体——即假设水体纵向面积不随水位变化，这与论文图1（左上角为库岸斜边坡实拍）显示的库岸斜坡明显不符，无法刻画库区水体的真实形态。 二是正演计算依赖近似。负荷格林函数积分法多集中于大尺度陆地水变化，难以用于局部小尺度重力扰动；而等效布格平板方法难以体现山区起伏地形对局部测点重力观测的实际影响——例如山谷测点会受到周围山体上方水负荷的向上引力效应。因此，要满足高精度地表重力监测的校正需求，就必须同时考虑动态水位下的库岸边坡淹盖，以及周边起伏地形含水层的变化。 方法：三角剖分 + 多面体外部引力场正演 针对上述问题，本文构建了库区蓄水负荷模型与区域水文负荷模型两套模型，核心思路是把复杂的水体表面与地形"如实"剖分后再做解析正演： 水体边界提取：选取 11 景拍摄于三峡水库最低（约 145 m）与最高（约 175 m）水位期间的高分辨率高分一号（GF-1）卫星影像，经预处理、影像融合后，通过归一化水体指数（NDWI）与 Otsu 最优阈值分割提取最低、最高水位边界，并结合人工目视解译消除干扰误差；再结合坝前实时水位数据构建库区动态蓄水负荷模型。 几何剖分：采用 Delaunay 三角剖分方法（符合空圆性、最大化最小角特性，网格规整均匀）将复杂水体表面、库岸边坡及起伏地表含水层剖分为三角形面，河道边界按 8 m 边界点距加密剖分，网格向河道中心递增。两套模型都由数百万个三角形面连接成外表面的复杂多面体。 正演计算：采用文献提出的均匀多面体外部引力场解析解算法（水体密度 ρ=1×10³ kg/m³），引力位、引力与引力梯度张量均以解析式给出，理论上不存在计算近似误差。本文关注水体负荷直接引力的垂向分量及重力梯度（设定重力向下为正）。 在正式建模前，作者还通过一个"棱柱体 vs 边坡六面体"的模型试验，量化了库岸边坡水体对局部重力场的影响：当测点靠近模型时，边坡水体与简化棱柱体的重力效应差异在模型边界 1 km 以内可达上百微伽（论文图5 红线），从机理上论证了用复杂边坡模型替代简化棱柱体的必要性。 图1　库首区 145–175 m 蓄水产生的重力、重力梯度效应空间分布（引自论文 Figure 6）模拟结果显示蓄水引起的重力响应主要受库区水体分布控制、集中于库区沿岸，并随距库区距离增大而迅速衰减；当水位变化 30 m 时，重力变化最大可达 200 μGal，重力垂直梯度最大约 10 E、水平梯度最大约 1 E。 结果与校验：动态蓄水模型显著优于静态模型 库区蓄水重力效应（图1）。对 145–175 m 水位变化激发的重力与重力梯度空间分布进行模拟，结果表明蓄水引起的重力响应主要受库区水体分布控制并集中在库区沿岸，随距离迅速衰减；当库区水位变化 30 m 时，重力变化最大可达 200 μGal，重力垂直梯度最大约 10 E，重力水平梯度最大约 1 E。 ...

论文信息：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）。 图1　二维地震走时反演中不同正则化约束的对比（引自论文 Figure 6）(b)(c)(d) 分别由光滑度约束（本文方法）、Key (2016) 公式与 Lelièvre & Farquharson (2013) 方法得到，结构相近；(e) 加入方向性光滑约束后，异常呈现与权重方向一致的优势延长；(f)(g)(h) 再叠加参考模型约束，含窄连接通道的异常体被完整恢复。 方法：非梯度正则化怎么在两个算子上实现 本文的核心思想是绕开轴向梯度，直接用非结构网格上单元的"邻居关系"来构造正则化项。 光滑度约束：粗糙度不再由梯度衡量，而是在相邻单元界面（三角形网格的公共边、四面体网格的公共面）上计算两侧单元的模型参数之差，并以两单元质心距离的倒数、公共边/公共面的长度（面积）加权。由于每一条公共边在粗糙度矩阵中单独成行，即使出现棋盘状分布也能被正确度量，因而天然免疫拼图问题，边界单元也无需特殊处理。 方向性光滑（向任意方向收紧/放松）：通过引入一个权重向量，把公共边/公共面的单位法向量与权重向量的点积（或叉积）作为权重（由指数 β 控制变化快慢）。当需要模型沿某个方向更平滑时，把权重向量转动到该方向即可——只需"旋转"一个向量，非常直观。这正弥补了传统逐界面差分无法施加各向异性约束的缺陷。 结构相似性约束（把参考模型的"形状"搬进来）：先获取每个单元与其邻居的若干参数差，将它们映射到参数空间。对三角形网格，三个差值构成三维参数空间中的向量，可构造准跨梯度（quasi-cross-gradient）公式（类比经典的跨梯度构造，但求叉积的坐标方向随每个三角形单元自身的三条边而定）；对四面体网格，四个差值构成四维向量，由于四维空间没有叉积，改用点积公式度量两个模型的相似性。这样即使算法已引入其他地球物理/地质数据的复杂约束，也能有效度量反演模型的结构信息。 整个反演仍是标准的 Tikhonov 框架，目标函数为数据拟合项与模型目标项之和、以正则化因子平衡；只是光滑度项与参考模型项都用上述非梯度算子构造。由于全程不出现轴向空间梯度，该方法对不同坐标系、甚至更复杂的混合（非一致）网格都易于推广。 算例与结果：二维走时、三维重力与三维走时 论文用三个合成算例检验算法。 二维地震首波走时（三角形网格）：地下埋有两个主要异常体与一条连接它们的窄通道。仅用基本光滑度约束时，两个主体恢复得不错，但窄连接通道丢失（高慢度小结构在反演中权重小、本就难恢复）；采用方向性光滑约束后，中间连接部分大致恢复，且反演的高慢度异常呈现与权重方向一致的优势延长方向；再叠加参考模型约束（准跨梯度、点积或经典跨梯度均可），异常体（含窄连接通道）被非常完整地恢复。不同方法得到的结构几乎相同，验证了新算子能有效度量模型的结构信息；在参考模型无结构的区域，等效于仅施加光滑度约束。 三维重力（四面体网格）：一个具明显延长方向的密度多面体埋于起伏地形之下。仅用最小结构+光滑度约束时基本看不出延长方向（三维重力反演常见现象，深度加权也补不了走向/倾角信息）；给定一个固定权重向量后，反演密度体的优势延长方向被"拨正"到权重方向，整体结构更贴近真实；再叠加参考模型约束，异常体轮廓基本恢复、延伸得到良好界定。 三维地震首波走时（四面体网格）：地下埋有三个异常体。仅用光滑度约束能恢复位置与大致取向，但形态（如表面起伏）不够准；Key (2016) 公式的结果相对发散，其中扁平异常体几乎无法从背景中区分；加入参考模型约束后，异常体边界基本跟随真实模型轮廓、显著改善形状分辨率。 图2　三维地震首波走时反演在不同约束下的对比（引自论文 Figure 13）(a) 仅光滑度约束；(b) Key (2016) 公式；(c) 光滑度叠加参考模型约束。加入参考模型约束后，异常体的边界与形态明显改善。 ...

目的 分析与测试带动量的LGD方法的搜索效率与算法特性。 算法分析 LGD算法的核心是步长随机的梯度下降法，其中步长取值符合莱维分布。这样做带来的好处主要有两个： 因为步长的最大取值为无限大，因此理论上算法具有全局收敛性。给定足够的搜索次数，算法一定会发现全局最优解； 搜索过程顺带对解空间的形态（起伏）进行了探索，因此可利用搜索路径的统计信息给出解的不确定度。 但也有一个缺点：由于在地球物理反演中一般会采用正则化的方法改善反演问题的奇异性，因此一些解空间的凹性质并不突出，这使得LGD在大多数迭代步骤中下降量不足。为了增加LGD的下降效率，同时保留其优点。考虑到Adam这类的最优化方法具有良好的下降效率和一定的全局最优性。因此计划在LGD算法中引入类似的动量概念，提升迭代效率。 下面是带有动量的LGD算法的简单分析： // 这里我们考虑一个二维最优化问题 // 首先计算归一化的迭代方向（这使得对动量和梯度模的直接统计没有了意义） // 迭代方向乘步长得到迭代矢量（这是我们真正想记录的量，它直接决定了迭代路径） direct_mod = sqrt(g2[0]*g2[0] + g2[1]*g2[1]); g2[0] = levy_length*g2[0]/direct_mod; g2[1] = levy_length*g2[1]/direct_mod; // 计算修正系数 mt *= beta_m; vt *= beta_v; // 利用指数平均对动量与梯度模进行估计（有偏的） m[0] = (beta_m*m[0] + (1.0 - beta_m)*g2[0]); m[1] = (beta_m*m[1] + (1.0 - beta_m)*g2[1]); v[0] = (beta_v*v[0] + (1.0 - beta_v)*g2[0]*g2[0]); v[1] = (beta_v*v[1] + (1.0 - beta_v)*g2[1]*g2[1]); // 对动量与梯度模进行校正 mhat[0] = m[0]/(1.0 - mt); mhat[1] = m[1]/(1.0 - mt); vhat[0] = v[0]/(1.0 - vt); vhat[1] = v[1]/(1.0 - vt); // 利用动量与梯度模进行迭代（此处不再需要乘以步长了） x2[0] = x2[0] - mhat[0]/(sqrt(vhat[0]) + 1e-8); x2[1] = x2[1] - mhat[1]/(sqrt(vhat[1]) + 1e-8); 模型试验 对比二维非凸问题的搜索过程（见下图，其中白色为LGD_Momentum的搜索路径，黑色为LGD搜索路径，红色为Adam算法的搜索路径），在实验过程中，两种LGD算法每次迭代的步长取值都是一样的，Adam算法的步长为LGD搜索步长的平均值，LGD算法的搜索次数为1000次。对比下图可见： Adam算法的搜索范围很小，直接陷入了局部极小值； LGD算法的搜索范围最大，路径总体呈放射状，搜索效率受到了一定的影响； LGD_Momentum算法的搜索范围较小一些，全局收敛性较LGD更弱，但更集中于极值附近，搜索的效率整体更高。 迭代路径对比图 问题 怎样更直观地对比LGD与LGD_Momentum的求解效率？ LGD_Momentum在地球物理反演中是否更优？

View on Github 线性共轭梯度算法库（C++ Library of Linear Conjugate Gradient，LIBLCG） 张壹（zhangyiss@icloud.com） 浙江大学地球科学学院·地球物理研究所 简介 liblcg是一个简单的C++线性共轭梯度算法库，包含了一般形式的共轭梯度与预优共轭梯度算法。可用于求解如下形式的线性方程组： Ax = B （可用共轭梯度求解）或者 PAx = PB（可用预优共轭梯度求解） 其中，A是一个N阶的对称矩阵、x为N*1的待求解的模型向量，B为N*1需拟合的目标向量。P则为预优算法中的预优矩阵，一般是一个N阶的对角阵。共轭梯度法广泛应用于无约束的线性最优化问题，拥有优良收敛与计算效率。 安装 GCC（MacOS） mkdir build cd build cmake .. make make install Windows 请自行拷贝代码新建项目并编译…嘿嘿😁！ 使用说明 自定义Ax计算函数 通常我们在使用共轭梯度法求解线性方程组Ax=B时A的维度可能会很大，直接储存A将消耗大量的内存空间，因此一般并不直接计算并储存A而是在需要的时候计算Ax的乘积。因此用户在使用liblcg时需要定义Ax的计算函数，同时提供初始解x与共轭梯度的B项（即拟合的对象）。如果使用预优方法还需要提供预优矩阵P项。Ax计算函数的形式必须满足算法库定义的一般形式： typedef void (*lcg_axfunc_ptr)(void* instance, const lcg_float* x, lcg_float* prod_Ax, const int n_size); 函数接受4个参数，分别为： void *instance 传入的实例对象； const lcg_float *input_array Ax计算中的x数组的指针； lcg_float *output_array Ax的乘积； const int n_size 矩阵的大小。 自定义进程监控函数 用户可以以下面的模版创建函数来显示共轭梯度迭代中的参数，并可以在适当的情况下停止迭代的进程。 typedef int (*lcg_progress_ptr)(void* instance, const lcg_float* m, const lcg_float converge, const lcg_para* param, const int n_size, const int k); 函数接收6个参数，分别为： ...

Program Propose 3D model construction and forward calculation of gravity and magnetic data using the Cartesian coordinates. density or magnetic models are built using elements of rectangular blocks. The program can either build 3D models or forward calculating gravity and magnetic data from a input model file. Some typical source types are supported by the program for fast model construction. File format of the 3D model used in this program is the 2.0 .msh file of the Gmsh software. ...