地震波全波形反演(FWI)是近30年来发展起来的高分辨率地震数据处理方法.它利用地震观测记录的波形信息提取地下介质的弹性参数(如速度、密度和弹性模量等).全波形反演既可以为偏移成像提供精确的背景速度(刘国峰等，2012)，也可以应用到储层定量描述等领域(单蕊等，2011；李国平等，2011)，在油气勘探开发中具有重要作用.

由于全波形反演利用了地震波场的全部信息，理论上来讲，应该尽可能多的考虑影响波传播的因素，特别是将大孔径、宽方位数据用于全波形反演时，AVO/AVA效应和方位各向异性效应等需要特别考虑，这就要求将声波近似单参数全波形反演扩展到多参数全波形反演，包括纵横波速度、密度、吸收衰减效应、各向异性参数等.

然而，各参数之间的相互耦合，增加了多参数全波形反演的非线性程度.分析不同参数之间的相互关系，研究不同参数化方式对反演的影响程度，是建立多参数反演策略的有效途径.Tarantola(1986)和Mora(1987)指出，在P波震源反射地震波形反演中，弹性模量是一种比较差的参数化方式；对近偏移距数据(短波长成分)来说，阻抗是相对较好的参数化方式，而对远偏移距数据(长波长成分)来说，速度是比较好的参数化方式，但是，不管利用哪种数据信息，密度都不能被很好地估计.Kohn等(2012)指出，在弹性波多参数波形反演中，不管利用哪种参数化方式，弹性模量、纵横波速度、阻抗都可以被较好地反演出来，而密度反演结果较差；但总的来说，利用速度、阻抗参数化方式得到的反演结果要优于模 量参数化方式得到的反演结果.Forgues和Lambaré (1997)通过辐射模式分析指出，密度很难利用全波形反演得到.因此，大部分学者在反演过程中，将密度固定成某一常数或利用速度通过经验关系换算得到密度.

密度是岩石物性参数的重要组成部分，如果能通过反演同时得到比较可靠的纵横波速度、密度资料，这对储层评价、岩性解释和油藏描述等都具有重要作用.弹性波多参数全波形反演的研究实例证明，虽然可以得到相对可靠的纵横波速度，但密度并不能被很好地估计(Choi et al.，2008).因此，需要对密度反演进行深入研究，提出有利于密度反演的多参数全波形反演策略.Jeong等(2012)提出了一种频率域密度反演策略，首先将密度固定为任意常数，通过单参数反演得到弹性模量，进而将换算得到的纵横波速度作为初始模型，再同时反演速度和 密度，有效提高了密度的反演精度.Prieux等(2013) 先利用大偏移距数据单参数反演速度参数，再利用全偏移距数据同时反演速度、密度参数，有效提高了密度反演的稳定性，但密度反演结果仍然不够理想.他们同时指出，利用速度-密度参数化方式反演得到的速度、密度和阻抗结果要好于利用速度-阻抗参数化方式得到的反演结果.

本文利用速度-密度参数化方式，通过辐射模式及目标函数敏感度分析，研究了速度、密度在反演过程中的相互影响，提出了一种有利于速度、密度同时反演的策略，并通过球状异常模型和Marmousi-2模型的数值实验验证了该策略的有效性.

全波形反演通过拟合模拟数据u^mod和观测数据u^obs，使波场残差δu 达到最小，进而获取地下介质的物性参数m.目标函数可以表达为

选定某个初始矩阵B ⁰_k，对(3)式做m次递推可得：

由于模量-密度并不是一种好的参数化方式(Tarantola， 1986；Mora，1987；Kohn et al.，2012)，本文选取速度-密度为模型参数.目标函数对其它参数的梯度可以利用链式法则(Mora，1987)求得：

利用模量、速度、密度之间的关系：

本文的研究都基于公式(9)进行，在具体实验过程中，选取慢度(速度的倒数)和密度的倒数作为模型参数来更新模型(Mulder and Plessix，2004).

Forgues和Lambaré(1997)在Ray+Born近似条件下导出了变密度声波方程各种参数化方式对应的辐射模式的数学表达式，并通过对Hessian矩阵的特征值分解分析指出，密度反演比较困难.在Ray+Born近似条件下，速度、密度参数化方式对应的辐射模式如图 1所示.

图 1

Fig. 1

图 1速度-密度参数化方式对应的辐射模式 Fig. 1Radiation pattern as a function of diffraction angle for velocity and density

从图 1可以看出：(1)速度对应的辐射模式是各向同性的，密度只有在中小角度情况下散射能量，在大角度情况下几乎不散射能量，这说明速度对各种角度信息的数据都比较敏感，但密度只对中小角度的数据敏感；(2)在中小角度情况下，速度、密度的辐射模式基本一致，在这种情况下，速度、密度是耦合在一起的.因此当观测数据中只含有中小角度信息时，要将速度和密度完全区分开来是比较困难的(尽管此时可以得到较好的阻抗信息).

从上面的分析可以看出，辐射模式分析和目标函数敏感度分析的结论具有较好的一致性.如果直接利用全偏移距、大孔径反射波数据对速度、密度进行同时反演，可以得到速度的各种尺度的摄动成份和密度的高频摄动成份，但由速度的中长波长摄动所引起的目标函数的强烈非线性会严重影响密度的反演结果.前人的研究结果也表明，直接对反射地震数据进行速度、密度同时反演，虽然可以得到比较可靠的速度反演结果，但密度反演结果较差(Choi et al.，2008；Prieux et al.，2013).因此，在进行速度、密度反演的过程中，可以先利用大偏移距数据单参数反演得到速度的中长波长成份，然后再同时对速度、密度进行反演，这是一种比较合理的反演思路.但是，反演过程中无法消除密度的高频摄动对速度反演造成的影响.而且，在具体的实现过程中，速度、密度波长变化的复杂度无法与偏移距的大小建立明确的对应方式，并且反射地震数据中偏移距的大小本来就是一个相对概念，如何有效地将大、小偏移距的数据区分开来成为影响反演结果的重要因素.

因此，本文提出如下的分步反演策略.第一步，利用速度-密度方程同时反演速度和密度，这样可以得到比较可靠的速度反演结果，而密度反演结果较差.第二步，用第一步反演得到的速度作为初始速度模型，舍弃第一步反演得到的密度，利用最初给定的初始密度作为第二步的密度初始模型，继续进行双参数同时反演，这样可以同时得到相对可靠的速度、密度反演结果.为了进一步提高反演精度，可以将第二步反演得到的速度、密度结果作为初始模型，进行下一轮的速度、密度联合反演.

为了进一步说明速度、密度在反演过程中的相互影响，验证速度、密度分步联合反演策略的有效性，进行如下数值实验.

本节数值实验基于球状异常模型.球状异常体位于均匀背景模型中央，异常体半径100 m.速度模 型均匀背景值为3000 m·s^-1，球状异常值为3300 m·s^-1， 密度模型均匀背景值为2000 kg·m^-3，球状异常值为2200 kg·m^-3.采用四周观测方式，使得目标体在各个方向的照明均匀，进而消除观测系统的影响. 模型纵横向大小为2 km×2 km，网格离散间距10 m. 数值模拟152炮数据，记录长度1.5 s，时间采样间隔1 ms，震源函数选用主频7 Hz的Ricker子波.

首先，进行速度(v)、密度(ρ)单参数全波形反演.合成观测数据时，需要反演的模型参数中含有球状异常体，另外一个模型参数为均匀背景模型，初始模型为均匀背景模型.反演结果如图 2所示，图 3显示的是穿过球状异常体中心的纵向剖面.从图中可以看出，在此种观测方式下，模型各个方向照明均匀，速度、密度反演结果与真实模型非常接近.

图 2

Fig. 2

图 2单参数全波形反演结果 (a)速度；(b)密度. Fig. 2Mono-parameter full-waveform inversion results for velocity (a) and density (b).

图 3

Fig. 3

图 3穿过球状异常体中心的单参数反演结果纵向剖面 (a)速度；(b)密度. Fig. 3Vertical profiles across the center of the inclusion for (a) P-wave velocity and (b) density

接下来，进行速度、密度同时反演，以考察反演过程中速度、密度之间的相互影响.合成观测数据时，其中一个模型参数中含有球状异常体，另外一个模型参数为均匀背景模型.反演过程中所用的初始速度、密度模型为均匀背景模型.

图 4显示的是速度含有球状异常体，密度为均匀背景模型时双参数联合反演的结果，图 5为穿过球状异常体中心的纵向剖面.很显然，在反演过程中，由速度异常造成的数据残差同时反投影到了速度、密度反演结果当中.从图 4和图 5中可以看出，速度反演结果与真实模型非常接近，但密度反演结果却远远偏离于真实模型.

图 4

Fig. 4

图 4双参数全波形反演结果 (a)速度；(b)密度. Fig. 4Multi-parameter full-waveform inversion results for velocity (a) and density (b)

图 5

Fig. 5

图 5穿过球状异常体中心的双参数反演结果纵向剖面 (a)速度；(b)密度. Fig. 5Vertical profiles across the center of the inclusion for (a) P-wave velocity and (b) density

图 6显示的是密度含有球状异常体、速度为均匀背景模型时双参数联合反演的结果，图 7为穿过球状异常体中心的纵向剖面.可以发现，在反演过程中，由密度异常造成的数据残差也同时反投影到了速度、密度反演结果当中.从图 6和图 7中可以看出，速度、密度反演结果都与真实模型比较接近，但密度反演结果分辨率较低，还需进一步的反演.将图 6、7与图 4、5进行对比，可知速度对密度反演的影响较大，而密度对速度反演的影响相对较小.密度反演依赖于比较可靠的速度信息，而速度的准确反演也离不开有效的密度信息.

图 6

Fig. 6

图 6双参数全波形反演结果 (a)速度；(b)密度. Fig. 6Multi-parameter full-waveform inversion results for velocity (a) and density (b)

图 7

Fig. 7

图 7穿过球状异常体中心的双参数反演结果纵向剖面 (a)速度；(b)密度. Fig. 7Vertical profiles across the center of the inclusion for (a) P-wave velocity and (b) density

本节数值实验基于Marmousi-2模型(Martin et al.，2006)，横波速度为0，仅反演纵波速度与密度.模型纵横向大小为8 km×3.5 km，网格离散间距20 m.数值模拟79炮数据，第一炮位于地表水平方向100 m处，炮间距100 m，检波点分布于地表所有网格点上，数据记录长度5.7 s，时间采样间隔2 ms，震源函数选用主频7 Hz的Ricker子波.为了提高计算效率，实验过程中利用了相位编码同步震源技术(Krebs et al.，2009).为了降低同步震源相位编码技术造成的干涉噪声，需要改变每个频率组循环中相位编码的随机数(Krebs et al.，2009)，以尽量减少同步震源波场的相干性.

首先将密度固定，研究密度在速度反演中所起的作用.真实速度模型和高斯平滑初始速度模型如图 8所示.

图 8

Fig. 8

图 8(a)真实速度；(b)高斯平滑初始速度 Fig. 8(a) The true velocity model; (b) the Gauss-smoothed initial velocity

分别利用以下四种不同密度模型进行单参数(速度)全波形反演：

(1)常密度，密度值为水层密度1010 kg·m^-3，如图 9a所示；

(2)水层密度1010 kg·m^-3，水层以下密度值为 2000 kg·m^-3，如图 9b所示；

(3)高斯平滑初始密度，如图 9c所示；

(4)真实密度，如图 9d所示；

图 9

Fig. 9

图 9密度模型 (a)常密度，密度值为水层密度1010 kg·m-3；(b)水层密度1010 kg·m-3，水层以下密度值 为2000 kg·m-3；(c)高斯平滑初始密度模型；(d)真实密度模型. Fig. 9The density model (a) Constant density model which equals 1010 kg·m-3; (b) In sediments, density equals 2000 kg·m-3, in water,density equals 1010 kg·m-3; (c) The Gauss-smoothed initial density model; (d) True density model.

速度反演结果如图 10所示.由图 10a可以看出，当密度信息完全不准时，速度反演结果中有许多高频假象，这是因为进行单参数反演时，由密度变化引起的波场残差全部投影到速度反演结果当中，导致速度反演结果较差，特别是有强阻抗界面存在时，这种影响更加明显.由图 10(b—d)可以看出，考虑 强阻抗界面的影响时，由密度变化引起的波场残差对速度反演的影响并不是特别明显，但这种影响会导致速度反演结果在某些局部地方过高或过低，进而影响速度的后续使用.图 11显示了图 10(b—d)的纵向速度剖面，可以更清楚的看到这种影响.由此我们可以得知，利用速度、密度参数化方式单参数反演得到的速度结果并不能满足需要，特别是将其用于后续的多参数全波形反演.

图 10

Fig. 10

图 10单参数速度反演结果 (a)利用图 9a密度模型单参数反演得到的速度结果；(b)利用图 9b密度模型单参数反演得到的速度结果； (c)利用图 9c密度模型单参数反演得到的密度结果；(d)利用真实密度模型单参数反演得到的速度结果. Fig. 10Mono-parameter full-waveform inversion results for velocity (a) The inverted velocity using the density model shown in Fig.9a; (b) The inverted velocity using the density model shown in Fig.9b; (c) The inverted density using the density model shown in Fig.9c; (d) The inverted velocity using the true density model as shown in Fig.9d.

图 11

Fig. 11

图 11水平方向3 km(a)、4 km(b)、5 km(c)处的纵向速度剖面 黑线代表真实速度模型，红线代表初始速度模型，蓝线代表用图 9b密度模型单参数反演得到的速度结果，绿线代表 用图 9c密度模型单参数反演得到的速度结果，粉线代表用真实密度模型单参数反演得到的速度结果. Fig. 11Vertical profiles at horizontal positions of 3 km (a), 4 km (b) and 5 km (c) of the inverted results shown in Fig.10 The black line indicates the true velocity, the red line denotes the initial velocity, the blue line indicates the reconstructed velocity using the density shown in Fig.9b, the green line denotes the reconstructed velocity using the density shown in Fig.9c, the pink line denotes the reconstructed velocity using the true density.

将速度固定，研究速度对密度反演的影响程度.真实密度模型和高斯平滑初始密度模型如图 12所示.

分别利用图 10(b—d)速度模型与图 8a所示真实速度模型进行单参数(密度)全波形反演，密度反 演结果如图 13所示.从图中可以看出，当速度信息准确时，可以得到比较准确的密度反演结果.但是，速度模型的微小变化会引起密度反演结果的巨大变化，并且有可能导致密度反演结果远远偏离真实情况，这说明密度反演依赖于相对可靠的速度信息.图 14显示了图 13(b—d)的纵向剖面.

由上面的分析可知，密度反演依赖于比较可靠的速度信息，而速度的准确反演也离不开有效的密度信息，单参数反演得到的速度结果并不能用于重建密度模型，需要进行速度、密度分步联合反演.接下来通过如下数值实验，证实本文第3节提出的反演策略的有效性.

图 12

Fig. 12

图 12(a)真实密度；(b)高斯平滑密度 Fig. 12(a) The true density model; (b) the Gauss-smoothed initial density

图 13

Fig. 13

图 13单参数密度反演结果 (a)利用图 10b速度模型单参数反演得到的密度结果；(b)利用图 10c速度模型单参数反演得到的密度结果； (c)利用图 10d速度模型单参数反演得到的密度结果；(d)利用真实速度模型单参数反演得到的密度结果. Fig. 13Mono-parameter full-waveform inversion results for density (a) The inverted density using the velocity model shown in Fig.10b; (b) The inverted density using the velocity model shown in Fig.10c; (c) The inverted density using the velocity model shown in Fig.10d; (d) The inverted density using the true velocity model.

图 14

Fig. 14

图 14水平方向3 km(a)、4 km(b)、5 km(c)处的纵向密度剖面 黑线代表真实密度模型，红线代表初始密度模型，绿线代表用图 10c速度模型单参数反演得到的密度结果，蓝线代表 用图 10d速度模型单参数反演得到的密度结果，粉线代表用真实速度模型单参数反演得到的密度结果. Fig. 14Vertical profiles at horizontal positions of 3 km (a), 4 km (b) and 5 km (c) of the models shown in Fig.13 The black line indicates the true density, the red line denotes the initial density, the green line denotes the reconstructed density using the velocity shown in Fig.10c, the blue line indicates the reconstructed density using the velocity shown in Fig.10d, the pink line denotes the reconstructed density using the true velocity.

本节数值实验基于球状异常模型，模型参数和观测系统设置与4.1.1节完全相同.合成观测数据时，真实速度、密度模型中均含有球状异常体.

首先，利用均匀背景模型作为初始速度、密度模型，对速度、密度进行同时反演.反演结果如图 15所示，图 16显示了反演结果的纵向剖面.从图中可以看到，速度反演结果与真实模型具有较好的一致性，但密度反演结果却远偏离于真实模型，误差较大.

图 15

Fig. 15

图 15第一阶段双参数全波形反演结果 (a)速度；(b)密度. Fig. 15The inverted results in the first stage for velocity (a) and density (b)

图 16

Fig. 16

图 16穿过球状异常体中心的第一阶段双参数反演结果纵向剖面 (a)速度；(b)密度. Fig. 16Vertical profiles across the center of the inclusion for (a) P-wave velocity and (b) density

接下来，将第一步反演得到的速度作为第二阶段的初始速度模型，初始密度模型依然利用均匀背景模型，再一次进行速度、密度联合反演.反演结果和纵向剖面如图 17和18所示.从图中可以看到，经过第二阶段的联合反演，可以得到相对比较准确的 速度、密度反演结果，但密度反演结果的分辨率不高.

图 17

Fig. 17

图 17第二阶段双参数全波形反演结果 (a)速度；(b)密度. Fig. 17The inverted results in the second stage for velocity (a) and density (b)

图 18

Fig. 18

图 18 穿过球状异常体中心的第二阶段双参数反演结果纵向剖面 (a)速度；(b)密度. Fig. 18Vertical profiles across the center of the inclusion for (a) P-wave velocity and (b) density

为了进一步提高密度反演精度，利用第二阶段反演得到的速度、密度作为初始模型，进行第三阶段的双参数联合反演.反演结果和纵向剖面如图 19和20所示，不难发现，密度反演结果的分辨率较第二阶段有了进一步的提升.

图 19

Fig. 19

图 19第三阶段双参数全波形反演结果 (a)速度；(b)密度. Fig. 19The inverted results in the third stage for velocity (a) and density (b)

图 20

Fig. 20

图 20穿过球状异常体中心的第三阶段双参数反演结果纵向剖面 (a)速度；(b)密度. Fig. 20Vertical profiles across the center of the inclusion for (a) P-wave velocity and (b) density

本节数值实验基于Marmousi-2模型(Martin et al.，2006)，模型参数和观测系统设置与4.1.2节完全相同.为了提高计算效率，实验过程中利用了相位编码同步震源技术(Krebs et al.，2009).为了降低同步震源相位编码技术造成的干涉噪声，需要改变每个频率组循环中相位编码的随机数(Krebs et al.，2009)，以尽量减少同步震源波场的相干性.

真实速度、密度模型如图 8a和图 12a所示，高斯平滑初始速度、密度模型如图 8b和图 12b所示.首先，利用给定的初始速度、密度模型，对速度、密度进行同时反演.反演结果如图 21所示，图 22显示了反演结果的纵向剖面.从图中可以看到，虽然密度结果远偏离于真实模型，但速度反演结果与真实模型具有较好的一致性.

图 21

Fig. 21

图 21第一阶段双参数反演结果 (a)速度；(b)密度. Fig. 21The inverted results in the first stage：(a) P-wave and (b) density

图 22

Fig. 22

图 22第一阶段双参数反演结果纵向剖面，分别位于水平方向3 km（a1,b1）、4 km（a2,b2）和5 km（a3,b3）处的速度反演结果（A）与密度反演结果（B） 黑线代表真实模型，红线代表初始模型，绿线代表反演结果. Fig. 22Vertical profiles at horizontal positions of 3 km (a1,b1), 4 km (a2,b2) and 5 km (a3,b3) of the models shown in Fig.21 for P-wave velocity (A) and density (B) The black line indicates the true model, the red line denotes the initial model and the green line denotes the reconstructed model.

接下来，将第一步反演得到的速度作为第二阶段的初始速度模型，初始密度模型依然利用最初给 定的初始模型，再一次进行速度、密度联合反演.反 演结果和纵向剖面如图 23和24所示.从图中可以看到，经过第二阶段的联合反演，可以得到比较准确的速度、密度反演结果，并且第二阶段速度反演结果的分辨率有进一步的提升.

图 23

Fig. 23

图 23第二阶段双参数反演结果 （a）速度；（b）密度. Fig. 23he inverted results in the second stage: (a) velocity; (b) density

图 24

Fig. 24

图 24第二阶段双参数反演结果纵向剖面，分别位于水平方向3 km（a1,b1）、4 km（a2,b2）和5 km（a3,b3）处的（A）速度反演结果与（B）密度反演结果 黑线代表真实模型，红线代表初始模型，绿线代表反演结果. Fig. 24Vertical profiles at horizontal positions of 3 km (a1,b1), 4 km (a2,b2) and 5 km (a3,b3) of the models shown in Fig.23 for (A) P-wave velocity and (B) density The black line indicates the true model, the red line denotes the initial model and the green line denotes the reconstructed model.

从图 21图 21—24可以看出，进行分步多参数全波形反演，可以得到比较好的速度、密度反演结果，接下来进行第三阶段的多参数联合反演，看是否可以进一步提高反演精度.此时，利用第二阶段反演得到的速度、密度结果作为第三阶段的初始模型.第三阶段的反演结果和纵向剖面如 图 25和图 26所示.不难发现，密度反演结果的分辨率较第二阶段有了进一步的提升，但提升空间有限，主要是由于第二阶段密度反演结果已非常接近真实模型.

反演策略对多参数全波形反演至关重要.本文通过数值实验证实，密度反演依赖于比较可靠的速度信息，而速度的准确反演也离不开有效的密度信息.固定密度，单参数反演得到的速度往往并不能满足后续对速度的进一步要求.因此，本文提出了有利于速度、密度同时反演的策略.该策略充分考虑了多参数全波形反演中目标函数对速度、密度的敏感性差异，分步实现对速度、密度的反演，在一定程度上降低了多参数全波形反演的非线性程度.第一步，利用给定的初始模型对速度、密度进行同时反演.此时，可以得到相对可靠的速度反演结果，而密度反演结果较差.第二步，将第一步反演得到的速度作为第二步的初始速度模型，舍弃第一步反演得到的密度，利用最初给定的密度作为第二步的初始密度模型，进行双参数同时反演，这样可以同时得到比较可靠的速度、密度反演结果.为了进一步提高反演精度，可以将第二步反演得到的速度、密度作为初始模型，继续进行下一轮的双参数同时反演.数值实验证实了该策略的有效性.尽管本文是基于变密度声波方程研究多参数全波形反演策略，但该策略同样可以应用到弹性波或各向异性的多参数全波形反演中.