交叉梯度函数的MATLAB实现及代码分享01

交叉梯度函数可用于反演成像中。作为一个连接不同物性参数的桥梁，交叉梯度函数可以实现不同物性参数的联合反演成像。

一、交叉梯度函数的定义

地球物理联合反演中，建立合适的目标函数将不同物性参数联合起来是非常重要的一步工作，然而并不是所有的物性参数之间都存在较为明显的显式关系。地球物理联合反演中用于耦合不同物性参数的新方法——交叉梯度法。

该方法充分利用交叉梯度函数的性质，从结构上对不同物性参数进行耦合。在联合反演算法的目标函数中引入交叉梯度项，无需不同物性间的其他耦合关系即可实现广义意义下的联合反演。理论密度数值模型验证了交叉梯度函数的相关性质，相关应用实例也表明了该方法在联合反演中的应用效果良好，具有一定的实际应用价值。

地球物理研究都是为了解决相应的地质任务，其包含了两个基本问题：一是获得可靠的反映地下地质情况的地球物理数据和资料；二是对获得的数据和资料进行正确的解释。解释的合理与否将影响到后续的一系列生产实践，因此必须受到高度重视。

地球物理资料解释中很重要的一个步骤就是通过获得的观测数据确定地下异常体的物性参数、形态、大小、产状等，这就是通常所说的反演问题。反演技术在过去几十年间经历了从单一反演到联合反演的发展。随着勘探深度的加深，勘探目标的进一步复杂，单凭一种物性参数已经很难精确刻画地下目标体，联合反演是反演技术发展的必然趋势，许多联合反演实例证明了联合反演效果的优越性。

根据目前国内外的研究现状，从参数耦合的角度进行分类，联合反演主要分为以下三大类：
（1） 利用用于联合反演的物性参数具有相同的岩石物理学基础进行联合反演。比如：Heincke B 和 Hobbs R 利用地震波速度、重力和电阻率空间传播分布的转换关系实现了地震、重力和大地电磁数据的联合反演。
（2） 利用用于联合反演的物性参数空间分布结构的关系进行联合反演。比如：Gallardo-Delgado利用重力场和磁场空间分布的相似性进行了重力数据和磁法数据的联合反演。
（3） 不同物性参数的之间同步反演。比如：陈晓等利用正则化约束实现了地震和大地电磁测深数据的自适应正则化同步联合反演。

然而上述几种方法都有其自身的缺陷：
第一种方法需要知道两种物性参数的岩石物理学联系关系式；
第二种方法需要预先知道异常体的分布情况；
第三种方法中如何解决数据的融合问题一直是制约这种方法发展的重大障碍，其发展程度相对较低，实现起来相对困难。

基于上述方法的不足，Gallardo 和 Meju 于2003年首次提出了交叉梯度函数的概念，并将其应用于直流电测深和地震数据的联合反演，取得了良好的效果。该方法在反演过程中，无需预先明确联合反演参数间的显式解析关系式，也无需对模型做各类假设，是一种比较灵活有效的方法，这种方法的提出为解决联合反演问题提供了新的思路。

三维情况下，交叉梯度函数定义为：

交叉梯度函数在二维情况下定义为：

采用中心差分的形式对交叉梯度函数进行离散：

根据上式即可编程实现交叉梯度函数的计算。

二、交叉梯度函数的性质

在数学中，我们知道：
1） 标量场中某一点上的梯度指向标量场在这一点处增长最快的方向，梯度的大小是标量场在这一点最大的变化率；
2） 两个向量叉乘的模等于两个向量的模的乘积乘以sinθ，其中θ是两个向量的夹角。如果两个向量平行，那么夹角是0°或者180°，sinθ等于零，此时两个向量的叉乘等于零。

将数学角度的理解推广到地球物理领域，可以比较容易得出交叉梯度函数用于地球物理联合反演的理论基础：
1） 当用于联合反演的两种物性参数变化方向平行（相同或者相反），或者其中某一种物性参数不发生变化时，两者的交叉梯度函数等于零；

2） 当用于联合反演的两种物性参数的变化方向不同时，两者的交叉梯度函数不等于零。

以上两条结论就是交叉梯度函数在联合反演中能够有所应用的原理基础。在联合反演的目标函数中，当满足交叉梯度函数等于0时，则表示参与运算的两种物性变化方向平行或者其中之一不变化，也就是说两种物性所反映的地质体存在着结构上的相似性，也即存在着共同的物性边界。
在联合反演算法中引入交叉梯度函数，其会从结构上构建起两种参数之间的关系，对反演起约束，使得最终反演结果向两种物性所反映的共同边界变化，有效提高反演结果的精确性。

三、模型算例

建立了几组数值模型，来充分的认识一下交叉梯度函数。分别建立两种不同物性参数（电磁波吸收系数模型、地震波速度模型），计算他们的交叉梯度函数值。

模型一：结构一致

电磁波吸收系数分布图：

地震波速度分布图：

交叉梯度函数值分布图：

当两种物性参数的结构完全一致时，交叉梯度函数等于0.

模型二：横向结构不一致

电磁波吸收系数分布图：

地震波速度分布图：

交叉梯度函数值分布图：

当两种物性参数结构不一致，其中一种物性参数横向变化时，交叉梯度函数不等于0，反映了异常体的上下边界。

模型三：纵向结构不一致

电磁波吸收系数分布图：

地震波速度分布图：

交叉梯度函数值分布图：

当两种物性参数结构不一致，其中一种物性参数纵向变化时，交叉梯度函数不等于0，反映了异常体的左右边界。

模型四：倾斜结构不一致

电磁波吸收系数分布图：

地震波速度分布图：

交叉梯度函数值分布图：

当两种物性参数结构不一致，其中一种物性倾斜变化时，交叉梯度函数不等于0，反应了异常体的边界结构。

四、MATLAB代码分享

% 计算交叉梯度函数

tao2d = zeros(nex, ney);

for i = 2:ney-1
    for j = 2:nex-1
        tao2d(i,j) = (alpha2d(i,j+1) - alpha2d(i,j-1))/2*dx(i)* ...
                     (vel2d(i+1,j) - vel2d(i-1,j))/2*dy(j) - ...
                     (alpha2d(i+1,j)- alpha2d(i-1,j))/2*dx(i)* ...
                     (vel2d(i,j+1) - vel2d(i,j-1))/2*dy(j);
    end
end