1 引 言

AVO 反演技术已被广泛应用在石油勘探开发中^[1-8]，但目前的AVO 反演几乎都是利用Zoeppritz方程近似式实现的，勘探过程中所遇到的许多问题很难满足Zoeppritz方程近似的成立条件，即满足小角度和弱反射界面^[9-19]，在解决大角度和强反射反演问题时遇到了困难^[9-19]，现有的基于Zoeppritz方程近似的AVO 反演方法难以满足勘探需要.在勘探阶段地层密度信息十分缺乏，而地层密度是反映地层岩性的重要的参数之一；又因在地震波速度公式中含有介质的密度，使得由Zoeppritz方程表示的地震波反射系数与介质密度的关系更加复杂，利用反射系数反演介质密度的难度增加，而再利用Zoeppritz方程近似实现地层密度(一般不考虑波速度中介质密度的影响)反演，方法本身就存在大的偏差，因此基于Zoeppritz偏导方程(Zoeppritz方程对介质密度的导数方程)精确解，考虑波速度中介质密度的影响，开展包括适合于大角度反射问题的地层密度反演将具有实际意义.另一方面在实际的地震记录中包含了丰富的AVO 信息，这也正是地震勘探方法的优势所在，如何充分利用多角度(偏移距)信息仍是需要研究的新问题.目前的反演中很少同时考虑多角度反射波的联合反演问题，多角度反射波信息并没有得到充分利用.为减少反演的多解性、提高反演精度和收敛性^[20]，需要充分利用多角度信息实现多角度的联合反演.文献[21]给出了反射系数对介质密度比导数的计算方法，但尚未涉及利用地震波反射系数对介质密度导数(称反射系数梯度)精确解实现地层密度的反演.本文利用Zoeppritz偏导方程精确解，由反射系数对介质密度比的导数给出了计算反射系数梯度的方法，利用反射系数梯度精确解实现了地层介质密度的多角度联合反演，算例获得了好的反演效果，拓展了传统AVO 反演技术的应用范围，克服了常规利用AVO 反演技术计算地层介质密度所遇到的困难.

2 Zoeppritz偏导方程精确解

文献[21]已经给出了Zoeppritz方程对反射界面两侧介质密度比偏导的计算方法及矩阵元表示.由于介质密度的反演需要的是反射系数梯度，所以需要在Zoeppritz方程矩阵元对介质密度比偏导表示^[21]的基础上进一步给出地震波反射系数梯度的精确计算方法，通过下面的方法可实现反射系数梯度的精确计算.

设入射波为纵波(p波)，在反射界面，反射系数所满足的Zoeppritz方程为^{[1, 2]}

(1)

其中

(2)

在(2)式中R_pp、R_ps、T_pp、T_ps分别为反射p(纵)波、sv(横)波和折射p、sv波的反射和折射系数, T 表示矩阵的转置，a₁₁ =sinα，a₁₂ =cosβ，a₁₃ =-sinα′，a₁₄ =cosβ′，a₂₁ =cosα，a₂₂ =-sinβ，a₂₃ =cosα′，a₂₄ =sinβ′，a₃₁ =cos2β，a₃₂ =-$\frac{{{\upsilon }_{s1}}}{{{\upsilon }_{p1}}}$sin2β，a₃₃ =-$\frac{{{\rho }_{2}}}{{{\rho }_{1}}}\frac{{{\upsilon }_{p2}}}{{{\upsilon }_{p1}}}$cos2β′，a₃₄ =-$\frac{{{\rho }_{2}}}{{{\rho }_{1}}}\frac{{{\upsilon }_{s2}}}{{{\upsilon }_{p1}}}$sin2β′，a₄₁ =$\frac{\upsilon _{s1}^{2}}{{{\upsilon }_{p1}}}$sin2α，a₄₂ =vs1cos2β，a₄₃ =$\frac{{{\rho }_{2}}}{{{\rho }_{1}}}\frac{\upsilon _{s2}^{2}}{{{\upsilon }_{p2}}}$sin2α′，a₄₄ =-$\frac{{{\rho }_{2}}}{{{\rho }_{1}}}$v_s2cos2β′ ，b₁ =-sinα ，b₂ =cosα ，b₃ =-cos2β ，b₄ =$\frac{\upsilon _{s1}^{2}}{{{\upsilon }_{p1}}}$sin2α .α 为p波入射角，β为转换sv波(横波)的反射角，α′，β′ 为透射p,sv波的折射角，v_p1，v_s1，ρ₁ ，v_p2，v_s2，ρ₂ 分别为波在介质1、2 中的p、sv波的速度和介质的密度.

设在反射界面两侧介质的密度比为ρ₂₁ =ρ₂/ρ₁，将Zoeppritz方程的第4行同除以vp1，然后将方程两边ρ₁、ρ₂ 求导得

(3)

(4)

方程(3)、(4)为Zoeppritz偏导方程，其中^[21]

(5)

(6)

(7)

求解方程(1)、(3)、(4)可计算反射系数对介质密度的偏导数(反射系数梯度)，其中速度比利用了关系$\frac{{{\upsilon }_{p2}}}{{{\upsilon }_{p1}}}$=$\sqrt{\frac{1}{{{\rho }_{21}}}}$ξ₂₁，$\frac{{{\upsilon }_{s2}}}{{{\upsilon }_{p1}}}$=$\sqrt{\frac{1}{{{\rho }_{21}}}}$η₂₁，在此ξ₂₁，η₂₁ 为与反射界面两侧与介质弹性参数有关的常数.由此得介质波速度比对介质密度比的导数为

利用这两关系式，实现了波速度比对介质密度比导数的计算，在此并不需要计算ξ₂₁和η₂₁.

3 地层介质密度的多角度反演

设反演地下介质密度的目标函数为

(8)

其中D为实际地震记录，S(t)= W*R，为模型响应(合成地震记录)，W为地震子波，R是m的非线性函数.设f(m)=S-D，则f(m)为参差向量函数，若设m^k为m的第k次迭代近似.使反演的目标函数极小得到迭代公式

(9)

(10)

其中G_k=$\nabla $F(m)，λ^k为第k次迭代计算的步长.

在式(9)中，为提高计算的速度和稳定性把正定对角矩阵加到G_k^TG_k上，使矩阵变成条件数较好的对称正定矩阵，反演迭代公式修改为

(11)

其中I为n阶单位矩阵，α_k为一个正实常数.计算步骤如下:

步骤1 给定初始点m¹，初始参数α₁，增长因子β，允许误差ε>0，计算F(m¹)，置α=α₁，k=1；

步骤2 置α=α/β，计算G_k；

步骤3 解方程组:

(G_k^TG_k+α_kI)dk=G_k^Tf^k求得方向d_k后，令m^k+1=m^k+d^k；

步骤4 计算F(m^k+1)，若‖F(m^k+1)‖＜ ‖F(m^k+1)‖则转向步骤6，否则转向步骤5；

步骤5 若‖G_kf(m^k)‖ ＜ε，则停止迭代，得m=m^k；否则，置α =αβ，转向步骤3；

步骤6 若‖G_kf(m^k)‖ ＜ε，则停止迭代，得m=m^k+1；否则置k=k+1，返回步骤2.

4 地震子波与反射系数对介质密度导数的褶积

对(8)式取梯度有

(12)

其中ΔW=0.由(12)式可知在地层介质密度迭代反演方程(11)中，反演矩阵G的计算最终归结为地震子波与反射系数对介质密度导数的褶积，因此研究地震子波与反射系数对介质密度导数的褶积计算是本文反演方法的关键，该褶积反映了合成地震记录的振幅对地层介质密度变化反映的敏感程度.

假设地震子波振幅为W= [w₁，w₂，…wNW]，Nw 为子波的长度，则合成地震记录可写为

(13)

由(12)、(13)式得地震子波与多角度反射系数对介质密度导数的褶积可写为

(14)

其中S_n为第n层反射界面反射波的合成地震记录.

设第j层反射界面有Nα 个入射角度的地震记录，入射角分别为α₁…α_i…α_Nα，第i个入射角p-波的反射系数为R_j=R_pjαi，因每个反射界面波的反射系数仅与反射界面两侧的介质密度、波速度及入射角有关，与其他层的介质密度等参数无关，因此第j个反射界面反射系数对介质密度的导数只有对ρ_j、ρ_j+1 的导数不为零.在只有纵波地震记录的情况下，利用多角度反射系数对介质密度导数可构造如下反射系数梯度矩阵:

(15)

同时有纵波和横波地震记录，设在第j层反射界面，纵波反射系数为R_p，jαi，横波反射系数为R_s，jαi，总反射系数为R_j= [R_p，j_αi，R_s，jαi]^T，则可构造如下梯度矩阵:

(16)

5 地层介质密度反演算例

下面选择泥岩—油砂反射界面，地层参数:波速度v_p1=2500m/s，v_p2=2900m/s; v_s1=1020m/s，v_s2=1550m/s; 密度ρ₁=2200kg/m³，ρ₂=2300kg/m³.地震子波为Ricker 子波，频率f=35 Hz，采样间隔2ms，子波长度80ms.

表 1为地层介质密度反演迭代20次的计算结果，N_α 为反演角度的个数，每隔5°选择一个反演角度，分别选取了N_α =3~9 个角度进行了反演.当Nα=9则表示选取α=0°，5°，10°，15°，20°，25°，30°，35°，40°共9个入射角的地震记录反演地层的密度.图 1明显表明随着反演角度的增多地层介质密度反演结果的计算精度在逐步提高.表 1 中算例显示反演结果的前4位有效数字是相同的，说明层间波不存在相干时该地层密度反演方法是可以达到较高的计算精度.

图 2和表 2考察了反演介质密度的计算精度随反演迭代次数(k)的变化关系，列出了几种不同反演初值的反演结果.图 2(a~c)地层模型数值为:ρ_1,mod=2200kg/m³，ρ_2,mod=2300kg/m³，图 2a的反演初值ρ_1,sta=ρ_2,sta=2250kg/m³，取反射界面两侧介质密度的平均值，图 2b的反演初值ρ_1,sta=ρ_2,sta=2240kg/m³，偏向接近第一层介质的密度，图 2c的反演初值ρ_1,sta=ρ_2,sta=2270kg/m³，接近第二层的密度；图 2(d~f)为改变地层模型密度的计算结果(图 2d为ρ_1,mod=2200kg/m³，ρ_2,mod=2400kg/m³；图 2e为ρ_1,mod=2300kg/m³，ρ_2,mod=2000kg/m³；图 2f为ρ_1,mod =2300kg/m³，ρ_2,mod=1800kg/m³)，在这三种情况下反演初值均取反射界面两侧介质密度的平均值.图 3考察了当两层介质密度初值取不同数值时对反演结果的影响，图 3 中所有模型的密度均为ρ_1,mod=2300kg/m³，ρ_2,mod=1800kg/m³，图 3(a~d)为选取不同反演初值的计算结果(图 3a为ρ_1,sta=2200kg/m³，ρ_2,sta=1700kg/m³；图 3b 为ρ_1,sta =2250 kg/m³，ρ_2,sta= 1900kg/m³；图 3c为ρ_1,sta =2400kg/m³，ρ_2,sta=1700kg/m³；图 3d 为ρ_1,sta =2350 kg/m³，ρ_2,sta=1900kg/m³).图 2、图 3 的所有计算都迭代了20次，均显示在这些迭代的计算过程中，4~8次之间的迭代校正量较大，当迭代次数超过9次后，反演结果趋向于一个较为稳定的数值.这两个算例还显示:反演初值对计算结果确实影响较大，当反射界面两侧介质密度初值取相同数值时(ρ_1,sta=ρ_2,sta)，如果初值能取两介质密度的平均值，计算结果能收敛到模型参数值(图 2a,d~f)，但当初值不是取两介质密度的平均值时(注意两个初值仍取ρ_1,sta=ρ_2,sta)，反演结果不一定收敛到地层模型值，而是趋近于一个确定的常数值，这表现出了反演的多解性(图 2(b,c))；当反射界面两侧介质密度初值取不相同数值时(ρ_1,sta≠ρ_2,sta，见图 3)，如果取反演初值分别接近两层介质的模型数值，反演结果也能较好地收敛到模型的数值.由表 2(对应于图 2a)还可以看出在迭代次数10~20之间，计算反演结果有一个小的波动，其计算结果收敛于一个固定的常数值.这两个算例从第11次迭代起目标函数模的平方和均已小于10^-12.

图 4和表 3 的计算设定了4 个反射界面(5 个介质层)，层间反射波存在相互干涉，地层的层速度依次为vp=1980，2200，2440，2660，1980m/s; vs=808，1260，1304，1478，808m/s.ρ_m、ρ_st、ρ_inv分别为地层模型(介质密度)、反演初值和反演结果.图 4地层模型存在4个反射界面时的反演结果，不同反射界面间的反射波存在相互干涉(层间距小于波的相干长度)，(两相邻反射界面之间分别存在10个采样点)，各不同层段反演结果能较好地逼近地层模型参数值，准确地计算结果见表 3的ρ_inv1，最大相对误差在1%之内.

表 3的ρ_inv2为假设地层反射界面共有4 个，在第一个反射界面上方为均匀介质，第4 个反射界面的下方也是均匀介质，并且这4 个连续反射界面的间距较小(即这4 个连续反射层的层间距取为采样间隔)时的计算结果，参数取值ρ_inv1完全相同(见表 3)，尽管反演初值与地层模型参数的差距有的反射层可能较大，但经过10次反演迭代后地层介质密度的反演结果与模型介质密度的最大相对误差在1% 之内，说明该方法同样适应于反演层间反射波存在多次干涉的情况，且具有较高的计算精度.

6 结 论

本文给出了由地震波反射系数对介质密度比导数(公式(5))计算反射系数梯度的计算方法(公式(3)、(4))，实现了基于Zoeppritz偏导方程精确解的多角度联合反演介质密度的计算(利用公式(3)、(4)、(11)和(15)).通过算例对计算方法进行了考察，结果显示:在反演条件相同的情况下无层间干涉的反演结果比存在层间相干时要高，但迭代10次后计算结果的最大相对误差能够控制在1% 之内；多角度反演结果随反演角度的增加而提高；计算反射系数及解Zoeppritz偏导方程采用的是方程组的直接解法，在求解这些参数时的计算机耗时较少.由于本方法实现反射系数梯度精确解的计算并没有附加任何新的条件，并充分考虑了波速度中介质密度的影响，所以方法的适应范围是全角度和不同变化的地层(即Zoeppritz方程成立的条件)，克服了常规VAO 反演仅适应于小角度弱反射的问题，使方法不受反演角度及界面反射强弱的限制.利用地震波反射系数梯度精确解多角度反演地层密度的实现，使反演能够充分利用多角度的地震记录信息，尤其是大反射角的信息，这对于实现高精度地震资料的反演将具有重要意义.