用球面三角剖分做重力界面反演:为月球地壳厚度"量体裁衣"

论文信息:Zhang, Y., Mooney, W. D., Chen, C., & Du, J. (2019). Interface inversion of gravitational data using spherical triangular tessellation: an application for the estimation of the Moon’s crustal thickness. Geophysical Journal International, 217(1), 703–713. https://doi.org/10.1093/gji/ggz026 从重力数据里反演密度界面的起伏(例如盆地的基底界面、Moho 起伏、月球地壳厚度),是重力勘探与行星地球物理中最经典的一类问题。规模一大,就必须考虑球面曲率,问题也就自然落到球坐标系下的正演与反演上来。本文以一作张壹为代表,提出用球面三角剖分(Spherical Triangular Tessellation, STT)作为 3D 密度模型的基本表示,在球坐标系下完成密度界面的正演与反演,并将其运用于月球地壳厚度的估计。下面围绕这篇 GJI 论文,梳理研究动机、方法要点与主要结果。 为什么需要"球面三角剖分"这类表示 在大区域尺度上用重力数据做界面反演,现有方法大致分两类:一类在球谐(SH)域开展工作(如 Wieczorek & Phillips, 1998 将 Parker 公式推广到球谐域并用于月球地壳厚度);另一类在空间域把球面划分成矩形区域,用**球面棱柱单元(tesseroid)**逼近界面起伏(如 Heck & Seitz, 2007;Grombein 等, 2013;Uieda & Barbosa, 2017,后者用于南美洲 Moho)。这两类做法都很有价值,但在一些场景下存在不便: 不规则的数据覆盖与高纬度区域:球谐方法在球面上是全局、等分辨率的,难以针对稀疏或非规则观测做局部加密; 变分辨率 / 多分辨率模型:tesseroid 的矩形剖分在球面上做局部分辨率调整并不自然,相邻网格大小与形状难以平滑过渡; 全球等分辨率模型:矩形剖分在极区会出现网格严重变形、单元形状不均匀的问题。 相比之下,STT——即对球面做三角剖分(可由 Delaunay 三角剖分、递归细分,或直接连接矩形格点对角线得到)——用连续三角面片逼近曲面,单元之间没有缝隙,且能自然地实现局部分辨率的加密。这正是本文选择 STT 作为模型构建基础的出发点。 ...

August 23, 2026 · 2 分钟 · 404 字 · 张壹

用电阻率“读出”上地幔黏度:位错蠕变与电导的微观点缺陷桥梁

论文信息:Li, M., Xu, Y., Liu, L., Yang, B., Zhang, Y., Zhu, Y., & Liu, S. (2024). Physical link between effective viscosity and electrical resistivity for dislocation creep in upper mantle and its application in Northwest Xinjiang, China. Geophysical Journal International, 240, 1295–1307. https://doi.org/10.1093/gji/ggae438 上地幔的有效黏度是理解岩石圈动力学与板块运动的核心参数,但在野外观测中直接"测量"它却极其困难。好在近二十年大地电磁(MT)方法给出了一条诱人的路径:既然电阻率剖面可以刻画地下的温度、含水量与部分熔融,那么它能否进一步"换算"出黏度分布?过去的研究大多建立在"电阻率与黏度正相关"的经验假设上,其物理根基仍不够扎实。本文(第一作者为李蔓,浙江大学地球科学学院,我们课题组的合作成果之一)从点缺陷化学与电中性原理出发,为位错蠕变控制的上地幔建立了一条电阻率—有效黏度之间的微观物理桥梁,并成功应用于新疆西北部的塔里木—天山—准噶尔地区。下面梳理这篇 GJI 论文的核心思想与结果。 为什么要把黏度与电阻率联系起来 有效黏度的经验公式(Karato & Wu, 1993;Xu et al., 2018)形如 η = σ / ε̇ = A · dᵐ · Cwʳ · fO₂^q · σⁿ⁻¹ · exp[(E + P·V)/RT], ...

August 23, 2026 · 2 分钟 · 380 字 · 张壹

用神经算子为频域电磁数据建模提速:一次正演快上百倍的代理模型

论文信息: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 展平后的输出相乘,使网络能直接预测任意测点/频率上的电场磁场(或视电阻率与相位)。 ...

August 23, 2026 · 2 分钟 · 349 字 · 张壹

用等效源实现位场异常的局部分离:月球 Von Kármán 撞击坑下地幔隆起的 3D 结构

论文信息: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) 指出的,以往绝大多数分离方法针对的是埋深明显不同的源体所产生、波长差异较大的异常;而当我们想要从观测中分离出波长相近的局部异常时(例如几处相邻源体的异常叠加),基于频率差异的方法便难以奏效,这正构成了本文要解决的"局部分离"问题。 ...

August 23, 2026 · 3 分钟 · 469 字 · 张壹

用自动微分与隐式微分算地球物理反演的灵敏度与梯度

论文信息: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 无法显式地对它的解求导。常见的两个补丁都不理想: ...

August 23, 2026 · 2 分钟 · 326 字 · 张壹

用自适应梯度下降求解电磁反演问题:向深度学习借一种更省时的优化器

论文信息: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 为例,模型更新形如 ...

August 23, 2026 · 2 分钟 · 338 字 · 张壹

让大地电磁反演"又快又准":物理引导自编码器与尺度定律外推

论文信息: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 正演算子。这样一来: ...

August 23, 2026 · 2 分钟 · 341 字 · 张壹

起伏地形下的三峡水文重力效应:从"简化蓄水"到高精度负荷建模

论文信息:朱明涛,张壹,马险,王林松(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。 ...

August 23, 2026 · 1 分钟 · 212 字 · 张壹

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

论文信息: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) 光滑度叠加参考模型约束。加入参考模型约束后,异常体的边界与形态明显改善。 ...

August 23, 2026 · 1 分钟 · 191 字 · 张壹

dsv_table: 一文看懂 C++ 高性能表格数据处理

引言 在科学计算、数据分析和工程应用中,表格数据(如 CSV、TSV、空格分隔文件)是最常见的数据格式之一。然而,C++ 标准库并未提供开箱即用的表格处理工具,开发者往往需要自己编写解析逻辑,处理分隔符、表头、注释、空值等琐碎问题。 dsv_table 是 utilsxx 库中的一个核心组件,它提供了一套完整的表格数据读写与操作方案,支持多种分隔格式、灵活的行列访问、数据过滤与排序,以及 JSON 互操作。本文将带你从零开始,一文看懂 dsv_table 的设计理念与使用技巧。 一、什么是 dsv_table? dsv_table 中的 DSV 是 Delimiter-Separated Values 的缩写,泛指所有以分隔符(逗号、空格、竖线等)组织的文本表格格式。dsv_table 不仅能处理 CSV,还能处理任意自定义分隔符的文本文件。 核心设计特点 特性 说明 统一存储 所有单元格内部以字符串存储,支持按需转换为 int、double、std::string 类型标记 每个单元格带有 String / Int / Float 类型标记,影响输出行为 输出控制 可单独控制行/列/单元格的输出开关(Enable/Disable),实现"软删除" 行列命名 支持通过名称(如 "SurfaceArea_n")或内置编号(如 R3、C5)访问 多格式支持 原生支持 .txt、.csv、.json 的读写 头信息/注释/标记 文件中的注释行、标记行、头信息行会被自动提取并保留 二、快速上手 2.1 基本加载与查看 假设我们有一个 CSV 文件 sample_data.csv: id,lon,lat,depth,den,sus,name,date,location ,32.23,65.23,,2.3,,gabbro,,"enshi, hubei" ,65.2,12.7,,2.5,,basalt,,"wuhan" ,70.23,10.6,,3.1,,volcanic,,"hangzhou, zhejiang" 加载并查看列信息: #include "dsv_table.hpp" int main() { utilsxx::dsv_table t("data/sample_data", ".csv"); t.info(utilsxx::ColInfo); // 打印列信息 return 0; } 输出示例: ...

May 10, 2026 · 5 分钟 · 884 字 · 张壹