专利名称：一种控制全声波方程反演过程的地震勘探数据处理方法
技术领域：
本发明涉及一种控制地震数据全声波方程反演过程，提高地震勘探 全声波方程反演精度和对地质目标(尤其是岩性地层油气藏)的勘探能 力的地震勘探数据处理方法。
背景技术：
地震反演是利用野外观测的地震数据，通过数学方法和计算技术求 取地层物理参数，如速度、密度和阻抗等。全声波方程反演直接采用了
声学介质中地震波的传播方程，与波阻抗反演、弹性阻抗反演和AV0反 演相比，它能处理任意复杂的非均匀介质情况，提高地震反演的精度和 对地质目标(尤其是岩性地层油气藏)的勘探能力。
全声波方程反演计算费用昂贵，限制了它在实际生产中的应用。解 决这一问题的关键 一是提高声波方程计算波场的效率，另一是提高反 演的收敛速度。Tarantola等(1988)提出的深度域延拓计算地震波场和 Pratt等(1999)提出的通过计算部分频率的波场值实现全波场的模拟都 显著提高了求解波动方程的速度。这两种方法减少了单次迭代反演所需 的时间，但由于地震反演基本上都是迭代进行的，迭代次数越多，反演 的时间越长，费用就越高。
对于波动方程这类巨大计算量的反演，目前采用的最好的算法是最 速下降法或共轭梯度法。这两种算法收敛速度比较快，且比较稳定，但仍需要一定的迭代次数。要将波动方程反演用于实际地震数据处理，常 规的最速下降法或共轭梯度法所需的迭代次数是不能被允许的，也就是 说迭代次数是目前波动方程反演工业化应用的一个瓶颈。
现有最速下降法或共轭梯度法中，人们采用线性搜索法或抛物线拟 合法可以在一定程度上降低迭代次数，但在每一迭代中至少增加了两次 额外的计算地震波场工作，显著增加了每一步的计算成本。
发明内容
本发明的目的是利用一种控制方法对野外实测地震数据进行声波方 程反演过程的自动控制，降低波动方程反演的计算成本，提高地震反演 的精度和对地层油气藏的判断能力。
本发明是一种基于声波理论对野外实测地震数据进行反演，并通过 对反演过程的控制，提高反演的速度和精度。具体步骤包括
(1) .野外激发、接收获得地震波数据，用常规处理方法对该数据进行 静校正、地表一致性振幅补偿和叠前去除噪音，获得炮集数据；
(2) .从步骤(1)获得的炮集数据中抽取共中心点道集数据，进行速度 分析，时间域层速度计算和时间-深度转换，获得深度域体积模量和密 度初始模型；
① 步骤(2)中速度分析为常规地震速度分析；
② 步骤(2)中时间域层速度用速度分析获得的均方根速度和Dix公式进
行计算，2di
V;=-
/ 一/ (1) 。
。-1
其中ti为第i层双程旅行时，Vi为第i层的速度；
③ 步骤(2)中时间-深度转换是按下式将时间域的层速度转换成深度域层
速度，
=会(^一,/—》/ (2)
其中ti为第i层双程旅行时，hi为第i层的厚度，Vi为第i层的速度;
④ 步骤(2)中密度(模型参数l)初始模型用Gardner公式计算，
P = 0.31/25 (3)
其中p为密度，v为纵波速度。
步骤(2)中体积模量(模型参数2)初始模型用下列方法计算
《="2 (4) 其中K为体积模量，p为密度，v为纵波速度。
(3).利用步骤2建立的初始体积模量和密度模型，进行地震波场正演
模拟；
① 步骤(3)中地震波场正演模拟利用的是声学介质中的波动方程；
② 步骤(3)中用常规的伪谱法进行波场模拟，即对空间变量进行付氏变
换，在付氏域中计算声波方程中波场对空间变量的偏导数；声波方 程中波场对时间变量的偏导数用二阶中心差分法计算；
③ 步骤(3)中震源函数采用零相位雷克子波；(4) .按下式计算步骤(3)模拟的炮集数据与步骤1获得测量数据的差 (残差数据)-
A""幽-《/ (5)
其中Ad为残差数据，d。b,为测量地震数据，d^为计算地震数据。
(5) .按下式计算误差能量，
E二丄AdHd 2
D—一 (6)
式中E为误差能量，Ad'为Ad的转置，CD为数据协方差矩阵； 当E〈別寸，输出反演结果，并停止反演，否则，继续步骤6 步骤9， 直到E〈S， 5为任意给的一个非常小的数， 一般取弘1.0e-4 1.0e-6;
步骤(5)中数据协方差矩阵用下列方法计算
2
C n =
1 ^
"X (7)
其中X为到测量点的距离；t为地震波到达的时间；(Td为数据方差； 指数p，对于二维问题，对于三维情况，取;^1;
(6) .将步骤(4)获得的残差数据作为源数据，用正演模拟类似的方法 进行逆时传播，获得残差数据的剩余地震波场；
① 步骤(6)中的逆时传播与步骤(3)的波场正演过程相反，即解波动
方程，按反时间方向进行波场传播；
② 与骤(3)正演模拟所用的震源相对应，步骤(6)中进行波场模拟时，
震源用的是步骤(4)计算的残差数据；
(7) .将步骤3获得的正演地震波场和步骤(6)获得剩余地震波场分别对时间变量进行求导，进行零延迟互相关，然后对时间变量进行积分，
对激发炮数进行累加，获得初始密度参数的共轭修改量；将步骤(3) 获得的正演地震波场和步骤(6)获得剩余地震波场分别对空间变量 求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对 激发炮数进行累加，获得初始体积参数的共轭修改量，即按下式计算 模型参数的共轭修改量；
x (x) s (8)
<5/5 (x) = ~^——5] [y伊"^ ^ (x,f;x丄grac 尸(x,Z;xJ P (x) " ■*>
其中P为正演波场，M/为回传的剩余波场，户或^表示波场对时间的一 阶偏导数，grad( *)为梯度运算，T为地震记录长度，K和p为模型参数， c^和^为模型参数共轭量；
(8).利用步骤(7)获得的模型参数的共轭修改量，按下列方法的迭代 修改模型参数，
m +1=m - g (9) 其中，m=(K， p)为模型，mnw和mn分别为第n次修改后的模型和当前 模型；an为第n次迭代修改步长；gn为第n次迭代的梯度； 或①步骤(8)中梯度gn用下列方法获得，
g" = (^"，^P") (10)
其中^^和c5^是第n步迭代的模型参数的共轭修改量，用骤(7)公
式计算。② 步骤(8)中修改步长Otn为反演过程的控制，用下列方法计算
= [Ad'C》(m" A: E (11) 其中，p(m)为映射函数，yt为反馈增益系数，J为地震数据对模型参数 的导数，E为误差能量；
③ 步骤(8)②中映射函数用下列方法计算
当n〉1时，函数p(m") = g — g —!，
当n=l时，—Ju 步骤(8)②中反馈增益系数取满足ytX10A,)的任一数(ts为系统固有
时间，ts取5 20范围任一值)； (9).将修改后的模型数据作为新的初始模型，重复步骤(3) 步骤(5)。 发明的效果
本发明通过对反演过程进行自动控制，使声波方程反演稳定、快速 收敛。理论模型测试结果表明控制反演收敛速度比常规方法显著提高了 近5倍。实例应用为我国西部某气田实际数据反演，经过两次迭代后， 获得的反演结果与实测数据基本一致。
-
图1是可控声波方程反演流程图2是理论模型合成的地震数据；
图3是反演结果与真实数据图4是反演误差能量随迭代次数变化曲线；
图5是野外测量地震数据图(左为频谱，右为炮集数据)；
图6是经过两次迭代反演的结果与实测数据图；实施例l:
理论模型为一个三层介质模型，图2中黑色曲线为三层介质模型参 数随深度变化曲线。图3为理论模型合成地震数据，对该数据进行反演， 具体实现步骤为-
1. 采集理论模型地震数据，抽取共中心点道集数据，进行速度分析，时 间域层速度计算和时间-深度转换，获得深度域体积模量和密度初始模 型；
2. 基于初始模型，用伪谱法模拟声学介质中地震波场；
3. 按式(5)计算理论模型与初始模型模拟数据的差；
4. 给定S-1.0e-6，按式(6)计算误差能量E，当E《J，停止反演，输出 反演结果，如果E〉5，继续下列步骤；
5. 步骤(3)计算的残差数据作为输入数据，用伪谱法计算声波方程的波 场，并按反时间方向传播；
6. 按式(8)计算体积模量和密度模型的共轭修改量；
7. 按式(11)计算修改步长，并用(10)式获得梯度；
8. 用(9)式对初始模型进行修改；
9. 将修改后的模型数据作为新的初始模型，重复步骤2 步骤4。
图2灰线显示了 5次迭代后的反演结果，与真实模型数据吻合的很 好。图4显示了分别用常规共轭梯度法和可控反演方法反演的误差能量 随迭代次数变化曲线，控制反演方法收敛速度快，约为常规方法的5倍， 且稳定。
实施例2:实例2是位于我国西部的某气田。该气田为陆相湖泊沉积，水下分 流河道发育。目标层埋藏深度大、埋藏时间长、成岩演化程度高、物性 差，属于典型的低孔隙度、低渗透率气田。储层纵向上厚度大、平面上 大面积分布，但有效砂岩单层厚度薄、横向变化大，非均质性强。沉积 微相、物性分析和试气成果证实，有效储层主要发育在三角洲平原和三 角洲前缘亚相中的高能分流河道中。
由于有效储层和围岩纵波速度和阻抗差异小，常规地震反演结果达
不到勘探的要求。2006年我们用可控声波方程叠前反演方法对该区地震 数据进行了反演。图5为野外测量地震数据(左为其频谱，右为炮集数 据)，用于地震反演的输入，图6显示了两次迭代反演后钻井处的反演结 果(灰线)和实际测井数据(黑线)，反演获得数据与实测数据基本一致。 具体实现步骤为
1. 采集地震数据，用常规方法对地震进行静校正、地表一致性振幅补偿 和叠前去噪，形成炮集数据；
2. 用常规方法抽取共中心点道集数据，进行速度分析，时间域层速度计 算和时间-深度转换，获得深度域体积模量和密度初始模型；
3. 10.同"实施例l"中步骤2 9。
1权利要求
1.一种控制全声波方程反演过程的地震勘探数据处理方法，其特征在于具体步骤包括(1).野外激发、接收获得地震波数据，用常规处理方法对该数据进行静校正、地表一致性振幅补偿和叠前去除噪音，获得炮集数据；(2).从步骤(1)获得的炮集数据中抽取共中心点道集数据，进行速度分析，时间域层速度计算和时间-深度转换，获得深度域体积模量和密度初始模型；①步骤(2)中速度分析为常规地震速度分析；②步骤(2)中时间域层速度用速度分析获得的均方根速度和Dix公式进行计算，<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><msubsup> <mi>v</mi> <mi>i</mi> <mn>2</mn></msubsup><mo>=</mo><mfrac> <mrow><msubsup> <mi>v</mi> <mrow><mi>r</mi><mo>,</mo><mi>i</mi> </mrow> <mn>2</mn></msubsup><msub> <mi>t</mi> <mi>i</mi></msub><mo>-</mo><msubsup> <mi>v</mi> <mrow><mi>r</mi><mo>,</mo><mi>i</mi><mo>-</mo><mn>1</mn> </mrow> <mn>2</mn></msubsup><msub> <mi>t</mi> <mrow><mi>i</mi><mo>-</mo><mn>1</mn> </mrow></msub> </mrow> <mrow><msub> <mi>t</mi> <mi>i</mi></msub><mo>-</mo><msub> <mi>t</mi> <mrow><mi>i</mi><mo>-</mo><mn>1</mn> </mrow></msub> </mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo></mrow> </mrow>]]></math></maths>其中ti为第i层双程旅行时，vi为第i层的速度；③步骤(2)中时间-深度转换是按下式将时间域的层速度转换成深度域层速度，<maths id="math0002" num="0002" ><math><![CDATA[ <mrow><msub> <mi>h</mi> <mi>i</mi></msub><mo>=</mo><mfrac> <mn>1</mn> <mn>2</mn></mfrac><mrow> <mo>(</mo> <msub><mi>t</mi><mi>i</mi> </msub> <mo>-</mo> <msub><mi>t</mi><mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn></mrow> </msub> <mo>)</mo></mrow><msub> <mi>v</mi> <mi>i</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow>]]></math></maths>其中ti为第i层双程旅行时，hi为第i层的厚度，vi为第i层的速度；④步骤(2)中密度初始模型用Gardner公式计算，ρ＝0.31v0.25 (3)其中ρ为密度，v为纵波速度。⑤步骤(2)中体积模量初始模型用下列方法计算K＝ρv2(4)其中K为体积模量，ρ为密度，v为纵波速度。(3).利用步骤2建立的初始体积模量和密度模型，进行地震波场正演模拟；①步骤(3)中地震波场正演模拟利用的是声学介质中的波动方程；②步骤(3)中用常规的伪谱法进行波场模拟，即对空间变量进行付氏变换，在付氏域中计算声波方程中波场对空间变量的偏导数；声波方程中波场对时间变量的偏导数用二阶中心差分法计算；③步骤(3)中震源函数采用零相位雷克子波；(4).按下式计算步骤(3)模拟的炮集数据与步骤1获得测量数据的差(残差数据)Δd＝dobs-dcal(5)其中Δd为残差数据，dobs为测量地震数据，dcal为计算地震数据。(5).按下式计算误差能量，<maths id="math0003" num="0003" ><math><![CDATA[ <mrow><mi>E</mi><mo>=</mo><mfrac> <mn>1</mn> <mn>2</mn></mfrac><msup> <mi>Δd</mi> <mi>t</mi></msup><msubsup> <mi>C</mi> <mi>D</mi> <mrow><mo>-</mo><mn>1</mn> </mrow></msubsup><mi>Δd</mi><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo></mrow> </mrow>]]></math></maths>式中E为误差能量，Δdt为Δd的转置，CD为数据协方差矩阵；当E＜δ时，输出反演结果，并停止反演，否则，继续步骤6～步骤9，直到E＜δ。δ为任意给的一个非常小的数，一般取δ＝1.0e-4～1.0e-6；步骤(5)中数据协方差矩阵用下列方法计算<maths id="math0004" num="0004" ><math><![CDATA[ <mrow><msub> <mi>C</mi> <mi>D</mi></msub><mo>=</mo><mfrac> <mn>1</mn> <mi>X</mi></mfrac><mo>·</mo><mfrac> <msubsup><mi>σ</mi><mi>d</mi><mn>2</mn> </msubsup> <msup><mi>t</mi><mrow> <mn>2</mn> <mi>p</mi></mrow> </msup></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo></mrow> </mrow>]]></math></maths>其中X为到测量点的距离；t为地震波到达的时间；σd为数据方差；指数p，对于二维问题，取p≥0.5，对于三维情况，取p≥1；(6).将步骤(4)获得的残差数据作为源数据，用正演模拟类似的方法进行逆时传播，获得残差数据的剩余地震波场；①步骤(6)中的逆时传播与步骤(3)的波场正演过程相反，即解波动方程，按反时间方向进行波场传播；②与骤(3)正演模拟所用的震源相对应，步骤(6)中进行波场模拟时，震源用的是步骤(4)计算的残差数据；(7).将步骤3获得的正演地震波场和步骤(6)获得剩余地震波场分别对时间变量进行求导，进行零延迟互相关，然后对时间变量进行积分，对激发炮数进行累加，获得初始密度参数的共轭修改量；将步骤(3)获得的正演地震波场和步骤(6)获得剩余地震波场分别对空间变量求导，并进行零延迟互相关，将互相关结果对时间变量进行积分，对激发炮数进行累加，获得初始体积参数的共轭修改量，即按下式计算模型参数的共轭修改量；<maths id="math0005" num="0005" ><math><![CDATA[ <mrow><mi>δ</mi><mover> <mi>K</mi> <mo>^</mo></mover><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <mn>1</mn> <mrow><msup> <mi>K</mi> <mn>2</mn></msup><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow> </mrow></mfrac><munder> <mi>Σ</mi> <mi>s</mi></munder><msubsup> <mo>∫</mo> <mn>0</mn> <mi>T</mi></msubsup><mi>dt</mi><mover> <mi>P</mi> <mo>·</mo></mover><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>t</mi> <mo>;</mo> <msub><mi>x</mi><mi>s</mi> </msub> <mo>)</mo></mrow><mover> <mi>ψ</mi> <mo>·</mo></mover><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>t</mi> <mo>;</mo> <msub><mi>x</mi><mi>s</mi> </msub> <mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo></mrow> </mrow>]]></math></maths><maths id="math0006" num="0006" ><math><![CDATA[ <mrow><mi>δ</mi><mover> <mi>ρ</mi> <mo>^</mo></mover><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <mn>1</mn> <mrow><msup> <mi>ρ</mi> <mn>2</mn></msup><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow> </mrow></mfrac><munder> <mi>Σ</mi> <mi>s</mi></munder><msubsup> <mo>∫</mo> <mn>0</mn> <mi>T</mi></msubsup><mi>dt</mi><mi>grad</mi><mi>ψ</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>t</mi> <mo>;</mo> <msub><mi>x</mi><mi>s</mi> </msub> <mo>)</mo></mrow><mo>·</mo><mi>grad</mi><mi>P</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>t</mi> <mo>;</mo> <msub><mi>x</mi><mi>s</mi> </msub> <mo>)</mo></mrow> </mrow>]]></math></maths>其中P为正演波场，ψ为回传的剩余波场， id="icf0007" file="A2008101167140004C3.tif" wi="2" he="3" top= "238" left = "121" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>或 id="icf0008" file="A2008101167140004C4.tif" wi="2" he="3" top= "239" left = "131" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>表示波场对时间的一阶偏导数，grad(·)为梯度运算，T为地震记录长度，K和ρ为模型参数， id="icf0009" file="A2008101167140004C5.tif" wi="5" he="4" top= "260" left = "27" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>和 id="icf0010" file="A2008101167140004C6.tif" wi="4" he="4" top= "261" left = "39" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>为模型参数共轭量；(8).利用步骤(7)获得的模型参数的共轭修改量，按下列方法的迭代修改模型参数，mn+1＝mn-αngn (9)其中，m＝(K，ρ)为模型，mn+1和mn分别为第n次修改后的模型和当前模型；αn为第n次迭代修改步长；gn为第n次迭代的梯度；或①步骤(8)中梯度gn用下列方法获得，<maths id="math0007" num="0007" ><math><![CDATA[ <mrow><msub> <mi>g</mi> <mi>n</mi></msub><mo>=</mo><mrow> <mo>(</mo> <mi>δ</mi> <msub><mover> <mi>K</mi> <mo>^</mo></mover><mi>n</mi> </msub> <mo>,</mo> <mi>δ</mi> <msub><mover> <mi>ρ</mi> <mo>^</mo></mover><mi>n</mi> </msub> <mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo></mrow> </mrow>]]></math></maths>其中 id="icf0012" file="A2008101167140005C2.tif" wi="6" he="5" top= "108" left = "44" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>和 id="icf0013" file="A2008101167140005C3.tif" wi="5" he="4" top= "109" left = "57" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>是第n步迭代的模型参数的共轭修改量，用骤(7)公式计算。②步骤(8)中修改步长αn为反演过程的控制，用下列方法计算其中， id="icf0015" file="A2008101167140005C5.tif" wi="9" he="4" top= "153" left = "43" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>为映射函数，k为反馈增益系数，J为地震数据对模型参数的导数，E为误差能量；③步骤(8)②中映射函数用下列方法计算当n＞1时，函数 id="icf0016" file="A2008101167140005C6.tif" wi="38" he="4" top= "186" left = "71" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>当n＝1时， id="icf0017" file="A2008101167140005C7.tif" wi="20" he="4" top= "197" left = "61" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/>④步骤(8)②中反馈增益系数取满足k＞(10/ts)的任一数，ts为系统固有时间，ts取5～20范围任一值；(9).将修改后的模型数据作为新的初始模型，重复步骤(3)～步骤(5)。
全文摘要
本发明涉及一种控制全声波方程反演过程的地震勘探数据处理方法，采集理论模型地震数据，获得深度域体积模量和密度初始模型；用伪谱法模拟声学介质中地震波场；计算理论模型与初始模型模拟数据的差；给定δ＝1.0e-6，计算误差能量E，当E≤δ，停止反演，输出反演结果，如果E＞δ，继续下列步骤；计算残差数据；计算体积模量和密度模型的共轭修改量；计算修改步长，获得梯度；对初始模型进行修改；将修改后的模型数据作为新的初始模型，通过对反演过程进行自动控制，使声波方程反演稳定、快速收敛，反演收敛速度比常规方法提高了近5倍。
文档编号G01V1/28GK101630018SQ20081011671
公开日2010年1月20日 申请日期2008年7月16日 优先权日2008年7月16日
发明者石玉梅 申请人:中国石油天然气股份有限公司