地球物理反演是指利用接收到的信号推测地下介质的物理性质及结构分布,不同信号能够反演出不同的物理参数信息,为地球物理解释提供证据,有助于寻找矿产、油藏等资源。近地表地质结构复杂,分布着沼泽、沙漠、戈壁、风化层等松软结构,而且存在山地等起伏地形,给地下介质反演带来很大困难。因此,利用有限的数据资料反演出近地表的地质结构,有利于地球物理深部成像。

联合反演在地球物理学的发展始于由Vozoff和Jupp(1975)提出的利用直流电磁测深数据和大地电磁测深数据联合反演一维层状电阻率模型,采用的反演算法为迭代二阶马奎特阻尼最小二乘法。随后,联合不同数据的反演逐渐发展,Gómez-Trevino和Edwards(1983)联合可控源电磁数据和直流电测深数据利用奇异值分解法(SVD)反演了地层的电阻率和厚度,Roth和Holliger(1998)利用高分辨率的瑞雷波和导波联合反演的P波速度和S波速度,Bosch和McGaughey(2001)在先验条件下利用重力数据和电磁数据联合反演二维密度结构和磁导率分布,Tryggvason等(2002)利用P波走时和S波走时信息联合反演纵横波速度和地震定位,Pieta和Bala(2008)联合直流电测深数据和电磁数据应用蒙特卡罗全局优化算法实现了并行计算反演视电阻率等。

地震数据和非地震数据之间的联合反演,也得到广泛发展。周辉等(1994)利用广义线性反演方法及非线性反演的预条件最速下降法开展了一维地震-大地电磁测深资料反演方法研究,得出顺序反演效果优于地震、电磁单独反演效果,而同时反演的效果最优。王西文(1997)采用剥离法对重力、地震资料联合反演目的层密度值,进而预测油气藏。Hering等(1995)利用一维直流电测深数据和地震面波数据联合反演得到近地表的电阻率模型和面波速度模型。Misiek等(1997)利用电法数据和勒夫波数据联合反演了Borsod地区近地表的速度和电阻率结构。Meju和Gallardo(2003)验证了近地表地区电阻率和纵波速度存在结构一致的现象,提出了交叉梯度算法,实现了二维直流电测深数据和地震走时数据的联合反演(Gallardo,Meju,2003,2004),以及二维大地电磁和地震走时数据的联合反演(Gallardo,Meju,2007)。交叉梯度在联合反演中得到极大地推广(Linde et al,2006; Fregoso,Gallardo,2009; Ogunbo,Zhang,2014)。

本文利用地震走时数据和航空电磁数据,应用交叉梯度对模型的结构限制,联合反演了近地表的速度和电阻率模型。地震走时正演是三维模拟,航空电磁正演是一维模拟,通过三维Tikhonov正则化(Tikhonov,Arsenin,1977),联合反演出三维速度模型和三维电阻率模型,再对合成棋盘模型进行了测试,并应用到一组实际数据中,检验反演结果。

为了对比联合反演和单独的走时反演、电磁反演的结果,本文对棋盘模型进行测试,真实模型如图1所示。图1a～d展示了不同深度速度模型的4个切面,可以看出速度模型内部共有9个异常块体,5个高速体和4个低速体交错放置。图1e～h展示了不同深度的电阻率模型的4个切面,同样可以看出,在电阻率模型内部,高速体位置对应高电阻率异常,低速体位置对应低电阻率异常。速度模型的观测系统,炮点在地下10 m,检波器在地表,炮间距为250 m,检波器的间距为50 m,共有169个炮点,5 929个检波器。电阻率模型的观测系统,

图1 不同深度的真实速度模型切面(a～d)及电阻率模型切面(e～h)
Fig.1 Slices of real velocity models(a～d)and real resistivity models(e～h)at different depth

发射接收装置离地表50 m,在X方向和Y方向的间距均是200 m,每个发射器伴随着一个接收器,二者之间的水平距离为7.9 m。

首先,分别进行单独的地震初至波走时层析成像和伪三维航空电磁反演,将没有异常的背景速度结构和背景电阻率结构作为反演的初始模型,得到三维速度模型和三维电阻率模型(图2)。然后,将单独反演的结果作为联合反演的初始输入模型。进行联合反演之前,我们首先需要讨论在目标函数中应用的权重参数(τ_r,τ_s,ω_r,ω_s,β和γ)的选择。根据数据和模型的二阶范数的L曲线,选定速度和电阻率模型的正则化项平滑因子τ_r和τ_s分别为0.001和0.5。设定参数ω_r=1,ω_s=1,β=1,然后测试了不同的交叉梯度的权重下的联合反演。图3为不同γ值反演的地震走时数据不匹配度、频率域航空电磁数据(Frequency Airborne Electromagnetic,简称FAEM)不匹配度、归一化的交叉梯度值以及目标函数值随迭代次数的变化,交叉梯度在第三次反演迭代时引入。图3c显示随着γ的增大,交叉梯度对反演的作用越来越大,电阻率模型和速度模型在结构上越来越相似,作为补偿,走时数据的不匹配度逐渐增加。在4组测试中,我们选择较合理的γ=3时的反演结果展示在图4中。

图2 单独走时反演的速度模型切面(a～d)及单独航空电磁反演的电阻率模型切面(e～h)
Fig.2 Slices of velocity model inversed by standalone travel time tomography(a～d)and slices of resistivity model inversed by standalone frequency airborne electromagnetic inversion(e～h)

从图4a～d中反演的速度模型中可以看出,Z为250 m,300 m和400 m的切面中有明显的低速体,说明联合反演中交叉梯度项起作用,使得电阻率模型的信息传递给了速度模型。Z为450 m的切面中低速体不明显,是由于电阻率反演在深度上的分辨率低,且航空磁测的探测深度有限,深部的电阻率反演结果和真实模型有差距,导致传递给速度模型的信息有差异,因此深部的速度模型也没有明显的低速结构。

通过以上合成模型测试,可知地震走时数据和航空电磁数据联合反演的方法能够得到比单独的地震走时反演近地表速度和单独的电磁数据反演近地表电阻率更准确的结果。通过交叉梯度的作用,速度模型和电阻率模型互相传递信息,使同一地区的2个物理特征量保持结构上的一致。

图3 不同的交叉梯度权重下反演的数据不匹配度随迭代次数的变化
Fig.3 The inversion data misfit varied with iteration number under different cross-gradient weightings

图4 交叉梯度权重γ=3时联合反演的速度模型(a～d)及电阻率模型(e～h)
Fig.4 The joint inverted velocity models(a～d)and resistivity models(e～h) with cross-gradient weighting of γ=3

我们将联合地震走时数据和航空电磁数据反演三维速度和电阻率结构的方法应用到野外采集的实际数据中。加拿大阿尔伯塔地区的地表起伏严重,走时数据采集有难度,因此走时数据的观测系统比电磁数据的观测系统较稀疏。地震数据观测系统在X方向覆盖10 000 m,Y方向覆盖7 000 m,有3 542个炮点,4 621个接收器,共接收了6 355 116条地震记录,炮间距和接收器间距都为50 m。电磁数据采集仪器为Fugro Airborne RESOLVE System,观测系统共有6 933个发射接收装置,每条测线的间距为100 m,发射器和接收器水平距离为7.9 m,高程随地表变化,平均在地表以上35 m的位置,接收到数据包括二次场与一次场的比值的实部和虚部,共有5个频率(0.390 kHz,1.787 kHz,8.391 kHz,40.740 kHz,131.630 kHz)。

首先根据地震记录在共炮域拾取初至波走时,利用拾取到的走时文件进行初至波走时反演。从图5a～d反演的速度模型4个Y方向切面(4 151 m,5 851 m,7 351 m和8 551 m处)可以看出,该研究地区地表起伏大,速度分层明显,从地表到400 m处大致为一层,400～600 m大致为第二层。在浅层分布着一些离散的、小的低速体,速度横向不均匀。利用航空电磁数据进行单独的伪三维反演出电阻率模型,从图5e～h可以看出,在浅层存在明显的低电阻率层,在Y=8 551 m的切面中,X从60 000 m到63 000 m之间凸起地表处有明显的高电阻率异常区。

图5 单独走时反演的速度模型切面(a～d)及单独航空电磁反演的电阻率模型切面(e～h)
Fig.5 Slices of velocity model inversed by standalone traveltime tomography(a～d)and slices of resistivity model inversed by standalone frequency airborne electromagnetic inversion(e～h)

进一步将单独的反演结果作为联合反演的初始速度模型和初始电阻率模型输入,设置不同的权重因子,比较反演结果。固定ω_r=1,ω_s=1,β=1,测试交叉梯度权重因子γ分别为0.1、0.5、1.0和10.0作用下电磁数据不匹配度、走时数据不匹配度、归一化的交叉梯度和目标函数随迭代次数的变化,如图6所示。由图6a看出,随着交叉梯度权重的增大,电磁数据的不匹配度没有明显的变化,仍然能够降到很小的值。由图6b看出,当γ<1时,走时数据的不匹配度仍然会随着迭代次数的增加而降低,但当γ=10时,在交叉梯度引入的第三次迭代后,走时不匹配度陡增,而且随着反演迭代次数增加而增加,这是由于交叉梯度的权重过大,强烈限制了速度模型和电阻率模型在结构上逐渐一致,而电阻率模型的分辨率较低,传递给速度模型的结构信息不准,所以在走时数据上会出现不匹配度增加的现象。由图6c看出,随着γ的增加,速度模型和电阻率模型在结构上越来越相似,交叉梯度的值越来越小,这也验证了交叉梯度在联合反演中对模型结构上约束的作用。当γ<1时,图6d中目标函数的值仍然随着迭代次数的增加而降低,γ>1时,其值随着迭代次数增加而增加,与走时数据的变化相同。进一步对比γ分别为0.1和0.5的情况,在第四次迭代之前二者数值几乎重叠,从第五次迭代开始,γ=0.5的曲线整体在γ=0.1曲线的下方,取值更小。也就是说,当γ=1时,交叉梯度的作用变大,联合反演中模型的匹配度开始以牺牲数据的拟合度来实现,速度模型和电阻率模型在结构上的相似度较高。因此,本文认为对这组地震走时数据和航空电磁数据的联合反演,γ=1时反演的速度模型和电阻率模型较为可靠。

联合反演的速度模型和电阻率模型如图7所示。反演的速度模型中,图7a中的浅层部分出现了明显的低速块体,在图7d中低速体的边缘更加清晰。反演的电阻率模型中,图7e中浅层出现横向不均匀现象,与速度模型有相似结构,图7f中在薄的低电阻率结构下分布着不均匀的高电阻率结构,图7g、7h中在X为60 000～62 000 m,凸起的地表处浅层的电阻率值偏低,这是由于在速度结构中此区域为低速区,交叉梯度的作用使得速度结构信息传递给了电阻率结构,使二者在结构上趋于一致。

图6 不同的交叉梯度权重下反演的数据不匹配度随迭代次数的变化
Fig.6 The inversion data misfit varied with iteration number under different cross-gradient weightings

图7 交叉梯度权重γ=1时联合反演的速度模型切面(a～d)及电阻率模型切面(e～h)
Fig.7 Slices of the joint inverted velocity models(a～d)and resistivity models(e～h) with cross-gradient weighting of γ=1

为了进一步对比单独反演算法和联合反演算法得到的速度结构的差异,我们分别在2个速度模型的基础上对地震数据进行静校正,然后实施叠加处理,如图8所示,箭头指示了叠加剖面变化明显的位置,尤其在图8b和d中,基于联合反演得到的速度结构计算出来的静校正处理得到的叠加剖面图显示了反射层较好的连续性和较少的起伏变化。从对比的结果可以看出,整体上联合反演得到的速度结构产生的叠加剖面要优于单独的走时反演得到的速度结构产生的叠加剖面。因此,联合航空电磁数据和地震走时数据反演近地表的电阻率和速度结构可以提供更准确的速度值来计算静校正,有利于我们进行后续的地震数据处理和地下深部成像。

图8 基于单独走时反演(a)、(c)、(e)、基于联合反演(b)、(d)、(f)的速度模型计算的静校正量得到的叠加剖面
Fig.8 Portion of CMP stacking results associated with velocity models based on travel time tomography result(a)、(c)、(e)and joint inversion result(b)、(d)、(f)

本文提出了利用地震初至波走时数据和航空电磁数据联合反演近地表三维速度和电阻率结构的方法,该方法应用了三维地震走时反演,一维电磁正演,结合三维Tikhonov正则化,并利用三维交叉梯度算子对结构的约束,实现了伪三维反演,得到近地表的三维速度结构和电阻率结构。这种伪三维反演算法大大提高了计算效率,节省了存储空间。电磁反演能对高电阻率区和低电阻率区成像,初至波走时反演对低速区成像有困难,而电磁数据和走时数据可以互相补充,联合反演同一区域的速度和电阻率结构,提高速度成像质量。通过对理论模型测试和实际数据应用,验证了这种联合反演方法的可行性和结果的可靠性。

本文仔细讨论了不同的权重因子在联合反演中的作用。交叉梯度的权重过小,联合反演的模型之间传递的结构信息不足以影响联合反演结果; 交叉梯度的权重过大,数据的不匹配度增加,反演结果出现假象。因此,在联合反演中,需要根据走时数据、电磁数据、交叉梯度以及目标函数随反演迭代次数变化的曲线判断权重因子的合理取值。在实际数据的应用中,本文对比了联合反演的速度结果和单独的走时反演的速度结果,然后分别基于这2种速度结构对数据做了静校正和叠加处理。叠加剖面显示,联合反演的速度结构优于单独反演结果。因此,电磁数据和走时数据的联合反演有助于得到更接近地下真实结构的模型。

航空电磁数据具有采集方便、数据量大的特点,能够反映地下电阻率的分布。航空电磁数据和地震走时数据联合反演近地表的电阻率结构和速度结构,为地球物理解释提供了更丰富的信息,也降低了反演的非唯一性问题。非地震数据和地震数据的结合有效地扬长避短,在勘探地球物理成像领域具有广泛的应用前景。因此,在未来联合更多的数据,比如重力数据、水文监测数据、波形数据等,对研究地球物理成像有极大的推动作用。