0 引言

全波形反演是利用数据拟合来估算地下参数的一种反演方法。Tarantola^[1]借鉴逆时偏移的思想，利用正传波场和反传波场的互相关进行迭代反演成像，开启了时间域全波形反演的先河。此后，Mora^[2]将其应用到弹性波时间域全波形反演。Pratt等^[3]将该思想应用到频率域全波形反演，在频率域用几个频率进行反演就可以得到较好的效果，引起广泛的兴趣。目前，全波形反演已经能够反演得到高分辨的模型参数，比如速度^[4]、密度^[5]、衰减^[6]和各向异性参数^[7-8]。这些模型参数可以描述地震波的传播，根据不同的问题选取合适的参数描述地下特征。这些模型参数可以通过单参数或者多参数全波形反演^[9-10]得到。

速度和密度在地震数据的定量解释和油气解释中扮演着重要角色。当速度不变、密度变化时也会引起反射率变化^[11]。因为全波形反演的目标函数考虑了地震数据的振幅和相位信息，而通常假设密度是常量，只反演速度，导致了反演的速度结果存在串扰^[12]。Plessix等^[9]做过一个研究测试，为了保持道集的平滑性，对密度进行了反演。但是，正如Forgues等^[13]所述，利用全波形反演得到密度很困难，因为很难将密度和速度的影响区分开来。速度扰动引起的波场变化在各个方向是一致的，而密度扰动引起的波场变化随着角度的变大而逐渐减小；所以，可以较好地重建速度成分，而只能重建密度的小尺度和部分中尺度成分。Jeong等^[14]提出了一种频率域全波形密度反演的分层策略，并将这种策略引用到声波频率全波形反演中。

本文对频率、初始速度模型、初始密度模型和密度速度同时反演四个方面对密度的影响进行了研究，详细展示了它们对密度反演的影响程度，并据此制定了比较稳定的反演策略：首先将密度固定，反演速度，此时波动方程中不再含有密度项，因而可以得到准确的速度模型；再将其作为初始速度模型进行速度、密度同时反演，可以较好地减小速度对密度的影响，从而得到稳定的密度模型。此后，用理论模型测试策略的有效性。

1 变密度声波全波形反演基本理论 1.1 频率域变密度声波方程正演计算

频率空间域二维声波波动方程为

式中：ρ是密度；U是波场；ω是角频率；v是速度；f是震源项；s是PML(完全匹配层)吸收边界参数，下标x，z分别表示x轴和z轴。

用最优化9点有限差分法将方程(1)离散，最终式(1)简单表示为

式中：U(ω)为频率域波场；f(ω)为频率域震源；S(ω)为阻抗矩阵。S(ω)是N×N(N为模型空间网格数，N= N_x×N_z)的大型稀疏矩阵，通常调用已有程序包来求解。目前存在的程序包包括UMFPACK^[15]，MUMPS，SUPERLU等，求解方法分为迭代法和直接分解法。在频率域中，对于多震源问题，直接分解法只需要将稀疏矩阵分解一次，其分解结果可以被各个震源利用，从而大大节约了计算量。本文采用了UMFPACK程序包中的直接分解法求解该大型稀疏方程。

1.2 频率域变密度声波全波形反演

正如大多数反演问题一样，我们最小化剩余误差函数。误差函数为

式中：E(m)为误差函数；m是模型参数；T表示转置；*表示共轭转置；δd为计算波场与观测波场的差值。

式中：d_i为第i个检波器的观测波场；u_i为模拟波场；n为检波器数。

我们通过误差函数的梯度方向来更新模型参数：

式中：k是迭代次数；α是步长，可以通过线性搜索等最优化方法得到；∇_mE为误差函数的梯度，可通过将式(3)两边对模型参数求导得到。

误差函数对速度和密度的梯度^[16]可分别表示为：

式中：U_f表示正传波场；U_b表示反传波场。梯度可以表示为正传波场和反传波场的互相关。速度和密度之间辐射模式的不同使梯度在具体表达方式上有所差异^[11]。

全波形反演遇到的挑战之一是解的不适定，它采用的是局部最优化算法，而其中较多的局部极小值容易使目标函数收敛到局部极值^[17-20]，导致反演结果不准确。针对这个问题，本文采用多尺度算法：先反演低频成分，然后将反演结果作为相邻较高频率反演的初始模型。因为低频能够反演模型的大尺度成分，高频反演小尺度成分，所以多尺度算法能够较好地减小局部极小值的影响，使目标函数较好地得到收敛，同时可以压制噪音和反演假象。

1.3 能量补偿

在全波形反演中，浅层地区能量较强，能够得到较好的成像；而随着深度的增加，由于透射损失和集合扩散等原因使能量减弱，极大地影响了深部成像效果。研究通过在梯度处理环节引入频率域补偿因子，并对其进行改造，使其对浅层梯度场进行噪音压制的同时对深层梯度场进行补偿，以提高深层成像效果。

梯度处理后的速度迭代更新公式为

补偿系数γ为

式中：f_s为频率；z为深度；C为常数。

2 密度影响因素

为了对频率、初始速度模型、初始密度模型和密度速度同时反演四个方面对密度的影响进行测试，本文对凹陷模型、SEG/EAGE推覆体模型和抽稀的Marmousi模型进行了密度反演，得到了一致的结论。受限于文章篇幅，本文只对于抽稀的Marmousi模型进行分析。抽稀的Marmousi模型水平坐标为0～3 110 m。地表放炮，共放311炮，炮间距10 m；地表接收，共311个检波器，检波器间距为10 m。

为了防止其他影响因素的干扰，在针对频率、初始速度模型、初始密度模型和密度速度同时反演的研究测试中，这四个方面之外的其他影响因素始终保证完全一致，保证了研究测试的科学性和可靠性。

2.1 频率

低频反演对应大尺度成分，高频反演对应小尺度成分。董良国等^[16]指出地震反射波对密度的较小尺度摄动具有很强的敏感度，但中长尺度的密度摄动对反射波基本没有影响。为了测试频率影响，本文用6个测试频率组进行单参数密度反演，频率组分别为5~9，10~14，15~19，20~24，25~29，30~34 Hz，每个频率组迭代20次。图 1a、b分别为真实密度模型和初始密度模型，图 2a、b分别为真实速度模型和初始速度模型。

图 3显示的是不同频率反演得到的密度模型对比。由图 3a和图 3b—f的对比结果可以看出，频率较低时反演的结果不准确；这是因为低频对应大尺度成分，而大尺度的密度摄动对反射波的振幅和时间走时都基本没有影响，因此无法利用低频反演密度大尺度成分。从图 3b—f可以看出，随着频率的增加，反演尺度越来越小，描述也越来越精细；但是一些较大尺度的结构却不能得到很好的反演，如图 3中黑色圆圈所示。因此，结合反演的初始模型合理选择频率是正确反演密度与提高反演精度的基础。

2.2 初始速度模型

多参数全波形反演的一个问题就是多参数带来的不稳定性，Forgues等^[13]指出速度和密度之间的相互耦合或者称为串扰加大了密度反演的难度。为了更清楚地认识速度对密度反演的影响，本文采取了不同的初始速度模型进行单参数密度反演。图 4a为偏离真实速度值较大的线性递增速度模型，简称偏离速度模型：

v=3 350+1.95z。

其中，v为速度，m/s。图 4b为线性递增速度模型：

v=1 850+1.95z。

光滑速度模型如图 2b所示。

本文采用多尺度算法，据2.1小节的反演结果，下文均选取3个频率组进行反演，分别为10~14，15~19，20~24 Hz，每个频率组迭代20次，低频反演结果作为相邻高频的初始模型。

图 5显示的是不同初始速度模型反演得到的密度模型对比(初始密度模型为图 1b)。图 5a得到错误的反演结果，可见速度对密度反演影响很大，一旦速度模型偏离真实模型较大，密度亦无法得到正确结果。而随着初始速度模型接近真实模型，反演得到的密度模型也越准确。可见，一个较为准确的初始速度模型是密度反演成败必不可少的因素。提高速度反演的准确度至关重要。

2.3 初始密度模型

全波形反演是反演精度高却“娇贵”的方法，因为它对初始模型的依赖性很强，尤其是多参数之间的相互影响更加剧了对初始模型的依赖。为了清楚地认识密度反演对初始模型的依赖，本文采取了不同的初始密度模型来进行单参数密度反演。图 6a为偏离真实密度值较大的线性递增密度模型，简称偏离密度模型：

v=2 500+0.724z。

图 6b为线性递增密度模型：

图 7显示的是由不同初始密度模型反演得到的密度模型对比(初始速度模型为图 2b)。图 7a得到的反演结果虽然反演出Marmousi模型的大致层位，但数值和真实模型差距较大。图 7b、c都能比较准确地得到反演结果。随机取出图 7b、c横向不同位置处的密度对比两种方法的效果，结果如图 8所示。由图 8可以看出，以光滑密度模型为初始密度模型的反演结果更接近实际模型。可见，随着初始密度模型接近真实模型，反演得到的密度模型也更加准确。

由图 5a和图 7a对比发现，初始速度模型对密度反演结果的影响甚大，甚至已经超过了初始密度模型对密度反演的影响。对比结果更加说明，得到较为准确的初始速度模型是密度反演成败的关键因素。

2.4 速度密度同时反演

Jianyong等^[14]指出，速度和密度的同时反演可以减少速度密度之间的耦合。为了更好地研究速度密度同时反演的效果，本文采用相同的速度和初始密度模型进行单参数密度反演和速度密度同时反演。

图 9为不同方法反演得到的密度模型对比(初始密度模型和初始速度模型分别为图 1b和图 2b)，图 10为图 9a、b横向不同位置处的密度对比。可以 看出：速度密度非同时反演和速度密度同时反演都能正确地得到反演结果；但是速度密度同时反演的结果更加稳定，与实际数据的数值更加接近。

速度密度同时反演减少速度密度之间的耦合，从而使反演更稳定；所以速度密度同时反演不失是提高密度反演精度的方法之一。

3 密度反演策略

通过以上的研究测试，本文提出了一种稳定有效的密度反演方法。首先，将密度固定，只反演速度。此时式(1)变为

由于式(10)中不再含有密度项，可以稳定反演得到速度。准确的速度模型是正确反演密度的基石。

然后，将反演得到的速度模型作为初始速度模型，进行速度和密度同时反演。本文简称该方法为交互。

为了验证本文提出的密度反演策略的有效性，进行三种实验方案的对比。方案一，交互反演；方案二，单参数密度反演；方案三，直接进行速度密度同时反演，只观察密度结果。本文对凹陷模型、SEG/EAGE推覆体模型和抽稀的Marmousi模型进行了反演策略测试，得到了一致的结论；受限于文章篇幅，笔者只对于抽稀的Marmousi模型进行研究分析。

图 11为不同方法反演得到的密度模型对比(初始密度模型和初始速度模型分别为图 1b和图 2b)，图 12为图 11b、c横向不同位置处的密度对比。由图 11可以看出，交互的反演结果更加稳定，没有出现像密度速度非同时反演和密度速度同时反演结果那样的异常值(如图 11a、b中绿色方框所示)，基本与实际模型完全相同，能够清晰准确地反演得到背斜下面深部的层位结构(如图 11c中绿色方框所示)。从图 12中更能直观地看到交互反演的反演结果与实际数据更接近。

测试表明，本文提出的全波形反演密度反演策略使密度反演更稳定、更准确。

4 结论与建议

对频率、初始速度模型、初始密度模型和密度速度同时反演四个方面对密度反演的影响进行了研究，结果表明：1) 频率的选择对密度反演至关重要，较低的频率得不到正确的反演结果，因为大尺度的密度扰动基本不会引起波场的变化；2) 初始速度模型对密度反演起着决定性作用，甚至已经超过了初始密度模型的影响，因此准确的反演速度是准确反演密度的基础和必需；3) 初始密度模型越接近实际模型，密度反演结果越准确；4) 密度速度同时反演会降低速度对密度反演的干扰，使得密度反演更加稳定。

根据研究测试结果制定了比较稳定的反演策略：先将密度固定，反演速度，此时波动方程中不再含有密度项，因而可以得到准确的速度模型；将其作为初始速度模型进行速度密度同时反演，可以较好地减小速度对密度的影响，从而得到稳定的密度模型。理论模型测试结果充分说明了策略的有效性。