在进行构造或沉积学研究时，地球介质的非均质性可分为大小两种尺度，其中大尺度非均质性描述的是介质的分层或者缓慢变化的特性，小尺度非均质性指附加在大尺度非均质性上的扰动(Aki，1969).由于小尺度非均质性的随机分布特征，常规的确定性建模方法无法对其进行完整描述.Ikelle等(1993)、姚姚和奚先(2002)、奚先和姚姚(2001，2005)等将地下介质视作随机变量，提出了应用统计学原理对地下介质的非均匀特性进行建模的方法，即随机介质建模方法.但是，要将该方法应用于地震勘探，需要事先了解地下介质的统计特性.

我们知道，当使用地震勘探方法进行沉积盆地研究时，沉积盆地中的地层界面为连续的强反射，而小尺度非均质性则为不连续的弱反射波.因此，如何定量地建立小尺度非均质性与不连续反射波的关系，并利用这个关系估计深部介质的非均质特征成为了研究人员关注的重点(Hurich，1996; Hurich et al., 2000).Gibson(1991)首先提出地下介质的横向非均质性可用叠后地震数据的横向自相关函数来近似.Poppeliers(2007)、秦臻(2008)建立了垂向上地下介质与叠后地震数据的函数关系式，并提出了应用反褶积方法估计纵向自相关函数及其特征参数的方法.虽然上述方法都曾被用于地震解释(Hurich et al., 2000; 秦臻，2008)，但都是从水平和垂向两个方向分别进行的，均是一维估计，因此对于二维或三维情况下的非均质性分析来说显然是有相当大的差距的.为了解决这个问题，Scholer 等(2010)后来建立了二维地下介质与叠后地震数据的函数关系式，并提出了应用蒙特卡洛反演法从叠后地震数据中同时估计纵、横向自相关函数的方法.但是该方法估计的自相关函数仍然是一维的，因此只能对复杂程度不高的二维介质进行估计.例如，当地层有倾角时，介质的扰动方向可能与水平方向成一定的夹角，这种情况下计算水平和垂向自相关函数，其结果必然是有较大误差的.为了解决上述问题，本文以褶积地震模型理论为基础，建立了叠后地震数据与随机介质波阻抗的功率谱关系，提出了二维平稳随机介质自相关函数的功率谱估计方法.该方法可同时估计二维自相关函数及其三个特征参数，即纵、横向自相关长度和自相关角度.以高斯型随机介质模型为例，本文分别针对不同特征参数的随机介质模型，以及不同主频的合成地震记录，对本文提出的方法进行了可行性验证，并针对随机介质参数估计中地震主频对估计结果的影响、模型的大小设置等问题进行了初步探讨.最后应用本文提出的方法，对实际地震数据进行了参数估计试验，并分析了不同成因相的介质随机特性. 2 随机介质的数学描述

随机介质模型是描述岩石物性在空间上具有随机变化特征的介质模型(Ikelle et al., 1993)，其基本特性可由统计量参数表征，包括均值、标准差和自相关函数.大尺度的非均质性用一阶统计量——均值来表征，小尺度的非均质性用二阶统计量——标准差和自相关函数来表征.以速度为例，随机介质速度模型的定义如下(Ikelle et al., 1993; 奚先和姚姚，2001)：

以上用于描述随机介质的统计量参数中，均值、标准差和自相关函数的类型可以直接通过声波、密度测井数据的统计得到(秦臻，2008)，不是本文讨论的重点；而本文中说提到的“随机介质参数”主要指的是自相关函数中的三个特征参数，即横向、纵向自相关长度a、b，以及自相关角度θ. 3 随机介质二维自相关函数的估计原理与流程 3.1 随机介质参数估计原理

对于时间域的叠后地震剖面，其一次反射波可用线性褶积模型表达，即地震道近似为反射系数与地震子波的褶积：

根据(1)式，设连续介质波阻抗由两部分组成：

该方法首先计算随机介质模型的功率谱，然后根据Wiener-Khintchine，由傅里叶反变换直接得到随机介质的自相关函数.与其他方法(Poppeliers，2007; Scholer et al., 2010)相比，该方法在估计自相关函数前既无需先求随机介质模型，也避免了最优化方法的搜索步骤，因此具有直观且高效的优点.

但是，实际工作中，在计算(5)式时，会遇到|S_f(ω)|²在某些频率处存在零值或接近于零值的情况，这时求得的|S_δZ(ω，k_x)|²值就变成无穷大，显然这是不合理的.所以，在实际工作时，可给|S_f(ω)|²加上一个白噪ε，使得求取过程变成一个稳定的过程.为此，我们把式(5)改写为

我们以一个实例来说明随机介质参数估计的具体流程.设有一个二维高斯型随机介质，其模型在时 间域(Time)和横向空间域(X)中的参数为： Nt=Nx=600，dt=1 ms，dx=1 m，速度均值=3000 m·s^-1，标准差=500 m·s^-1，自相关长度a=50 m，b=20 ms，自相关角度θ=30°，其自相关 函数如图 1a所示.假 设介质密度ρ(t，x)可由Gardner公式估计得到(Gardner et al., 1974)：

图 1

Fig. 1

图 1 二维随机介质的理论自相关函数(a)、波阻抗模型(b)、Ricker地震子波(c)和合成地震记录(d) Fig. 1 Theoretical autocorrelation function(a) and impedance model(b)of 2D r and om medium，Ricker wavelet(c)， and synthetic seismic records(d)

具体的估计步骤如下：

(1)从地震记录(图 1d)中估计出地震子波(图 2a);

图 2

Fig. 2

图 2 地震数据的子波估计结果(a)，二维汉明窗(b)和加窗后的地震数据(c) Fig. 2 Wavelet estimation from seismic data(a)，2D Hamming window(b) and windowed seismic records(c)

(2)为了避免地震数据的频谱能量向旁瓣泄漏，需要对地震数据做加窗(图 2b)预处理，本文采用二维汉明窗，处理结果如图 2c所示；

(3)计算子波的功率谱|S_f(ω)|²(图 3a)和地震数据的功率谱|S_s(ω，k_x)|²(图 3b)；

图 3

Fig. 3

图 3 子波功率谱(a)和地震数据的功率谱(b) Fig. 3 Power spectra of wavelet(a) and seismic records(b)

(4)根据(6)式计算随机介质的功率谱|S_δZ(ω，k_x)|²(图 4a)；

图 4

Fig. 4

图 4 估计得到的随机介质功率谱(a)和二维自相关函数(b) Fig. 4 Power spectrum of r and om medium(a) and 2D autocorrelation function(b)

(5)对随机介质的功率谱|S_δZ(ω，k_x)|²做傅里叶反变换得随机介质的二维自相关函数R(t′，x′)(图 4b).

(6)将图 4b中R(t′，x′)≥exp(－1)的点设为1，R(t′，x′)<exp －1 的点设为0(图 5a).图中取值为1的点构成一个椭圆；

图 5

Fig. 5

图 5 特征向量法估计自相关函数的特征参数步骤(1)的结果(a)和步骤(2)的结果(b) Fig. 5 Characteristic parameter estimation of autocorrelation function by eigenvector method in step(1)(a) and (2)(b)

(7)由图 5a中取值为1的M个点的坐标，计算方向矩阵，得

(8)提取图 5a中图形的边界(图 5b)；

(9)在边界上搜索与特征向量 α和β 相交的点，与O点距离最大的点设为A点，为横向自相关长度a，搜索与O点距离最小的点设为B点，为纵向自相关长度b(图 5b)；由此得到横向和纵向自相关长度分别为a=50 m和b=12 ms；

(10)计算自相关角度即OA与X轴的夹角θ=

以上流程中，为了方便，本文省去了时深转换的环节，所构建的随机介质为时间域的模型，其估计出的随机介质参数也是一个时间域中的相对值.要估计深度域的随机介质参数，理论上需要应用深度域地震剖面来进行参数估计工作. 4 二维平稳随机介质参数估计数值实验

为了进一步研究随机介质参数估计算法的可行性和稳定性，本文设计了三个试验，来研究不同随机介质参数模型、地震数据的主频和地震数据的规模对估计误差的影响. 4.1 试验1

本试验主要测试不同随机介质参数所构成的模型对估计误差的影响.为此构建了一系列不同自相关长度和角度的高斯型随机介质波阻抗模型(图 6)，用于测试估计结果的误差.这组模型大小Nx=Nt=300，网格间距dx=1 m，dt=1 ms，速度均值为5000 m·s^-1，标准差为500 m·s^-1，密度值根据Gardner公式估算.表 1列出了构建的随机介质模型自相关函数的特征参数，图 6(a1、b1、c1、d1)为对应的自相关函数图，图 6(a2、b2、c2、d2)为构建的二维平稳随机波阻抗模型.图 6(a3、b3、c3、d3)为主频40 Hz Ricker子波的合成地震记录.

图 6

Fig. 6

图 6 随机介质模型、合成地震记录和估计结果(从上到下为模型1—4) Fig. 6 Stationary r and om medium models，synthetic seismic records and estimated results

表 1(Table 1)

表 1 随机介质模型参数 Table 1 Parameters of r and om medium models 理论模型 理论参数 a/m b/ms θ/(°) 自相关函数 随机介质模型 合成地震记录 (主频40 Hz) 模型1 20 5 0 图 6(a1) 图 6(a2) 图 6(a3) 模型2 20 5 30 图 6(b1) 图 6(b2) 图 6(b3) 模型3 30 5 30 图 6(c1) 图 6(c2) 图 6(c3) 模型4 30 5 45 图 6(d1) 图 6(d2) 图 6(d3)

表 1 随机介质模型参数 Table 1 Parameters of r and om medium models

应用功率谱法对合成地震记录(图 6(a3、b3、c3、d3))分别进行了随机介质参数的估计，得到了如图 6(a4，b4，c4，d4)所示的二维自相关函数图像，估计出的特征参数如表 2所示.对比理论参数和估计结果，发现存在一定的误差，但还是比较接近理论值的，因此基于功率谱法的随机介质参数估计方法是可行的.

表 2(Table 2)

表 2 二维自相关函数参数估计结果对比 Table 2 Comparison of estimated results with different primary frequencies 模型 理论参数 估计结果 a/m b/ms θ/(°) a/m b/ms θ/(°) 二维自 相关函数 模型1 20 5 0 31.1 7 -1.4 图 6(a4) 模型2 20 5 30 21 7 25.2 图 6(b4) 模型3 30 5 30 33 7 26.9 图 6(c4) 模型4 30 5 45 24 6 44.3 图 6(d4)

表 2 二维自相关函数参数估计结果对比 Table 2 Comparison of estimated results with different primary frequencies

本试验的目的是测试不同主频的合成地震记录中随机介质参数估计的精度.为此构建了大小为Nx=Nt=300的高斯型随机介质波阻抗模型，其中dx=1 m，dt=1 ms，速度均值为5000 m·s^-1，速度标准差为500 m·s^-1，密度值根据Gardner公式计 算，随机介质参数为a=50 m，b=10 ms，θ=30°.图 7(a1)和(a2)分别为随机介质模型的自相关函数及其对应的随机介质波阻抗模型.根据波阻抗模型分别生成主频为20 Hz(图 7(b1))、30 Hz(图 7(c1))、40 Hz(图 7(d1))和50 Hz(图 7(e1))的合成地震记录，用于二维随机介质参数估计试验.

图 7

Fig. 7

图 7 不同主频地震数据对估计结果影响对比 Fig. 7 Comparison of estimated results with different primary frequencies

图 7(b2)、(c2)、(d2)、(e2)分别为计算得到的地震数据的自相关函数.图 7(b3)、(c3)、(d3)、(e3)为功率谱法估计得到的二维自相关函数，具体结果见表 3.以上结果与所设计的随机介质模型的二维自相关函数(图 7(a1))对比，其估计结果都比较接近理论值.但是，地震数据中子波的主频越高，估计的横向和纵向自相关长度相对真实值偏小的程度越大，Bean 等(1999)认为这是由于地震数据的频带宽度引起的.

表 3(Table 3)

表 3 不同主频地震数据的随机模型参数估计结果对比 Table 3 Comparison of estimated results of 2D r and om medium parameters 自相关函数 参数估计结果 a/m b/ms θ/(°) 理论模型 图 7(a2) 50 10 30 估计结果 20 Hz合成地震记录(图 7(b2)) 44 11 27 30 Hz合成地震记录(图 7(c2)) 41 8 26 40 Hz合成地震记录(图 7(d2)) 37 7 26.8 50 Hz合成地震记录(图 7(e2)) 33 7 26.7

表 3 不同主频地震数据的随机模型参数估计结果对比 Table 3 Comparison of estimated results of 2D r and om medium parameters

本试验的目的是测试在相同的随机介质参数下，地震数据的规模对于参数估计精度的影响.为此，我们构建了模型大小分别为Nx=Nt=200、300、400、500四种情况，而其他参数相同的高斯型随机介质波阻抗模型各10个.其中网格间距dx=1 m，dt=1 ms，速度均值为5000 m·s^-1，速度标准差为500 m·s^-1，a=50 m，b=20 ms，θ=30°.然后，对主频40 Hz Ricker子波合成地震记录进行参数估计，结果列在表 4—7中.网格数为Nx=Nt=200、300、400和500对应模型的估计结果，其综合平均误差分别为37.1%、19%、22.1%和22.1%.

表 4(Table 4)

表 4 随机介质参数估计结果(Nx=Nt=200) Table 4 Estimation of parameters of r and om media with the size: Nx=Nt=200 from Test 1 参数 理论参数 估计结果 1 2 3 4 5 6 7 8 9 10 误差/(%) a/m 50 37.2 39.2 37.7 37.7 33 37.8 28.2 19.7 40.8 40.5 29.5 b/ms 20 17.2 19.4 12.3 13.4 12.5 50.9 25.7 13.5 81.4 20.2 64.9 θ/(°) 30 27.5 26.7 21.4 25.2 21.8 31 26.0 26.8 25.6 19.4 16.7

表 4 随机介质参数估计结果(Nx=Nt=200) Table 4 Estimation of parameters of r and om media with the size: Nx=Nt=200 from Test 1

表 5(Table 5)

表 5 随机介质参数估计结果(Nx=Nt=300) Table 5 Estimation of parameters of r and om media with the size:Nx=Nt=300 from Test 2 参数 理论参数 估计结果 1 2 3 4 5 6 7 8 9 10 误差/(%) a/m 50 41.7 37.2 39 35.2 44.1 35.2 44.7 36.7 34.2 49.2 20.5 b/ms 20 14.4 17.5 13.9 19.2 15.2 11.8 18.4 16.2 17.6 14.7 20.2 θ/(°) 30 24.8 24.0 25.3 26.2 26.8 27.9 25.9 24.1 23 22.9 16.2

表 5 随机介质参数估计结果(Nx=Nt=300) Table 5 Estimation of parameters of r and om media with the size:Nx=Nt=300 from Test 2

表 6(Table 6)

表 6 随机介质参数估计结果(Nx=Nt=400) Table 6 Estimation of parameters of r and om media with the size: Nx=Nt=400 from Test 3 参数 理论参数 估计结果 1 2 3 4 5 6 7 8 9 10 误差/(%) a/m 50 42.2 40.2 44.7 37.2 38.7 44.7 40.7 34.2 36.0 45 19.2 b/ms 20 11.2 11.9 16.9 13.0 15.9 16.9 13.9 17.1 15.7 13 27.1 θ/(°) 30 26.3 25.3 22.2 22.2 22.9 22.2 26.4 23.8 23.4 24.6 20.1

表 6 随机介质参数估计结果(Nx=Nt=400) Table 6 Estimation of parameters of r and om media with the size: Nx=Nt=400 from Test 3

表 7(Table 7)

表 7 随机介质参数估计结果(Nx=Nt=500) Table 7 Estimation of parameters of r and om media with the size: Nx=Nt=500 from Test 4 参数 理论参数 估计结果 1 2 3 4 5 6 7 8 9 10 误差/(%) a/m 50 38.0 44.9 39.5 38.0 42 40.2 39.2 40.2 38.7 41.2 19.5 b/ms 20 11.9 12.7 14.5 11.9 12.8 11.1 11.4 12.5 13.1 19.3 34.1 θ/(°) 30 27.8 27.3 26.2 27.8 25.6 25.4 25.1 22.4 26.4 26.9 12.8

表 7 随机介质参数估计结果(Nx=Nt=500) Table 7 Estimation of parameters of r and om media with the size: Nx=Nt=500 from Test 4

总结上述三个试验的结果，可以发现自相关角度的估计误差较小，自相关长短轴的估计误差较大.总的来说，地震数据规模越大，随机介质参数估计误差越小.当地震数据的规模在200×200时的误差比较大，平均误差为37.1%；而当地震数据的规模在300×300以上，其估计平均误差均减小到20%左右.因此，假设a>b时，用于随机介质参数估计的地震数据的单边规模至少是a的5倍以上，可得到较为可信的结果.

我们认为，误差的主要来源可能是：(1)子波估计、功率谱估计中的加窗等多个步骤中都存在计算误差造成的；(2)在计算过程中，为了防止子波的某些频率的功率谱会出现零值，加白噪处理产生的误差；(3)更重要的因素，可能是随机介质地震记录不满足线性褶积模型的假设，从而产生计算误差. 5 实际地震数据的参数估计

图 8为歧口凹陷东一下亚段地震数据体中的一个典型剖面.从剖面可以看出，东一下亚段地层的底界内部发育有进积型地震相，以北西南东向分布于歧口滩海过渡区，地震相整体振幅为中等，连续性为中等或中、差连续，也有部分连续有顶超和底超终止现象，在终止点由于反射波复合干涉呈现弱振幅，一 般为弱-中强振幅.通过井震对比，我们将该进积型地震相解释为三角洲(图 8)，并划分为三角洲平原(1号区域)、三角洲前缘(2号区域)和前三角洲(3号区域)三个成因相.我们选取这三个典型区域的地震数据，进行随机介质参数的估计试验，估计结果如下：

图 8

Fig. 8

图 8 Line 2430典型地震相分析 Fig. 8 Analysis of typical seismic facies for the field seismic data in Line 2430

三角洲平原：自相关长度a为200.1 m；自相关长度b为10.3 ms；自相关角度θ约为-3.8°；

三角洲前缘：自相关长度a为348.5 m；自相关长度b为6.2 ms；自相关角度θ约为-1.7°；

前三角洲：自相关长度a为805.6 m；自相关长度b为15.8 ms；自相关角度θ约为-0.2°.

对比三角洲平原、三角洲前缘和前三角洲相，其介质的横向非均质性依次由强、中强过度到弱；而纵向非均质性，三角洲前缘的非均质性最强，三角洲平原次之，前三角洲最弱；而自相关角度的总体趋势为：其绝对值由三角洲平原向前三角洲逐渐变小，直至趋近于0.

从以上的分析结果可以看出，在三角洲沉积体系中，由于各成因相所对应的沉积环境不同，非均质体的分布特性必然有较大的差别.这种差别可以反映到随机介质的参数上，因此，随机介质参数估计结果可以用来识别地震相和沉积相. 6 结论和建议

随机介质参数估计是随机介质理论应用于地震勘探的关键.本文首先基于地震记录的褶积公式，推导出了二维情况下叠后地震数据与随机介质波阻抗模型的功率谱关系式，提出了功率谱法进行随机介质参数估计的方法，并结合实例详细介绍了估计的步骤.最后，通过多个理论模型和实际地震数据的估计试验，验证了方法是可行的和正确的.主要结论如下：

(1)随机介质参数估计的功率谱法，依据Wiener-Khintchine定理，可直接求得随机介质模型的自相关函数，不需要通过反演迭代计算过程，具有直观且高效率的优点；

(2)该方法是从二维叠后地震数据中估计平稳随机介质参数的方法，可同时估计出随机介质自相关函数在横、纵方向上的自相关长度以及自相关角度三个特征参数.

(3)通过地震主频与估计误差的分析表明，估计结果的偏差与地震主频具有一定的规律.对于同一个地震数据，由于各部分系统子波的主频差异较小，其估计结果的偏差也较为一致；

(4)为保证随机介质参数估计的精确度，本文提出了用于计算的地震数据规模的选取原则.一般情况下，地震数据单边的规模至少5倍于自相关长度以上，可得到20%左右误差的估计结果.如果用于多个差异较大(参数值的差异大于40%)的随机介质模型对比，估计结果还是较为可信的；

(5)本文针对实际地震数据中的三角洲平原、三角洲前缘和前三角洲三个成因相进行了随机介质的参数估计试验.结果表明，三角洲各成因相的随机介质参数具有较大的差异，进而可以为三角洲成因相的划分提供参考，反映了该方法拥有较好的应用前景.

本文仅对随机介质的参数估计方法进行了理论上的探索，模型试验的估计结果也还存在着较大的误差.我们认为，误差的主要来源可能是：子波估计、功率谱估计中的加窗和为了防止某些频率的子波功率谱出现零值进行的白噪处理，以及随机介质地震记录不完全满足线性褶积模型假设等原因造成的计算误差.

在应用于实际地震数据时，还需要在以下四个方面继续进行研究：

(1)寻求精度更高的子波估计、功率谱估计算法，以期从多个角度约束二维功率谱法随机介质参数估计的结果，减小估计误差；

(2)研究随机介质地震记录是否满足线性褶积模型的假设，改进或建立新的以波阻抗为基础的地震记录数学物理模型，提高功率谱法随机介质参数估计的精度；

(3)本文中的波阻抗模型是由速度与密度的经验公式计算得到的.由于公式中的参数是人为给定的，这显然与实际情况有一定的差距，因此是估计结果误差的来源之一.未来还将根据不同的沉积环境，对速度与密度的经验公式进行标定和校正，以减小估计误差；

(4)本文所提出的随机介质参数估计算法，主要是基于平稳介质的假设.实际介质的随机特性在空间上是变化的非平稳随机介质.因此，发展非平稳随机介质参数的估计方法将更接近于实际情况.