根据岩石物性参数的不同，地球物理方法可分为重力、磁法、地震和电法等，每种方法都只能从单一角度评价地下介质的岩石物理特性，单凭一种地球物理方法很难准确勘探地下结构，为了得到更准确的地下信息，利用多种地球物理参数对同一地下地质体进行综合解释已经成为当今发展的趋势(杨辉，2002).但是每种方法本身都存在一定的多解性，利用单一方法反演的单一参数模型结果进行综合解释往往很难获得结构一致的地质-地球物理模型.为了解决上述这一问题，国内外学者、专家对多物性方法多物性参数的联合反演方法进行了广泛的研究(Linde, 2008; Fregoso, 2009; Moorkamp, 2011, 2013; Lochbühler, 2013).其中包括不同物性参数之间通过经验关系式进行参数耦合，比如建立电阻率和速度之间的关系函数，Tillmann和Stöcker(2000)；Jones等(2009)开展了地震波速度和电阻率联合反演研究；(Heincke and Hobbs, 2006)和(Colombo et al., 2007)利用密度、电阻率和速度之间的先验经验函数关系，实现了电磁法、重力和地震数据的二维联合反演，但是该方法具有一定的局限性，在地下复杂的地区，很难找到正确的岩石物理关系式，因此，基于岩石物性经验关系的联合反演存在的缺陷，制约了联合反演的发展.同时，近年来采用几何空间结构耦合的联合反演方法已经成为国际上的趋势(Moorkamp et al., 2007; 于鹏等，2009; 彭淼等，2013; 王俊等，2013; Bennington, 2015; 李桐林等，2016)该方法是以不同物性在地下的岩石结构和岩石边界保持完全相同或者部分相同为前提的联合反演方法，依靠不同物性模型的空间结构相似性来耦合参数，不依赖于岩石物性关系.Gallardo和Meju(2003, 2004)首先提出了交叉梯度函数，通过对不同物性梯度之间叉乘来识别结构边界，并开展了地震走时数据和直流电阻率的二维交叉梯度联合反演研究；随后，又将交叉梯度约束加入到大地电磁和地震走时数据联合反演当中(Gallardo，2007)，以上研究中交叉梯度函数只针对两种物性参数进行约束，因此被限制在两种物理量之间的分析研究中，而在地球物理勘探中，同一地区需要进行多种物性参数的测量，综合分析研究，最终给出综合地质解释，所以基于交叉梯度的结构耦合不应该局限于两种物性参数，应该发展多种物性参数之间的相互约束.Gallardo等(2007)首次对墨西哥某一地区进行了地震初至P波和S波、直流电测深和磁法实测数据的多物性交叉梯度联合反演研究；Gallardo等(2012)又一次对巴西某地区长剖面上的海洋地震反射、海洋大地电磁、重力和磁法实测数据进行了多物性交叉梯度联合反演研究.以上联合反演适当增加物性参数，均有效降低了反演结果的非唯一性并得到了一个结构更相似的反演结果，但在上述方法中，只对野外实测数据进行处理，并没有针对结构不一致模型和复杂模型进行反演结果的验证.相比国外，国内已有的成果更多地针对二种物性之间的结构耦合，没有考虑更多不同物性参数之间的耦合问题.针对上述问题，我们建立了结构不一致模型和复杂模型，对大地电磁、重力、磁法和地震初至波走时方法的交叉梯度联合反演算法进行正确性的验证(Constable, 1987; De Groot-Hedlin, 1990; Singh, 2002; Vidale, 1990).

本文首先研究了两种地球物理方法的交叉梯度联合反演算法，在此基础上，推导并实现了多种地球物理方法(MT，重力，磁法，地震初至波走时)的多交叉梯度约束的联合反演目标函数，采用高斯牛顿方法对目标函数进行了模型改变量的公式推导，实现了多物性交叉梯度联合反演算法；其次，我们设计了结构不一致模型和复杂模型，针对多物性联合反演算法的正确性和有效性进行了模拟试算，并对单独反演结果和联合反演结果进行了对比; 最后，本文将成熟的卫星资料多光谱综合分析技术应用到复杂模型的联合反演中，将多物性参数反演模型结果图通过RGB(红-绿-蓝)模式进行合成，得到融合的RGB合成图，RGB合成图可以更好的了解地下模型情况，得到更多有价值的信息.

传统的交叉梯度联合反演中(Gallardo, 2003, 2004, 2007)只针对两种物性参数进行约束，目标函数通常包含了数据拟合项和模型平滑约束项两部分.

(1)

(2)

(3)

约束条件：τ(m₁, m₂)=∇m₁(x, y, z)×∇m₂(x, y, z)=0

其中，C_d1和C_d2分别为观测数据d₁和d₂的数据协方差矩阵；C_m1和C_m2分别为模型参数m₁和m₂的模型协方差矩阵；m₀₁和m₀₂为先验模型参数，α₁和α₂为阻尼参数，f(m₁)和f(m₂)表示为正演响应，∇为梯度，τ为交叉梯度，上述目标函数中，只有一个交叉梯度约束条件，即m₁和m₂模型参数的变化方向平行或者相反时，交叉梯度值为零.

Gallardo等(2004)和李桐林等(2015)已经给出了两种物性参数之间进行交叉梯度联合反演的公式推导，本文不再陈述.上述研究中交叉梯度函数只针对两种物性参数进行约束，因此被限制在两种物理量之间的分析研究中，而在地球物理勘探中，同一地区需要进行多种物性参数的测量，综合分析研究，最终给出综合地质解释，所以基于交叉梯度的结构耦合不应该局限于两种物性参数，应该发展多种物性参数之间的相互约束，因此本文进行了多物性(电阻率、密度、磁化率、速度)多交叉梯度约束的联合反演算法公式推导，目标函数如式(4)：

(4)

约束条件：τ(m)=0.

(5)

其中m_j为第j个物理模型的参数向量，m_i为第i个物理模型的参数向量，多物性联合反演中，不仅包含一种交叉梯度约束条件，而且需要多种交叉梯度函数的同时约束，这种思想可以很大程度的降低反演结果的非唯一性并得到一个结构更相似的反演结果.

本文采用高斯牛顿算法对交叉梯度约束下的联合反演目标函数进行数值求解，首先需要将非线性的目标函数式(4)和约束条件式(5)通过泰勒级数展开，转换成线性等式，然后使用拉格朗日算子法(Tarantola, 1987)将约束条件加入到目标函数式(4)中，得

(6)

(7)

(8)

其中，A为正演响应f(m)雅克比矩阵；B为交叉梯度函数τ(m)的导数，Λ为拉格朗日算子.

在模型空间下，本文直接对加入交叉梯度约束的目标函数式(6)求极值∂Ψ/∂m=0，得到模型更新量：

(9)

(10)

(11)

将式(9)式代入到式(4)的约束条件中可以得到拉格朗日算子Λ的表达式：

(12)

最终得到模型空间下的反演结果

(13)

本文设计了结构不一致的模型，即电阻率、密度、磁化率、速度模型异常体的边界并不完全一致，如图 1所示，单异常体电阻率模型，异常体电阻率值为10⁴Ωm，异常体大小为600 m×900 m，均匀半空间电阻率值为10^2.5Ωm；双异常体密度、磁化率模型，左侧异常体剩余密度为-0.50 g·cm^-3，磁化率为-0.1SI，异常体大小600 m×900 m；右侧异常体剩余密度为0.50 g·cm^-3，磁化率为0.1SI，异常体大小为600 m×900 m，均匀半空间剩余密度为0 g·cm^-3，磁化率为0.0001SI；单异常体速度模型，异常体速度为5000 m·s^-1，异常体大小为600 m×900 m，均匀半空间速度为4000 m·s^-1；全部异常体顶面埋深均为500 m.大地电磁数据包含了TM、TE两种极化模式下的视相位和视电阻率，频点个数为10个，频率范围在1~1000 Hz之间，9个观测点等距离分布在测线0~6 km之间；重力、磁法观测点各有30个；地震震源位于地下50 m处，共9个，接收点位于两口井中，水平位置分别为1.5 km和4.5 km，每口井内有10个接收点，间隔为0.2 km.重力、磁法、地震走时地下二维等间隔网格剖分都为140×60，MT地下二维网格剖分为174×74，需要在重磁震等间距网格剖分的基础上向外延展，延展网格以不等间隔剖分，重磁震网格剖分区域即为联合反演区域，反演网格剖分为70×30.

图 1

Fig. 1

图 1 结构不一致的理论模型 Fig. 1 Theoretical model with inconsistent structures

反演方法初始模型均采用均匀半空间模型，电阻率值为10^2.5Ωm，剩余密度为0 g·cm^-3，磁化率为0.0001SI，速度为4000 m·s^-1.结构不一致模型的单独反演结果如图 2A所示，(a1—a4)分别为电阻率、密度、磁化率、速度单独反演结果，单独反演迭代拟合差分别为0.49, 0.37, 0.56, 0.34，都收敛到均方根拟合差阈值1以下，如图 3a.在电阻率模型单独反演结果中，异常体大小明显大于真实模型，存在异常向四周发散现象；在密度和磁化率的单独反演结果中，异常体上下界面分辨能力较差，并且异常体的中心相对于真实模型偏上，没有较好的恢复异常体的空间几何形态结构，高密度和高磁异常体下部出现大面积高密度和高磁化率现象，并且围岩出现假异常，这些现象说明了重磁单独反演方法在纵向分辨能力方面较差；在地震单独反演结果中，异常体的空间几何形态有一定效果的恢复，但是恢复的结构形态相对于真实结构形态存在一定的发散现象，并且高速体左侧出现低速异常体的假异常现象，主要原因是地震射线分布不均，低速体区域射线分布相对较少.

图 2

Fig. 2

图 2 在结构不一致模型中单独反演(A)和联合反演(B)的电阻率(a1, b1)、密度(a2, b2)、磁化率(a3, b3)和速度模型(a4, b4)模型 Fig. 2 Resistivity (a1, b1), density (a2, b2), magnetization (a3, b3) and velocity (a4, b4) models obtained from separate inversion (A) and joint inversion (B) in structural inconsistency model

图 3

Fig. 3

图 3 在结构不一致模型中单独反演(a)和联合反演(b)的拟合差迭代曲线 Fig. 3 The iterative curve of the RMS misfit of the separate inversion (a) and the joint inversion (b) in structural inconsistency model

联合反演迭代拟合差曲线如图 3b所示，经过7次迭代，电阻率、密度、磁化率和速度模型最终迭代拟合差分别为0.47, 0.36, 0.23, 0.32.联合反演(图 2B)相比于单独反演，异常体的空间几何形态恢复的更好，重、磁联合反演结果中异常体下部高密度和高磁化率发散现象得到明显改善，异常体上下界面清晰可见，这是由于交叉梯度结构相似性的约束，重磁模型结构受到了地震和大地电磁模型结构的制约，同时电阻率和地震速度模型也由于重磁方法的结构相似性的影响，改善了水平分辨能力.通过不一致模型试算，可以得出，交叉梯度联合反演在地下不同物性参数结构不完全一致时，仍然可以准确的反演出不同物性参数的地下异常结构.

相比于简单的模型而言，本文又设计了一个类似于真实地下结构的复杂模型.复杂模型中包含了一些典型的地质单元体和断层结构，如图 4所示，在均匀半空间D中镶嵌了三个异常体，异常体A出露地表，下部为尖锐的阶梯状断裂结构.异常体B为孤立的地壳岩体；在模型深部区域，高阻异常体A和异常体B下面埋藏着低阻异常体C，异常体C左边出现了阶梯状的断裂构造.根据上述设计的模型我们可以用来检验不同反演方法是否具有划分横纵界面和断层、识别孤立的地层岩体和深部基底结构的能力.各方法的网格剖分和观测信息与简单不一致模型相同.

图 4

Fig. 4

图 4 复杂理论模型 Fig. 4 Theoretical model with complex structures

首先，我们对MT、重力、磁法、地震走时数据分别进行没有交叉梯度约束的传统单独反演；然后，对相同的数据进行交叉梯度约束的联合反演，所有的试验方法都是为了寻找达到期望拟合差的最优模型解.反演方法初始模型均采用均匀半空间模型，电阻率值为100 Ωm，剩余密度为0 g·cm^-3，磁化率为0.0001SI，速度为4000 m·s^-1.图 5为理论模型、单独反演结果模型和联合反演结果模型的正演响应，图 5(a1, a2, a3)分别为MT理论模型，MT单独反演和联合反演TE模式下的视电阻率拟断面图；图 5(b1, b2, b3)分别为MT理论模型，MT单独反演和联合反演TM模式下的视电阻率拟断面图；图 5(a4, a5, a6)分别为MT理论模型，MT单独反演和联合反演TE模式下的视相位拟断面图；图 5(b4, b5, b6)分别为MT理论模型，MT单独反演和联合反演TM模式下的视相位拟断面图.可以看出以上结果基本吻合一致.图 6a为重力方法的理论模型、单独反演和联合反演的正演响应曲线；图 6b为磁法的理论模型、单独反演和联合反演的正演响应曲线；图 6c为地震走时的理论模型、单独反演和联合反演的正演走时响应.

图 5

Fig. 5

图 5 MT正演响应视电阻率和视相位拟断面图 Fig. 5 Pseudo-sections illustrating forward modeling responses of apparent resistivity and apparent phase obtained by MT

图 6

Fig. 6

图 6 重力、磁法和地震走时正演响应 Fig. 6 Forward modeling responses of gravity, magnetic and seismic traveltime methods

图 7为复杂模型的单独反演结果(A)和联合反演结果(B)，反演结果由上到下分别为电阻率、密度、磁化率、速度模型，单独反演经历了6次迭代最终达到收敛，MT、重力、磁法和地震走时方法的拟合差分别为0.63，0.38，0.54，0.53，单独反演迭代曲线如图 8a所示；联合反演同样经历了6次迭代达到收敛，最终的拟合差分别为0.46，0.38，0.51，0.37，联合反演的迭代曲线如图 8b所示.在单独反演中，四种地球物理方法都有各自的优势，大地电磁方法相比于重磁法具有较高的纵向分辨能力，单独反演图 7A中，电阻率模型基本上恢复了真实模型结构，特别在纵向方面，可以不受高阻层的屏蔽，灵敏的反映出基底深部的低阻层的位置，但是在横向上，MT反演结果没有识别深部异常体C左侧阶梯状的断裂边界；重磁方法横向分辨率要高于MT和地震方法，在重磁单独反演结果中，可以分辨出水平方向的异常界面，但是纵向方面分辨能力很差，容易将深部区域的高密度、高磁化率异常体恢复在浅地表区域得到假的地下反演模型，从而对地质解释产生错误；地震方法相比于非震方法具有更好的分辨能力，但是地震射线不可能覆盖地下全部的勘探区域，图 9为地震射线分布图，在地下2~3 km处，射线分布相比于浅部区域较少，通过得到的反演结果可知，不能识别出异常体C左侧的阶梯结构边界.以上所有方法都存在各自的优缺点，如果对差异很大的单独反演结果进行地质解释，研究者很难准确的划分地质结构和判断地下地层单元.

图 7

Fig. 7

图 7 在复杂模型中单独反演(A)和联合反演(B)的电阻率(a1, b1)、密度(a2, b2)、磁化率(a3, b3)和速度模型(a4, b4)模型 Fig. 7 Resistivity (a1, b1), density (a2, b2), magnetization (a3, b3) and velocity (a4, b4) models obtained from separate inversion (A) and joint inversion (B) in complex model

图 8

Fig. 8

图 8 在复杂模型中单独反演(a)和联合反演(b)的拟合差迭代曲线 Fig. 8 The iterative curve of the RMS misfit of the separate inversion (a) and the joint inversion (b) in complex model

图 9

Fig. 9

图 9 地震速度射线分布图 Fig. 9 Distribution of rays in seismic velocity

加入交叉梯度约束的联合反演结果图 7B和单独反演结果之间具有很大的结构差异，由于结构相似性的约束，联合反演得到了结构一致性更高的反演模型结果.其中，MT反演结果在异常体C与围岩的界面处划分清晰明显，提高了深部区域横向分辨率；同时，在射线分布少的深部区域，地震反演结果也明显改善，准确的识别了深部异常体C与围岩的界面，得到了与真实模型基本吻合的地震模型结果；在重磁反演中，交叉梯度约束后的重磁反演很大程度改善了纵向分辨率不高和深部异常体上移的现象，同时可以准确的识别出孤立岩体的位置和大小.综上所述，相比于单独反演，联合反演得到的结果更接近于真实模型，无论在模型几何形态上，还是在异常体物性参数值上，都得到了有效的改善和提高.

在联合反演算法中，内部循环的目的是为了寻找结构相似性的模型，图 10展示了在每次内部循环中，交叉梯度RMS值均随着迭代次数的增加而收敛.为了进一步分析对比联合反演和单独反演结果的准确性，本文对单独反演和联合反演的最终反演模型进行了交叉梯度值的计算，如表 1所示，联合反演的交叉梯度值明显小于单独反演1~2个数量级，说明了联合反演结果的相似性程度远高于单独反演，得到的一致性模型可以更准确、更简单的进行地质解释.

图 10

Fig. 10

图 10 联合反演交叉梯度RMS值的迭代曲线 Fig. 10 Iterative curves of cross-gradient RMS value of joint inversion

表 1 (Table 1)

表 1 单独反演和联合反演的交叉梯度RMS值 Table 1 Cross-gradients RMS values for separate and joint inversion ∇mMT ∇mGrv ∇mMag ∇mSeis ∇mMT 0 (8.38×10-9) (1.79×10-9) (2.10×10-9) ∇mGrv 2.33×10-7 0 (1.51×10-9) (1.81×10-9) ∇mMag 3.36×10-8 1.15×10-8 0 (3.51×10-10) ∇mSeis 1.15×10-8 1.97×10-8 2.15×10-9 0 注：括号内的数值表示联合反演中每对物性参数之间的交叉梯度RMS值.

表 1 单独反演和联合反演的交叉梯度RMS值 Table 1 Cross-gradients RMS values for separate and joint inversion

根据2.2节复杂模型的联合反演和单独反演结果，我们绘制了单独反演图 11a和联合反演图 11b中电阻率和其他物性参数(密度、磁化率、速度)两两之间的交会图，可以发现，单独反演模型的物性参数交会图分散严重，物性之间没有明显的相互关系，无法推断出两种物性参数之间的关系表达式；对于联合反演模型的交会图，不同物性参数之间具有明显的相互关系，主要表现为线性分布结构.本文选取电阻率和密度之间的交会图作为例子，并进行详细的分析，通过电阻率和密度交会图可以得到，电阻率参数和密度参数之间表现为负相关性，即随着电阻率的增大，密度值逐渐减小，物性参数值主要沿着两个方向变化，红色方向最终延伸到电阻率1000 Ωm，密度为-0.85 g·cm^-3处，蓝色方向延伸到电阻率10 Ωm，密度为1.2 g·cm^-3处，结合本文设计的理论模型，发现此物性值基本与理论模型中异常体A、B和C物性值吻合，并且可以认为两个方向交叉点聚集处的物性值即为均匀半空间的物性值.根据上述得到的现象，我们对其他物性参数之间进行了验证，发现本文模拟的实验中，电阻率与密度、磁化率、速度具有负相关性，密度、磁化率和速度具有正相关性，并且可以通过交会图反向推断出地下异常体个数、异常体和均匀半空间的物性参数值等等，具有反向验证反演算法是否准确可信的能力.

图 11

Fig. 11

图 11 不同物性参数之间的单独反演(a)和联合反演(b)交会图 Fig. 11 Crossplots of pairs of different physical properties for models from separation inversion (a) and joint inversion (b)

为了便于反演结果的综合解释分析，我们将反演得到的电阻率、密度和速度模型合并成一个(红-绿-蓝)模式的RGB合成图.RGB颜色空间模型是一种加色混色模型，以红、绿、蓝三色互相叠加实现混色的方法.将电阻率赋予红色，色调越红表示电阻率越大；将密度赋予绿色，色调越绿表示为密度越大；将速度赋予蓝色，色调越蓝表示速度越大.通过对单独反演模型图 7A和联合反演结果图 7B，重新组合编码，生成单独反演图 12a和联合反演图 12b RGB合成图，可以发现，单独反演合成图中主要的颜色区域包括六部分，和本文设计的理论模型存在很大的差异，主要是由于不同方法的单独反演模型表现出不一致性，导致生成的RGB合成图中出现假异常现象；相比于单独反演，联合反演RGB合成图中的颜色区域明显少于单独反演，基本和理论模型相吻合，验证了联合反演方法可以得到更准确的反演结果.RGB模式在地球物理联合反演中起到了重要的作用，可以更直观的解释和分析多物性参数反演的结果，在对复杂的野外实际数据处理中的应用，更具有一定实际性和有效性.

图 12

Fig. 12

图 12 单独反演(a)和联合反演(b)的RGB合成图 Fig. 12 RGB images composed from the separate inversion (a) and joint inversion (b)

本文推导和实现了多种地球物理方法(MT，重力，磁法，地震初至波走时)的多交叉梯度约束的二维联合反演算法，并设计了结构不一致模型和复杂模型，通过对比合成数据的单独反演、联合反演结果，得到如下结论：

(1) 多物性交叉梯度联合反演结果相比于单独反演结果，无论在物性数值上，还在空间几何形态上都更接近真实模型，尤其对电阻率的横向分辨能力和密度、磁化率模型的纵向分辨能力以及地震方法由于低速层和深部区域射线分布较少导致的反演结果分辨能力差等问题的处理，均得到明显改善.

(2) 对比复杂模型联合反演和单独反演物性交会图，单独反演交会图表现出严重的发散现象，物性之间没有明显的相互关系；联合反演交会图，主要表现为线性分布结构，不同物性参数之间具有明显的相互关系，通过联合反演交会图可以反向推断出地下模型一系列相关的物性参数信息，反向证明了联合反演算法相比于单独反演更加准确和有效.

(3) 相比于单独反演，联合反演生成的RGB合成图更接近真实模型，RGB合成图中包含了全部物性参数的信息，可以更直观的分析和解释反演结果，该方法应用到综合地球物理解释中，将会成为一种有效的解释工具.