固体地球物理学是基于物理学原理，利用仪器开展多物理场观测，发展和应用地球物理正反演技术，揭示地球内部物质和结构三维分布信息，通过多种观测手段与现有物理模型对地表和地球内部的复杂物理过程进行解释和预测的一门学科。随着科技的进步和对地球认知的深入，产生了海量的多元地球观测数据；同时，反映地球运动和动力过程的物理模型也在不断更新与完善。

数据同化是一套成熟的调和理性主义与经验主义的方法论，其以贝叶斯理论、最小二乘法、变分法和控制论为基石，通过对模型与观测误差进行估计，并考虑数据时空分布，在物理模型的动态运行中不断融合新的观测数据，权衡模型与观测的不确定性，以减少系统误差，从而生成具有时间、空间和物理一致性的数据集，为系统提供更接近真实的模拟与预测^[1]。数据同化最初是为数值天气预报提供准确的初始参数，通过整合所有信息或数据来定义一个最大可能的大气状态^[2-4]，现已渗入到海洋^[5-6]、农业^[7]、现代控制等多个领域。

数据同化方法在固体地球物理学中的应用相对于大气科学一直进展缓慢，但目前的研究已经显示出其巨大的应用潜力和前景。为更好地推进数据同化方法在固体地球物理学中的应用，本文系统总结数据同化方法在固体地球物理学中的应用进展，并探讨当前亟待解决的关键问题。

1 数据同化方法概述

随着仪器精度的提升，现代固体地球物理学研究能够依赖高质量、多属性的客观观测和实验测量数据。同时，随着对固体地球认识的加深，地球物理模型也日益复杂化，但数值模拟与地球观测之间仍存在一定的差异。

在固体地球物理学领域，数值模拟是一种重要的模拟和研究地球物理过程的方法。其通常基于复杂的数学方程，这些方程依赖于初始假设与驱动数据描述地球物理系统中各状态变量动力学行为，如地震波传播、岩石圈应力分布以及地幔对流。然而，数值模拟仅为实际物理过程的简化和抽象，并受到数值计算、过程理解以及地球内部结构和物理性质认识的限制，加之初值、驱动数据以及参数的不确定性，使得误差会随模型的运行不断累积。与此同时，受到时空不连续和仪器误差的限制，观测仅能从理论上得到所测量对象在观测时刻和所代表空间上的“真值”。

为更精确地预测地球物理过程和状态，并估计地球内部性质，迫切需要通过数据同化方法将观测与模型进行融合。数据同化通过一种自洽的方法，充分利用地面实测数据在精度上和模型模拟在时间上的优势，在模型运行过程中不断加入时空不规则分布的观测数据，以不断调整模型的运行轨迹，更新模型状态，持续确定和优化模型参数，可以释放数值模拟过程中积累的误差，得到更精确的分析结果。

数据同化包括4个基本要素：模拟真实过程的物理模型、驱动模型所需的基本参数、数据同化算法、直接或间接观测的状态量。其中，数据同化算法是连接观测数据与模型预测数据的桥梁，其将观测数据不断融入物理模型，以校正模型参数和模型预测轨迹。早期数据同化算法包括多项式插值、优化插值等方法。20世纪90年代后，随着对地观测技术和地球系统科学的发展，现代数据同化技术得到更深入的应用，可大致分为连续数据同化与顺序数据同化两大类。连续数据同化是指在一个同化窗口中，利用最优算法将窗口中的全部观测数据和模型得到的状态进行最优估计，迭代调整模型初始场，最终目的是将模型轨迹拟合到同化窗口周期内的所有观测上，如三维变分(3DVar)同化和四维变分(4DVar)同化。

3DVar同化将模型模拟轨迹拟合至某一同化时刻t周围一定时间范围内的所有观测值上，利用该时间范围内的所有观测数据，并借助代价函数描述分析值与真值之间的差异。其代价函数J(X)可表示为：

$ \begin{aligned} J(\boldsymbol{X}) & =\left(\boldsymbol{X}-\boldsymbol{X}^b\right)^{\mathrm{T}} \boldsymbol{E}^{-1}\left(\boldsymbol{X}-\boldsymbol{X}^b\right)+ \\ & (\boldsymbol{Y}-\boldsymbol{G} \boldsymbol{X})^{\mathrm{T}} \boldsymbol{R}^{-1}(\boldsymbol{Y}-\boldsymbol{G} \boldsymbol{X}) \end{aligned} $

(1)

式中，X为状态量，X^b为背景场，E为背景场误差协方差矩阵，Y为观测数据，G为观测算子，R为观测场误差协方差矩阵。变分思想求解数据同化问题，就是最小化预测值与观测值之间的“距离”，以得到最优状态估计。该问题可转化为求解J(X)最小的状态量，即求解J(X)导数为0的状态值。

4DVar同化除了考虑状态的三维分布，还考虑时间的一维分布，可在一定程度上弥补3DVar同化丢失部分观测信息和模型初值与模型不协调的缺点。其代价函数J(X)为：

$ \begin{gathered} J(\boldsymbol{X})=\left(\boldsymbol{X}-\boldsymbol{X}^b\right)^{\mathrm{T}} \boldsymbol{B}^{-1}\left(\boldsymbol{X}-\boldsymbol{X}^b\right)+ \\ \sum\limits_{t=0}^{t_{\text {total }}}\left(\boldsymbol{Y}_t-\boldsymbol{G}_t\left(\mathbf{N}_t\left(\mathbf{N}_{t-1}\left(\cdots\left(\mathbf{N}_1(\boldsymbol{X})\right)\right)\right)\right)\right))^{\mathrm{T}} \\ \boldsymbol{R}^{-1}\left(\boldsymbol{Y}_t-\boldsymbol{G}_t\left(\mathbf{N}_t\left(\mathbf{N}_{t-1}\left(\cdots\left(\mathbf{N}_1(X)\right)\right)\right)\right)\right). \end{gathered} $

(2)

式中，t为时间，t_total为总体时间，B为背景场误差协方差矩阵，N为状态量随时间的变化关系，其余字母意义同式(1)。可以看出，4DVar同化引入了背景场提供的演变信息N，使其在处理非线性问题时具有较高的准确性。

顺序数据同化方法(又称滤波算法)在模型积分过程中不断通过新的观测数据更新预报场，下一时刻的状态最优估计为当前时刻观测值与预测值的加权平均，以使估计值尽可能接近观测值，并保持对模型动力学约束的一致性。卡尔曼滤波(Kalman filter, KF)是最早的顺序数据同化方法，其每步同化可分为预测与更新两个部分。对于预测部分，t+1时刻的状态预测值X_t+1^f及其误差协方差矩阵P_t+1^f可表示为：

$ \boldsymbol{X}_{t+1}^f=\mathbf{N}_{t, t+1} \boldsymbol{X}_t^a $

(3)

$ \boldsymbol{P}_{t+1}^f=\mathbf{N}_{t, t+1} \boldsymbol{P}_t^a \mathbf{N}_{t, t+1}^{\mathrm{T}}+\boldsymbol{V}_t $

(4)

式中，f表示预测，a表示分析，X_t^a为t时刻状态分析值，N_{t, t+1}为t时刻到t+1时刻的线性状态变化，P_t^a与V_t分别为k时刻分析误差协方差矩阵与模型误差协方差矩阵。当t+1时刻存在观测数据时，t+1时刻的状态分析X_t+1^a及其误差协方差矩阵P_t+1^a可通过下式进行更新：

$ \boldsymbol{X}_{t+1}^a=\boldsymbol{X}_{t+1}^f+\boldsymbol{Z}_{t+1}\left(\boldsymbol{Y}_{t+1}^O-\boldsymbol{G}_{t+1} \boldsymbol{X}_{t+1}^f\right) $

(5)

$ \boldsymbol{P}_{t+1}^a=\left(\boldsymbol{I}-\boldsymbol{Z}_{t+1} \boldsymbol{G}_{t+1}\right) \boldsymbol{P}_{t+1}^f $

(6)

式中，Y_t+1^O为t+1时刻的观测数据，G_k+1为k+1时刻的观测算子，Z_t+1为t+1时刻的增益矩阵，其公式中包含k+1时刻观测误差协方差矩阵R_k+1，具体可表示为：

$ \begin{gathered} \boldsymbol{Z}_{k+1}= \\ \left(\boldsymbol{G}_{k+1} \boldsymbol{P}_{k+1}^f\right)^{\mathrm{T}}\left[\boldsymbol{G}_{k+1}\left(\boldsymbol{G}_{k+1} \boldsymbol{P}_{k+1}^f\right)^{\mathrm{T}}+\boldsymbol{R}_{k+1}\right]^{-1} \end{gathered} $

(7)

为克服KF在计算协方差矩阵时的复杂性，采用集合思想，通过蒙特卡洛法计算预测误差协方差，即集合卡尔曼滤波方法(ensemble Kalman filter, EnKF)，可极大降低同化的计算量。另外，基于贝叶斯采样估计的顺序重要采样(sequential importance sampling, SIS)滤波思想提出的粒子滤波(particle filter, PF)，可以克服KF及其衍生算法受系统非线性误差高斯分布假设的限制，能更好地对非线性系统进行建模和估计。

2 全波形反演的数据同化方法

全波形反演(full-waveform inversion, FWI)是一种利用地震波来推导地球内部结构参数的地球物理勘查技术，对于二氧化碳封存、地热资源探查和提高石油回收率至关重要。鉴于传统FWI的计算时间较长，部分学者尝试将数据同化中的卡尔曼滤波方法融入到FWI中。顺序数据同化被视为解决地球物理逆问题的一种高效策略，能够整合之前所有时间步的数据到当前时间步的参数估计，并为此提供不确定性评估，从而减少对初始模型的过度依赖。

为降低FWI的计算成本，Li等^[8]将线性行程时间层析成像中的分层矩阵驱动KF方法应用到实际的CO₂注入模型。考虑到EnKF对非线性问题的优越适配性和FWI包含多参数的特性，Jin等^[9]使用EnKF解决一维叠前地震波形反演问题。Gineste等^[10-11]也使用该同化方法来反演一维速度剖面。为进一步适应FWI的需求，Thurin等^[12]提出将集合变换卡尔曼滤波器(ensemble transform Kalman filter, ETKF)与二维频率域粘声波全波形反演方案相结合，可获得解的后验协方差矩阵的低秩近似值。Wang等^[13]将EnKF与均匀采样技术同时应用于FWI，提出基于无放回均匀采样的集合Kalman滤波全波形反演方法(GEKUS)。传统FWI方法在模型复杂的情况下易落入局部极小值陷阱。若初始模型与真实模型之间存在较大误差，可能反演出完全错误的结果。此外，反演结果的质量在很大程度上取决于初始模型的配置。相比之下，融合数据同化算法的GEKUS方法通过整体优化反演模型参数和理论波场，依靠比较模型预测和实际观测数据来不断调整和更新模型参数。这种基于观测数据的持续校正机制使算法能够更好地探索参数空间，并逐渐引导模型集合朝着全局最优解的方向演化，从而有助于避开在复杂参数空间中的局部极小值问题。同时，这种自适应策略能使模型根据新的信息不断进化，减少对初始设定的依赖。

在该方法中，参数模型样本由无放回均匀采样方法依据初始参数模型生成，然后将参数模型与状态同时进行估计。这也体现了数据同化应用于FWI的一般思路，可表示为：

$ \begin{gathered} \boldsymbol{s}^a(x, t)= \\ \boldsymbol{s}^f(x, t)+\left(\boldsymbol{P}_{s u}^f\left(x, x_{r_1}, t\right), \cdots, \boldsymbol{P}_{s u}^f\left(x, x_{r_R}, t\right)\right) \\ \left(\boldsymbol{P}(t)+\boldsymbol{P}_{r r}(t)\right)^{-1}\left(\mathrm{~d}(t)-\left(\begin{array}{c} \boldsymbol{u}^f\left(x_{r_1}, t\right) \\ \vdots \\ \boldsymbol{u}^f\left(x_{r_R}, t\right) \end{array}\right)\right) \end{gathered} $

(8)

式中，s^a(x, t)为t时刻得到的最优解，s为包含波场信息u(距离、时间等)和参数信息p(速度、密度等)的组合，a表示分析值，s^f(x, t)为s的先验估计，f表示预测值，d(t)为波场u在接收器x_rk, k=1, …, R处观测值的组合，R为接收器数目，u^f(x_rk, t)为各接收器处的预测值，P_rr(t)为d(t)的误差协方差矩阵，P_su^f为包含波场和参数的预测误差信息的矩阵，其中下标s表示波场u与参数p的组合，su表示这个组合与波场u的协方差，P(t)为各接收器间的预测波场误差协方差矩阵的组合。式(8)得到的模型参数最优解最终通过泛函分析判定迭代是否收敛。可以看出，由于数据同化算法本身的设计特点，GEKUS方法可避免依靠伴随方程对梯度进行显式计算，而是通过EnKF框架来更新模型参数，这可以在很大程度上减少计算量。

扩展卡尔曼滤波(extended Kalman filter, EKF)是卡尔曼滤波中最基本的非线性版本，其利用当前估计的一阶泰勒展开逼近非线性，Eikrem等^[14]对迭代EKF评估二维各向同性时移全波形反演(time-lapse full waveform inversion, TLFWI)不确定性的能力进行分析，但每次迭代时数据协方差矩阵的庞大计算量会限制其在大规模问题中的应用。Huang等^[15]提出一种新的EKF变体形式——使用层次矩阵驱动的扩展卡尔曼滤波(hierarchical matrix powered extended Kalman filter, HiEKF)，该方法具有有效的时移数据的时间约束以及TLFWI结果的不确定性估计，在无需昂贵的存储和大量计算的条件下，能够获得高时空分辨率的速度变化，可用于CO₂监测模型。在此基础上，Huang等^[16]又构建出多参数HiEKF，可使速度估计更加精确。

3 数据同化在灾害防治中的应用 3.1 数据同化在火山危险性评估中的应用

为有效减轻火山灾害，火山学研究者收集火山监测数据并创建高精度数值模型，借助数据同化技术对火山的潜在危险性进行实时评估。在火山危险性评估中，需要处理大量的监测数据并应对岩浆系统的强非线性特性，EnKF由于其较低的计算成本而常被优先采用。Gregg等^[17]首次将GPS和InSAR数据通过EnKF同化到岩浆系统的粘弹性有限元模型中，实现火山活动动态模型的预测和误差估计。由于EnKF可充分估计不确定模型参数的时间演变，因此该方法后续也被应用到估算岩浆压力演变、约束储层规模、基底岩浆流入^[18]以及利用各种地面变形观测数据提供有关火山爆发前岩浆腔演变的信息^[19]。Bato等^[18]在预测岩浆房超压变化实验中以图表的形式直观展示自由模型运行与引入数据同化方法的数值模型的差异，通过采用Reverso等^[20]提出的双岩浆储层模型来模拟火山系统的动力学过程，并利用EnKF将模型与来自GNSS和InSAR数据的地表位移的合成观测进行集成，分析得到岩浆库超压的演化过程。此外，将数据同化分析结果与模型在无观测值校正的情况下自由向前传播的结果进行比对，发现模型自由运行得到的超压轨迹显著偏离真实值，而在数据同化情况下，分析集合的扩展范围显著缩小。这种差异表明，与自由模式运行相比，EnKF数据同化可提供更准确的超压估计。

此外，Albright等^[21]在EnKF中依靠GPS和InSAR数据成功预测了2008年阿拉斯加奥克莫克火山喷发，并利用热力学有限元模型(finite element model, FEM)成功预报了2018年厄瓜多尔加拉帕戈斯群岛谢拉·内格拉火山喷发^[22]。近年来，随着EnKF在火山学领域的应用不断深化，已形成多种EnKF与岩浆系统的适配机制，EnKF已经成为研究火山压力源及火山喷发机制的重要工具^[23]。

3.2 数据同化在地震强地面运动预测中的应用

地震预警(earthquake early warning, EEW)系统是防震减灾的重要手段之一，其通过对地震强地面运动进行实时监测，从而快速估算出地震响应范围，在破坏性的S波到达前发布地震预警信息。目前EEW系统大都基于源的方法进行预警，即首先快速确定地震位置和震级等震源参数，然后将震中距等参数代入强地面运动预测方程(ground motion prediction equation, GMPE)来预测强地面运动强度，但该方法存在高估或低估地面运动的问题。2011-03-11日本9.0级大地震虽然提前15 s发出预警，但由于断层破裂规模大，造成远场强地面运动被低估；此外，随后发生的多次余震又使预警系统高估了某些区域的地面运动强度^[24-25]。

为克服上述地震预警系统的问题，部分学者通过数据同化将强地面运动观测和弹性波场传播模式相结合，在无需确定震源参数的情况下快速、高精度地预测强地面运动强度。Hoshiba等^[26]基于二维弹性波传播的辐射传递模式，首先提出基于最优插值数据同化技术的数值地震动预测方法，并利用2011年日本9.0级地震以及2014年新潟县中部地震的观测数据验证了该方法可提前20 s预测强地面运动。为进一步提高强地面运动的预测精度，Wang等^[27]引入三维弹性波传播的辐射传递模式，采用Hoshiba等^[26]的数据同化方法对2011年日本9.0级地震的强地面运动进行实时预测；王勇胜等^[28]在对波场进行实时估计的同时，也对波动能量的传播进行数值模拟，从而得到目标地区的烈度信息。为了减少建筑物破坏和人员伤亡的潜在风险，Furumura等^[29]将基于快速三维有限差分方法(fast 3D finite difference method)的三维地震波传播模式与日本强震地震台网(K-NET和KiK-net)的观测数据相结合，利用数据同化的最优插值技术在超级计算机上实现了超过实际地震传播速度数倍的长周期地面运动预测。Oba等^[30]进一步采用基于格林函数的三维地震波传播模式，利用同样的同化技术实现了更快的长周期地面运动预测。

总之，基于弹性波场传播模式和强地面运动观测的数据同化技术可以实现快速、高精度的强地面运动预测，有望成为新一代的地震预警技术^[31]。

3.3 数据同化在地震预测探索中的应用

准确预测地震的发震时间、震中、震级是尚未攻克的科学难题。由于数据同化具有将观测数据与先验模型信息相结合获得更优状态估计的能力，部分学者试图将其引入地震预测研究中。Aalsburg等^[32]首次将古地震资料与Virtual California地震断层系统模型(简称VC)进行数据同化，评估约束VC模拟结果，并对加利福尼亚未来地震的发震时间和震中进行预测。为克服传统摩擦特性估算方法在处理大变量模型时计算效率较低的问题，Kano等^[33]将数据同化方法应用于高维系统，通过基于伴随的数据同化手段^[34]，利用断层的合成滑移速度数据来优化三自由度断层模型的摩擦参数，并证实该方法进行地震模拟的计算效率远高于马尔科夫链蒙特卡洛(Markov chain Monte Carlo, MCMC)方法与序列重要性采样(sequential importance sampling, SIS)。Werner等^[35]在大量数值模拟实验中结合顺序蒙特卡洛方法和EnKF用于后验分布的递归估计，证明了这两种数据同化方法相较于简单的确定性KF具有更高的可行性。Kano等^[36]进一步将此数据同化方法应用于震后余滑的时空演化分析，假设所有摩擦参数在研究区都是均匀的，将板块界面上的合成滑移速度观测值同化到描述震后余滑的基础模型中，验证了该流程用于优化摩擦参数的可行性。为进一步提高预测精度，Kano等^[37]使用2003年Tokachi-oki地震后基于GNSS数据估算的滑移速度替换原有的观测数据，再次应用于该地震的震后变形模拟，发现震后时间序列的预测能力显著提高。Hori等^[38]将地壳变形数据与俯冲带板块边界地震模型结合，利用SIS数据同化方法实现对非线性系统观测误差的校正，在输出大地震的时间、位置和震级信息的同时，还提供引起地震物理状态演变的信息，并利用合成数据估算日本南海海槽东段和西段地震复发间隔。Dinther等^[39]基于俯冲带地震序列模拟所涉及的非线性动力学，利用EnKF同化来自地表附近某点的速度和应力合成噪声数据，得到断层应力和动态强度演化的概率估计。

上述研究中，断层滑动模型的摩擦参数被纳入数据同化进行估计，或被视为已知参数，这可能导致参数估计的偏差。而数据同化在估计断层状态方面的能力受控于参数偏差，为此，Banerjee等^[40]采用顺序重要性重采样粒子滤波器同时进行状态和参数估计，以分离断层摩擦和剪应力的误差贡献，从而实现当前和未来的剪应力和滑移率的正确估计。值得注意的是，摩擦特性的估算对于监测长期慢滑移事件(long-term slow slip events, LSSEs)也至关重要。Hirahara等^[41]基于EnKF成功估算了摩擦参数和几个周期后的滑移率的演变。Diab-Montero等^[42]基于EnKF在理想模型实验中同化了断层附近处获得的噪声剪应力与速度合成值，对LSSEs中受速率与状态摩擦力控制的断层点状态进行估计。

4 数据同化在地球动力学中的应用 4.1 数据同化在地磁学中的应用

地磁数据同化(geomagnetic data assimilation, GDA)的目标是将地磁观测数据与地球动力学数值模型相结合，以更准确地估计地球外核的状态并预测地磁的周期性变化。Fournier等^[43]首次探讨了数据同化在地磁学中的应用，其利用一维磁场和速度场的简化动力学方程结合变分数据同化方法进行概念性分析，并提议将数据同化引入地磁学研究^[44]。同期，Sun等^[45]利用一维非线性磁流体动力学方程结合顺序最优插值同化方法进行相关的理论验证。此外，Liu等^[46]成功采用最优插值方案将地磁观测数据同化到地球动力学模型中。

这些简化或一维的模拟实验为后续的GDA研究打下了坚实基础。最优插值法^[47]、标准KF^[48-50]、SV同化(secular variation assimilation)^[51]等同化方法开始应用于地核动力学的观测约束上。自Kuang等^[52]提出首个地磁数据同化框架以来，GDA误差统计和模型开发工作逐渐展开^{[49, 53-60]}。Kuang等^[61]进一步尝试使用数据同化预测地磁场的长期变化，促使多种地磁预测的GDA系统被开发^[62-65]。可以看出，近十几年来，数据同化在地磁学的各个领域都取得显著进步。尽管由于其较短的发展历程，GDA在模型和观测偏差校正、预测协方差矩阵的收敛性以及非线性数据同化算法等方面仍存在挑战，但这并不会阻碍GDA的快速发展，GDA将在地球内部研究中得到更广泛的应用。

4.2 数据同化在地球动力演化过程中的应用

地球动力学专注于研究地球内部的物理过程，其目的是深入理解如地幔热对流和羽流、板块构造扩张和俯冲、沉积盆地演化、造山运动、火山和地震活动等地质现象。数据同化在地球动力学领域的核心目标是利用来自地球物理、地球化学和大地测量等多个领域的观测数据，估算地球动力模型参数，并深入理解地球动力演化过程^[66]。这一目标通常通过将与地球动力学过程相关的数学模型与现有的观测数据相结合来实现。

在该领域中，有3种常见的数据同化方法，分别为逆向平流法(backward advection, BAD)、变分法和拟逆法(quasi-reversibility method, QRV)。BAD相对简单，由于其忽略了热方程中的热扩散项，使得热平流方程可以在时间上进行反向求解，特别适用于扩散不显著的平流问题。BAD在恢复早期的底辟结构研究中应用广泛^[67-71]。20世纪末以来，该方法常用于基于当前地幔密度异质性来重建过去的地幔流动^[72-74]。由于地幔柱的平流被认为是造成地球表面热点的原因，Steinberger等^[75]通过BAD将其与大规模的地幔流场相结合，并对地幔粘度结构进行约束，以探求热点在深部地幔的运动。Conrad等^[76]在模拟地幔流动历史实验中，将BAD得到的热结构演化模型作为前向地幔对流模型的起点。Moucha等^[77]通过BAD模拟地幔对流，重建了非洲过去3 000万a的动态地形演变。

变分法旨在通过最小化目标函数来拟合预测模型与观测数据。该方法最初由Bunge等^[78]和Ismail-Zadeh等^[79-80]应用于地幔动力学模型，用于估算地质历史时期的地幔温度和流动。值得注意的是，由于Ismail-Zadeh等^[80]并未将斯托克斯方程纳入直接问题与伴随问题的迭代过程，使得变分法具有较低的计算成本。变分数据同化最常用于恢复过去阶段的地幔柱模型数值^[81-82]。Ismail-Zadeh等^[83]利用恢复模型从现在因热扩散而减弱的地幔柱中恢复了过去突出的地幔柱结构，同时探讨了热扩散和温度依赖粘度对地幔柱演化的影响。Liu等^[84]建立了同时反演地幔性质和初始条件的方法，并应用于重建Farallon板块俯冲^[85-87]和南美洲北部^[88]演化过程。Horbach等^[89]通过恢复过去4 000万a的全球地震横波研究，得到现今地幔非均质性，证明了变分法在高分辨率地幔环流模型中的实用性。Worthen等^[90]利用伴随方法解决了从地表速度观测和瞬时非线性地幔流动模型推断地幔流变参数场的问题。Ratnaswamy等^[91]提出一种基于邻接法的方法，推断地幔流动模型中的板块边界强度和流变参数。Li等^[92]建立一种基于伴随的反演方法，能够同时从当前地幔温度和历史板块运动中恢复时间相关地幔对流的初始温度条件和粘度参数。综上所述，变分法在地球动力学中具有广泛应用，并且研究者会根据各自需求为模型定制和推导适当的方法^[93-94]。

QRV的核心是在后向热力学方程中加入一个附加项，该项是正则化参数与高阶温度导数的乘积^[95]。通过最小化该正则化参数，可以实现预报模型状态与观测数据之间的优化拟合。Ismail-Zadeh等^[96]首次将QRV用于地球动力学建模，并进一步用于东喀尔巴阡山脉的地幔温度模型同化^[97]和日本群岛地幔的热状态动态恢复^[98]。此外，QRV也在重构以太平洋板块、印度板块、北大西洋板块和北大西洋区域为重点的全球地幔动力学研究中发挥着重要作用^{[74, 99-101]}。

5 结语

本文对数据同化在固体地球物理学中的应用进行初步总结。通过分析可知，数据同化在固体地球物理学中具有很大的应用潜力，但仍存在一些问题：1)观测数据的误差与不确定性。观测所携带的代表性误差难以估计，许多数据难以通过观测直接获得。尽管当前研究可以访问大量数据，但数据间可能的耦合性和观测数据所固有的代表性误差使得多源数据与模型的整合变得具有挑战性。2)物理模型只能作为真实情况的理想近似。当某一物理实体未被充分理解时，所构建的模型可能会存在固有偏差。当这种偏差过大时，即便进行误差估计，也难以实现有效的数据同化。因此，对地球内部动态的深入研究，以及对系统性偏差的纠正或适当处理至关重要。3)数据同化算法的选择与优化。数据同化起源于大气科学，其诸多方法都是为解决大气科学领域问题而设计的。鉴于固体地球物理学涉及的分支众多，对不同的物理实体选择不同的同化方法往往会得到不同的效果，为得到更优的结果，需要对多种方法进行尝试，并根据具体研究领域的特性对这些方法进行调整与优化。

虽然数据同化方法在固体地球物理学中的应用与大气和海洋科学尚有差距，但固体地球物理学各子领域的学者已对数据同化方法进行了创新和调整，并取得一系列具有独特价值的成果。这些成果充分展现了数据同化在固体地球物理学中的潜在应用价值。然而，固体地球物理学领域的数据同化发展仍处于起步阶段，尚缺乏专为固体地球物理特性而定制的数据同化方法，且由于观测数据的分辨率和精度尚待提高，模型的约束力仍有限。为增强预测能力，未来需要在固体地球物理系统中引入更先进的观测技术、方法和模型。