近年来，国内外金属矿地震勘探技术发展较为迅速，已经成为寻找深部矿产的重要方法之一.我国在金属矿地震理论和应用方面也开展了一些探索，取得了一定进展(吕庆田等，2010，2014；廉玉广等，2011).金属矿区人文环境与地质构造复杂、目标地质体与围岩波阻抗差较小，造成了金属矿反射地震信号能量弱和呈不连续的非层状反射，或者散射特征，往往还会受到强背景噪声的干扰，导致利用常规的速度建模技术构建的金属矿地震速度模型的精度较低，偏移成像效果差.此外，金属矿目标地质体与围岩之间波阻抗差异较小、或者存在渐变等特征，常规的纵波勘探难以取得预期效果.而基于弹性波多参数的多波多分量勘探方法有望克服传统的纵波地震勘探反射信息弱的问题(李庆春等，2007)，但是常规的速度建模和成像方法的精度和分辨率对于大倾角等复杂介质都会有很大的误差(韩建光等，2013)，为此需要发展高精度地震波速度和密度等多参数模型重建和偏移成像的方法.

地震全波形反演(FWI)理论于20世纪80年代建立(Tarantola，1984; Mora，1987)并逐步得到完善和改进，近年来更是取得了迅速发展(Virieux and Operto，2009；韩淼，2014).同时，弹性波全波形反演(EFWI)的计算速度和精度也有了很大的提高(Castellanos et al.，2011; Wang et al.，2012; Xiong et al.，2013; Son et al.，2014).继Pratt等(1998)将声波近似的全波形反演方法引入到频率域后，Brossier等(2009)在频率域实现了弹性波FWI，并采用多尺度策略提高反演精度.而时间域的多尺度反演方法能灵活选取时窗分离波场，便于降低反演问题的非线性性(Boonyasiriwat et al.，2009; Fichtner et al.，2013)，更适用于处理波场复杂的实际资料.由于波场中各属性参数相互耦合，因此利用弹性波全波形反演方法重建地下介质波场参数分布更加困难.Köhn等(2012)详细讨论了时域EFWI选用不同的参数组合，如拉梅系数，地震波速度或波阻抗作为反演参数时，参数间的相互影响.Xu和McMechan(2014)综合考虑各个参数的梯度变化，在得到多参数相互独立的搜索方向和步长后，沿着多参数合成的搜索方向寻找一个最优步长，以较快的收敛速度实现任意参数组合的反演；Groos等(2014)引入黏弹性介质的品质因子同时反演浅部体波和瑞雷波；Ren和Liu(2015)提出联合变差分阶数正演和基于子波的多尺度反演策略.这些方法的提出使得弹性波全波形反演更加成熟，已经成功应用到了海洋(Prieux et al.，2013；Raknes et al.，2015)和陆地地震勘探实际数据的多参数建模中(Plessix and Pérez，2015；Borisov and Singh，2015).

地震照明技术通过建立地质模型人工模拟野外采集方式，分析地震波在介质中传播过程中能量的分布，评价观测系统参数(谢小碧等，2013)，提高AVO/AVA分析精度(Askim et al.，2010)和优化偏移成像剖面(Zhou et al.，2014).Yang等(2013)将照明能量作为加权因子引入到成像域声波全波形反演中，消除建立速度模型过程中照明不均的影响.Bian等(2015)利用照明技术估计全波形反演中观测系统参数对反演结果的作用.陈永芮等(2013)利用双向照明优化观测系统并对声波全波形反演梯度进行能量补偿.应用全波形反演处理实际数据时，可以根据先验模型进行可视性分析找到最优炮点分布和最佳接收范围，并在对应位置加密震源和检波器，使观测数据包含更多的有效信息，提高反演精度和稳定性.由于时间域全波形反演的计算复杂度与模型大小及观测系统有直接关系，随着震源点数的增加，计算量成正比增长.Wang和Rao(2009)在声波全波形反演中采用层剥离的反演方法，解决了有限偏移距的波场效应和反演误差传递等问题，逐层反演以重构深部纵波速度模型，但这种方法的计算量较大，正比于对底层的人为分层数目.本文提出一种基于可视性分析与自适应能量补偿的弹性波全波形反演方法，并应用于复杂的金属矿地震勘探多参数建模处理.首先将观测系统中单个震源和检波器对目标体的总照明强度定义为该炮检对对目标体的可视性，对弹性波全波形反演中的观测系统进行可视性分析并计算出可视性分析平面图，将可视性作为加权因子引入到全波形反演中目标函数的计算，给出基于可视性的目标函数，以增加观测波场和模拟波场相互匹配过程中与目的层相关性较高的炮检对的残差权重.其次，对于金属矿复杂地区存在的照明 阴影区而导致弹性波场能量不均匀分布、成像效果差的问题，根据给定的观测系统计算波场对地下介质的双向照明强度，构建自适应的加权梯度函数，使得地震波场能量近似均衡分布，从而改善全波形反演效果.

全波形反演中的目标函数E可以表示为

目标函数E对模型参数m求导可得到梯度方向为

其中， u _cal和u _obs分别是模拟波场和观测波场，δ d 是两者间的残差.若雅克比矩阵∂ u /∂ m 已知，则对模型空间中所有扰动进行积分可计算数据空间波场的总变化值. 同理，对数据空间的波场扰动求积分可以得到模型空间的变化量为(Tarantola，2005)

其中，*代表相应的共轭矩阵.由于这种反问题的解并非唯一，即δ u ≠δ u ′，δ m ≠δ m ′，但由于[∂ u(m)/∂ m ]^*=∂ u(m)/∂ m，因此若将数据空间的扰动看作残差，则从数据空间到模型空间的映射就等同于目标函数的梯度：

因此，根据一阶扰动形式下的弹性波方程和格林函数，可以得到纵波速度v_p、横波速度v_s和密度ρ的梯度表达式为(Mora，1987)

其中，u_x和u_z分别是在震源点激发的弹性波场的水平分量和垂直分量，φ_x和φ_z分别是在检波点位置激发的反传残差波场的水平分量和垂直分量.选用合适的优化方法，沿着目标函数的负梯度方向，利用迭代公式(9)对初始模型 m ₀迭代更新，即可得到目标函数的全局最优解，公式(9)为

其中，α_n为第n 次迭代的步长， P 为预处理算子.若选用牛顿法优化参数模型，则 P 为海森矩阵的逆，若选用高斯牛顿法优化参数模型，则 P 为对角海森矩阵的逆.

金属矿地震勘探由于地表条件的复杂和地下介质的不均匀分布等因素的干扰，只有部分位置的炮点和检波点能激发并接收到与探测地质体有关的地震波信息，通过波动方程照明分析，可以得到观测系统中单个震源和检波点对目的层的归一化总照明强度即观测系统的可视性，它能直观地描述观测系统中任意炮捡对的观测数据对于成像或反演的有效作用程度.计算基于可视性的目标函数来优化全波形反演中对应炮检对的模拟地震记录与观测记录之间的残差，提高了与目的层有关的地震记录在目标函数中的权重，重点反演深部目标体的构造参数信息.对于二维观测系统中的任意炮检对(x _s，x _r)，其对地下介质的照明强度为

其中，I_s和I_r分别是单炮震源 x _s对空间某一点(x，z) 的照明强度及其散射到检波点 x _r处的能量，I表示(x，z)处的双向照明强度，G和G′分别表示下行波场和上行波场.对于弹性波多波照明，选用波场的水平分量和垂直分量之和作为总的波场值可以使照明效果达到最优.应用数学统计分析法，对于给定的目标层或目标范围，可以计算该炮检对关于该目的层的总照明强度，此处我们提出用归一化地震波场的总照明强度来定义观测系统的可视性，公式为

其中，vis表示地震波由 x _s处的震源激发 x _r处接收 时成像目标范围内整体的照明强度，即炮检对 (x _s，x _r)的可视性.

为便于估计任意一炮检对所得到的地震记录在波形反演计算时波场匹配过程中的作用，我们用可视性的平均值作为判别标准，大于平均值说明该炮检对接收到的地震记录与目的层的相关性较高，因此经预处理后可以得到该观测系统对目的层的相对可视性vis′为

其中，n_s和n_r分别为炮点和检波点的个数.φ为vis的平均值的倒数.当相对可视性的值大于1时，即该炮检对的可视性大于观测系统对目的层的可视性平均值，将其引入到波形反演中目标函数的计算，可以得到基于可视性的弹性波全波形反演的目标函数为

此时，与目标体相关性较高的地震数据在波场匹配时将会占有更多的权重.

图 1a和c分别为观测系统中某炮在MarmousiII 模型及反演初始模型上正演的水平分量地震记录，图 1b和d为相应的垂直分量地震记录.该炮对应的所有检波点对目标体的相对可视性曲线如图 2所示，该观测系统中共有400个检波器.可以看出，位于第120到190号和第280到400号之间的检波器与目标体的相关性较高，可以接收到更多来自于反演目标体的波场信息.以此可视性作用于弹性波全波形反演中波场水平分量和垂直分量的残差(图 3a和c，对应的位置残差能量增强(图 3b和d)可以使弹性波全波形反演计算时充分利用更多的有效数据，提高反演精度和稳定性.

图 1

Fig. 1

图 1 MarmousiII模型上正演的弹性波场(a)x分量、(b)z分量和初始模型上正演的弹性波场(c)x分量、(d)z分量 Fig. 1 Elastic wavefields from forward modeling on MarmousiII model(a and b) and initial model(c and d). (a) and (c)are x component.(b) and (d)are z component

图 2

Fig. 2

图 2 观测系统中某炮及其对应所有检波点的相对可视性曲线 Fig. 2 Visibility curve of one shot and all the receivers in acquisition system

图 3

Fig. 3

图 3 x分量观测与模拟记录的(a)原始残差和(b)经能量加权后的残差以及z分量观测与模拟记录的(c)原始残差和(d)经能量加权后的残差 Fig. 3 xcomponent wavefield residuals between the observed and calculated data:(a)initial residuals，(b)residuals weighted by energy. z component wavefield residuals between the observed and calculated data:(c)initial residuals，(d)residuals weighted by energy

反射能量较强的地震记录在波形反演数据拟合时占有更大的比重，而能量较弱的波场对应的速度结构在反演过程中得不到较好的恢复，即使其速度结构与真实模型相差很大，因此，这种波场能量的不均匀分布会影响反演结果的稳定性.两种原因会导致地震波能量在地下介质中不均匀分布，一是由于地震波的几何扩散现象，波场能量随着传播深度的增加而减小；二是由于地下介质的复杂性，使得波在其传播方向上存在照明阴影区.利用波动方程理论得到的双向照明强度自适应加权梯度，可以同时补偿因观测系统中偏移距不足导致深部波场能量弱的问题和复杂介质中存在照明盲区而引起的能量损失，加权因子取决于在地震波的传播方向对应位置 处的照明强度，这在数学上类似于利用海森矩阵优 化梯度，得到牛顿类优化方法中的模型更新量.地震波在弹性介质分界面处的能量分配情况可以用Zoeppritz方程描述. 波场由高速介质向低速介质传播时，透射系数T>1，地震波的透射能量强于反射能量，界面下部的照明强度高于界面上部，此时可以用双向照明强度反比例加权于相应参数的梯度；相反，地震波场由低速向高速介质中传播时，透射系数T ＜1，则地震波的反射能量强于透射能量，界面上部的照明强度高于界面下部，此时可以用归一化后的照明强度正比例加权于相应参数的梯度.由此可得到弹性波反演中经自适应能量补偿后的模型更新量δ m ′为

其中， I 是弹性波场照明强度.

采用MarmousiII模型(图 4)来测试本文的反演方 法.模型网格大小为174×500，网格间距为20 m，横向长10 km，纵向深3.48 km.观测系统模拟水中激发纵波，以海底多波观测方式采集接收纵波和转换波，100个震源水平分布于海水面以下40 m处，400个三分量的检波器位于海底.为了分析观测系统的可视性，设目的层为位于深部断层下的含油藏的圈闭构造，如图 4中黑框所示.采样率为2 ms，地震记录长度为6s.震源为主频10 Hz的雷克子波.反演选用LBFGS优化方法，最大迭代次数为100次.

图 4

Fig. 4

图 4 MarmousiII模型(a)vp;(b)vs;(c)ρ和线性初始模型;(d)vp;(e)vs;(f)ρ. Fig. 4 MarmousiII model(a)vp;(b)vs;(c)ρ Initial model;(d)vp;(e)vs;(f)ρ.

对于给定的目标范围，经照明统计可以得到该观测系统对地下介质的相对可视性如图 5a所示，图中横纵坐标分别代表该观测系统中检波器和震源的 序号，各点的数值代表当地震波由纵坐标对应位置 处的震源激发，横坐标对应位置处的检波器接收时，该炮检对记录对目的层的相对可视性.可以看出该观测系统中，第40到第80个震源，第100到第230个检波器，以及第280到第400个检波器所对应的 相对可视性均大于1，说明这个范围内的炮检对接 收到的波场对目的层的成像起到主要作用.图 5b为该观测系统对地下介质的弹性波场照明强度，由于观测系统中震源和检波器主要分布在模型中部，因此模型两侧的照明较弱，中上部照明较强.此外由于地震波的几何扩散现象以及复杂介质引起的照明不均，模型深部以及复杂断裂带下方的照明能量较弱.我们应用可视性分析方法提高对目的层成像作用较大的炮检点对应的水平分量和垂直分量的残差在波场匹配中的权重，同时采用自适应算法进行能量补偿.

图 5

Fig. 5

图 5 (a)观测系统对指定目的层的可视性和(b)弹性波场照明强度 Fig. 5 (a)Visibility of the acquisition system for specified target zone and (b)elastic wavefield illumination intensity

图 6为两种反演方法均迭代100次后的反演结果.纵波速度、横波速度和密度三种参数均得到了较好的恢复，可以看出，基于可视性分析和能量补偿的EFWI的反演结果要明显优于常规EFWI.从反演 结果、真实模型和初始模型中提取模型中水平距离为5.2 km处的参数变化曲线，如图 7所示.相比于常规EFWI，本文提出的EFWI方法反演的纵波速度、横波速度和密度曲线更加趋近于真实模型，对复杂构造的有效重构同时加快了反演的收敛速度，使本文方法在相同的迭代次数条件下获得更加精确的反演结果.

图 6

Fig. 6

图 6 反演结果 常规EFWI:(a)vp;(b)vs;(c)ρ和基于可视性分析和能量补偿的EFWI;(d)vp;(e)vs;(f)ρ. Fig. 6 Inversion results Conventional EFWI:(a)vp;(b)vs;(c)ρ. EFWI with visibility analysis and energy compensation:(d)vp;(e)vs;(f)ρ.

图 7

Fig. 7

图 7 不同方法反演的参数剖面对比(x=5.2 km处)：真实模型(黑实线)、初始模型(黑虚线)、常规EFWI反演结果(蓝线)、基于可视性分析和能量补偿的EFWI反演结果(红线) Fig. 7 Comparison of parameter(P-wave，S-wave and density)profiles from models recovered by conventional EFWI(blue line) and EFWI with visibility analysis and energy compensation method(red line)，initial models (black dashed line) and true models(black solid line)at the location of x=5.2 km.

对安徽省庐枞盆地某矿床的复杂地球物理模型(图 8)采用弹性波全波形反演方法重建速度和密度多参数模型.该模型总体呈复杂的似层状、平缓透镜状，空间上表现为穹隆状，中心以侵染状磁铁矿为主，富、厚矿多环于四周.矿床由多个矿体组成，其中规模较大的只有3个矿体，其余均为小矿体(吕庆田等，2010；廉玉广等，2011).对该模型进行重采样后得到的纵波速度模型如图 9a所示.模型网格点数为120×400，网格间距为20 m，模型横向长8 km，纵向深2.4 km.主矿体埋深距地表 1.3~1.8 km间，且厚度变化较大.横波速度和密度可以在纵波速度基础上由经验公式得出(Brocher，2005; Xu and McMechan，2014)，如图 9b和c所示，公式为

图 8

Fig. 8

图 8 S1线过罗河矿体地球物理模型(据廉玉广等，2011) Fig. 8 The complex geophysical model for Luohe deposits along S1 seismic line in Luzong Basin(Lian et al.，2011)

应用于弹性波全波形反演的初始模型如图 9d、e和f所示.采样率为2 ms，地震记录长度为2 s.震源子波为主频10 Hz的雷克子波.反演选用LBFGS 优化方法，最大迭代次数为100次.观测系统中100 个震源水平分布于地表下40 m处，起始震源位于水平距离40 m处，震源间距为80 m.400个检波 器深度均为20 m，起始坐标和间距分别为0 m和20 m.

图 9

Fig. 9

图 9 金属矿多参数模型 (a)vp;(b)vs;(c)ρ和初始模型;(d)vp;(e)vs;(f)ρ. Fig. 9 Multi-parameter models for metallic deposit (a)vp;(b)vs;(c)ρ Initial models;(d)vp;(e)vs;(f)ρ.

选定的目标范围即为矿体所在的位置，对该观测系统进行可视性分析的结果如图 10a所示，弹性波场照明强度如图 10b所示，弹性波的能量除了在地表较强，在矿体部位也有所聚集.将照明分析的结 果引入到全波形反演中，得到基于地震照明的弹性 波全波形反演结果如图 11d、e和f所示.与常规弹性波全波形反演的结果(图 11a、b和c)对比可以看出，矿体轮廓的反演更加清晰，说明在反演过程中引入可视性分析和照明补偿，可以有效地改善反演效果.

图 10

Fig. 10

图 10 (a)观测系统对指定矿体的可视性和(b)弹性波场照明强度 Fig. 10 (a)Visibility of the acquisition system for specific metallic deposit，(b)Elastic wavefield illumination intensity

图 11

Fig. 11

图 11 多参数模型反演结果常规EFWI:(a)vp;(b)vs;(c)ρ和基于可视性分析和能量补偿的EFWI;(d)vp;(e)vs;(f)ρ. Fig. 11 Inversion results of muti-parameter models for metallic deposit Conventional EFWI:(a)vp;(b)vs;(c)ρ. EFWI with visibility analysis and energy compensation:(d)vp;(e)vs;(f)ρ.

(1)基于可视性分析的弹性波全波形反演策略可以适用于金属矿等复杂介质构造的多参数建模.通过照明统计分析可以得到给定观测系统对反演目的层的可视性，将其应用于波场匹配过程中对应炮检点的残差能量计算.构建基于可视性的反演目标函数，可以充分利用来自于目的层的有效地震信息，提高反演精度.

(2)由于观测系统对地下介质的照明强度可以直观地描述地震波能量的分布情况，利用Zoeppritz方程将弹性波场照明强度引入到全波形反演中，使得地震波场的能量自适应地均衡分布和优化梯度，改善反演质量.

(3)弹性波全波形反演方法可以高精度地重构深部金属矿地质体的多参数模型，包括纵横波速度和密度，可以实现金属矿区地下介质构造的高精度偏移成像.由于弹性波全波形反演可以利用多波多分量全波场信息，克服了单一的纵波地震勘探时由于弱阻抗差可能带来的不利影响，为有效地应用地震勘探技术开展深部矿产资源探测提供了新途径.

致谢 感谢国家自然科学基金(41374115)和国家科技重大专项(2011ZX0502500104)的支持，感谢评审专家的宝贵意见，在论文研究过程得到了王者江教授、巩向博副教授、张凤蛟博士的有益帮助，在此一并致谢.