Question

模型参数的先验信息

Answer 1

3.4.2.1 最小长度模型先验信息

在前面第3.2节解欠定问题时，利用了模型长度最小作为先验信息。但是一般来说我们不希望模型参数接近零的最小，而是希望与一个背景模型的差的向量长度最小。此时有如下先验信息:

地球物理反演教程

其中:

为背景模型参数列向量。这样反演的目标函数可以写为

地球物理反演教程

令

，对式(3.55)极小化，并利用和第3.2节相同的推导可得反演的解:

地球物理反演教程

解出

后，利用

就可以计算出模型参数向量。

当然我们也可以建立另一种形式的目标函数:

地球物理反演教程

注意:这里μ为拉格朗日因子(是一个数字)，用来确定模型最小化和数据拟合差最小化这两个条件在目标函数中的作用的大小。

3.4.2.2 最平模型和最光滑模型先验信息

除了用最小长度模型作为先验信息外，还可以用模型的平度、光滑度为先验信息^[2]。在数学上平度用一阶导数来量度，光滑度用二阶导数来量度。

所以我们可以在目标函数中加入模型一阶导数的模最小(平度项)或二阶导数的模最小(光滑度项)的先验信息。通常用有限差分法计算导数。

模型向量的平度项可以写为如下形式^[2]:

L=(Dm)^T(Dm)=m^TD^TDm=m^TWm (3.58)

其中:W为模型参数平度权系数矩阵;D为平度矩阵，具体形式如下:

地球物理反演教程

若模型参数的个数为N，则D为(N-1)×N矩阵。式(3.59)计算的是模型参数依次的差值向量，如果认为模型参数的间距为1，等价于各个参数对距离的一阶导数向量。

模型向量的光滑度项可以写为如下形式:

L'=(D'm)^T(D'm)=m^TW'm (3.60)

其中:W'为模型参数光滑度权系数矩阵;D'为光滑度矩阵，具体形式如下:

地球物理反演教程

若模型参数的个数为N，则D'为(N-2)×N矩阵。参看图3.5可以很容易理解光滑度矩阵。

图3.5 模型参数二阶导数有限差分公式示意图

m₂处的二阶导数有限差分公式为

地球物理反演教程

如果令Δx=1，则式(3.62)就是式(3.61)所计算向量的第一个元素。因此式(3.61)的向量D'm就是从m₂点到m_N-1点的二阶导数向量。

因此最平模型目标函数可以写为

ψ=(d-Gm)^T(d-Gm)+μ·m^TWm (3.63)

当然在欠定问题时最平模型目标函数可以写为

ψ=λ^T(d-Gm)+m^TWm (3.64)

加入光滑度限制条件只要把式(3.63)和式(3.64)中W换为W'即可。

对于二维情况，N×N平度矩阵C按下式定义^[9]:

地球物理反演教程

上式中:N为二维模型网格中节点的数量，即模型参数个数，每个节点有一个模型参数值;r_ij为i节点到j节点的距离。二维平度项写为

L_c=(Cm)^T(Cm)=m^TC^TCm=m^TW_cm (3.66)

二维模型参数网格示意图如图3.6所示。

图3.6 二维模型参数网格示意图

对于图3.6模型网格，网格x、y间距分别为Δx和Δy，则Cm向量为

地球物理反演教程

从式(3.67)可知Cm向量的第一个元素为

地球物理反演教程

式(3.68)实际上是m₁的x、y一阶导数之和。所以按照本文的定义，W_c实际上是二维模型参数平度权系数矩阵。