CN109541681A - 一种拖缆地震数据和少量obs数据联合的波形反演方法 - Google Patents
一种拖缆地震数据和少量obs数据联合的波形反演方法 Download PDFInfo
- Publication number
- CN109541681A CN109541681A CN201710864481.4A CN201710864481A CN109541681A CN 109541681 A CN109541681 A CN 109541681A CN 201710864481 A CN201710864481 A CN 201710864481A CN 109541681 A CN109541681 A CN 109541681A
- Authority
- CN
- China
- Prior art keywords
- obs
- data
- streamer
- wave field
- model
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 37
- 238000007619 statistical method Methods 0.000 claims abstract description 3
- 238000004364 calculation method Methods 0.000 claims description 30
- 238000001914 filtration Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 230000017105 transposition Effects 0.000 claims description 3
- 238000007689 inspection Methods 0.000 claims description 2
- 238000004088 simulation Methods 0.000 description 7
- 238000005070 sampling Methods 0.000 description 4
- 238000004587 chromatography analysis Methods 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 2
- 239000000284 extract Substances 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 238000003786 synthesis reaction Methods 0.000 description 2
- 241001269238 Data Species 0.000 description 1
- 240000007762 Ficus drupacea Species 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 244000145845 chattering Species 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004513 sizing Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/61—Analysis by combining or comparing a seismic data set with other data
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及一种拖缆地震数据和少量OBS数据联合的波形反演方法,属于海洋地震勘探的全波形反演速度建模领域。本发明主要包括如下步骤:1)获取观测的拖缆和OBS地震记录数据,并对其进行预处理;2)采用高阶统计方法提取地震子波;3)分别计算拖缆和OBS地震记录对应的正传波场和反传波场;4)基于伴随状态法计算目标函数的梯度场;5)计算共轭梯度场和速度模型的修正量;6)迭代更新速度模型;7)判断是否满足迭代终止条件,如果不满足迭代终止条件,则返回步骤3,否则,输出结果。本方法解决了只利用拖缆地震记录进行波形反演,由于缺少低频地震信息,反演结果不准确的问题。
Description
技术领域
本发明涉及一种拖缆地震数据和少量OBS数据联合的波形反演方法,属于地震勘探的全波形反演速度建模领域。
背景技术
目前速度建模的方法有很多种,比如走时层析、斜率层析、菲尼尔体层析、偏移速度分析、全波形反演等。在这些方法中,全波形反演是近年来的研究热点之一,也是最有可能得到高分辨率、高保真度的速度建模方法之一。
然而,全波形反演严重依赖于初始速度模型的准确程度和观测地震记录中的低频地震信息。对于地质结构复杂的区域,利用现有的手段为全波形反演建立一个足够准确的初始速度模型通常比较困难。因此,为全波形反演提供足够低频的地震数据就显得尤为重要。但是,由于常用的采集技术的限制,由水听器记录的拖缆地震数据中通常缺少低频成份(< 5 Hz)。利用拖缆地震数据进行全波形反演,其反演通常收敛于局部解,反演结果不可信。
为此,许多学者提出了多种方法来构建低频地震信息。一种方法是利用高频地震资料通过非线性插值计算低频地震资料。该方法虽然在一定程度上可以抑制周期跳跃现象的发生,然而,插值得到的低频地震资料不满足波动方程。Beat tone 全波形反演利用数学变换从频率相近的两个高频信号中得到低频地震资料,然而构建的低频地震资料的质量严重依赖于高频地震数据。拉普拉斯域全波形反演在缺少低频地震资料时,可以得到足够平滑的速度模型,然而该方法需要大偏移距的地震数据。包络反演也可以构建低频信息,然而数学变换严重地破坏了地震波的波形。
这些方法主要通过数学变换来得到低频的地震信息,从而满足全波形反演的需要。然而,这些人工合成的低频地震记录与真实的低频地震信息存在一定程度的差异,使得反演结果中可能存在虚假的速度结构,降低了反演结果的可靠性。
发明内容
针对现有技术存在的上述缺陷,本发明提出了一种拖缆地震数据和少量OBS数据联合的波形反演方法,利用少量OBS观测获取实际的地震低频分量,解决了只使用拖缆数据进行全波形反演时,因缺乏低频分量使反演结果不可靠的问题。
本发明是采用以下的技术方案实现的:本发明所述一种拖缆地震数据和少量OBS数据联合的波形反演方法,包括如下步骤:
1)构建拖缆地震数据和少量OBS数据联合的混合域全波形反演目标函数,构建的目标函数如下:
, (1)
式中,E是目标函数值,α是权重系数,nr是每一炮水听器的个数,s是炮号,r是水听器号,t是时间,xs是模型空间中炮点s的位置,no是每一炮OBS的个数,o指OBS号,dcal(r,t|xs)是预测的拖缆地震记录,drea(r,t|xs)是观测的拖缆地震记录,dcal(o,t|xs)是预测的OBS地震记录,drea(o,t|xs)是观测的OBS地震记录;本发明通过最小化该目标函数来得到最优的反演结果;
2)获取野外观测的拖缆和OBS地震记录;
3)对观测数据进行预处理;通过带通滤波、F-K滤波等数据处理手段从观测数据中去除声波方程无法正演模拟的面波、转换波、随机噪声等;
4)采用高阶统计方法,从观测资料中提取地震子波;
5)计算权重系数;在反演迭代的初期,主要恢复模型的背景速度场,所以OBS资料应该起主导作用,而在反演后期,主要提高模型的分辨率,所以拖缆数据应该占主要地位;因此,本发明采用下式计算权重系数:
, (2)
式中,k是当前的迭代次数,nite是反演总共的迭代次数;
6)计算预测的拖缆地震记录及其对应的正传波场;利用给定的初始速度模型,计算正传波场,正传波场的计算公式如下:
, (3)
式中,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,v(x)是模型空间x点处的速度,Pfr(x,t|xs)是拖缆地震记录对应的时间域正传波场,sfr(x,t|xs)是从预处理后观测的拖缆地震记录中提取的子波;预测的拖缆地震记录的计算公式如下:
, (4)
式中,r表示水听器号,t表示时间,xs表示模型空间中炮点位置,x表示模型空间位置,xr是水听器在模型空间的位置,Pfr(x,t|xs)表示时间域正传波场,R表示将检波点处的波场值提取出来作为该时刻检波点的地震记录,dcal(r,t|xs)表示预测的拖缆地震记录;
7)计算预测的OBS地震记录及其对应的正传波场;利用给定的初始速度模型,计算正传波场,正传波场的计算公式如下:
, (5)
式中,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,v(x)是模型空间x点处的速度,Pfo(x,t|xs)是OBS地震记录对应的时间域正传波场,sfo(x,t|xs)是从预处理后观测的OBS地震记录中提取的子波;预测的OBS地震记录的计算公式如下:
, (6)
式中,o表示OBS号,t表示时间,xs表示模型空间中炮点位置,x表示模型空间位置,xo是OBS在模型空间的位置,Pfo(x,t|xs)表示时间域正传波场,R表示将OBS点处的波场值提取出来作为该时刻OBS的地震记录,dcal(o,t|xs)表示预测的OBS地震记录;
8)将预测地震记录和观测地震记录带入公式1,计算目标函数值;
9)计算拖缆地震记录的残差和拖缆地震记录对应的反传波场;拖缆地震记录残差的计算公式如下:
, (7)
式中,s是炮号,r是水听器号,t是时间,xs表示模型空间中炮点位置;拖缆地震记录对应的反传波场的计算公式如下:
, (8)
式中,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,v(x)是模型空间x点处的速度,Pbr(x,t|xs)是拖缆地震记录对应的时间域反传波场,sbr(x,t|xr)是第s炮拖缆地震记录的残差;
10)计算OBS地震记录的残差和OBS地震记录对应的反传波场;OBS地震记录残差的计算公式如下:
, (9)
式中,s是炮号,o是OBS号,t是时间,xs表示模型空间中炮点位置;OBS地震记录对应的反传波场的计算公式如下:
, (10)
式中,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,v(x)是模型空间x点处的速度,Pbo(x,t|xs)是OBS地震记录对应的时间域反传波场,sbo(x,t|xr)是第s炮OBS地震记录的残差;
11)将时间域拖缆和OBS地震记录对应的正传和反传波场变换到频率域,计算公式如下:
, (11)
式中,T是采样时间长度,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,ω表示角频率;
12)计算目标函数的梯度;梯度的计算公式如下:
, (12)
式中,α 是权重系数,nr是水听器的个数,no是OBS的个数,ωr是使用拖缆数据时选择的角频率,ωo是使用OBS数据时选择的角频率,s是炮号,x是模型空间的位置,xs是炮点在模型空间的位置,Pfr(x,ωr|xs)是拖缆数据对应的正传波场,Pfo(x,ωo|xs)是OBS数据对应的正传波场,Pbr(x,ωr|xs)是拖缆数据对应的反传波场,Pbo(x,ωo|xs)是OBS数据对应的反传波场,上标*表示复共轭转置;
13)计算共轭梯度场;共轭梯度场的计算公式如下:
, (13)
式中,dk-1是第k-1次的共轭梯度,βk是使相邻两次共轭梯度正交的系数,其计算公式如下:
, (14)
其中,gk是梯度向量,dk-1是共轭梯度向量,符号<,>表示内积运算;
14)迭代更新速度模型;速度模型更新的计算公式如下:
, (15)
式中,vk+1是第k次迭代更新后的速度,λ是迭代更新的步长,可以通过线性搜索的方法得到;
15)判断是否满足迭代终止条件;如果满足终止条件,则输出计算结果,否则,将更新后的速度模型作为新的初始速度模型,返回步骤5)。
进一步地,所述步骤4)中,提取的地震子波有两个,分别是从观测的拖缆地震记录提取的子波sfr和从观测的OBS地震记录中提取的地震子波sfo。
进一步地,公式3、5、8和10都采用时间2阶、空间12阶规则网格有限差分方法求解,同时采用完全匹配层吸收边界条件来压制模型边界处的反射波。
本发明的有益效果是:海洋拖缆地震数据增加少量OBS地震数据的联合地震波形反演,在反演初期,主要利用少量OBS地震数据来反演模型的背景速度场,而在反演后期,主要利用拖缆地震数据来反演模型的高频成份;克服了只利用拖缆数据进行反演,由于缺少低频地震信息,导致反演结果不准的问题;其次,相对于现有方法所利用的人工计算的低频地震记录, OBS观测的地震记录的低频信息,更真实、更可靠,采用OBS地震记录进行反演,反演结果的可靠性高。
附图说明
图1为本发明的流程图。
图2为本发明建立的真实速度模型图。
图3为本发明所用的初始速度模型图。
图4为第17炮拖缆地震记录图。
图5为第17炮OBS地震记录图。
图6为常规全波形反演方法只使用拖缆数据的反演结果图。
图7为常规全波形反演方法只使用OBS地震记录的反演结果图。
图8为本发明的反演结果图。
具体实施方式
为了使本发明的目的、技术方案更加清楚明白,下面以理论模型为例并结合附图,对本发明作进一步详细说明。
本发明的流程图,如图1所示,主要包括如下步骤:
1)获取观测地震记录;给定一个真实速度模型,采用高频的雷克子波作为震源,利用常密度声波方程正演模拟,将得到的地震记录作为观测的拖缆地震记录;采用低频的雷克子波作为震源,利用常密度声波方程正演模拟,将得到的地震记录作为观测的OBS地震记录;正演模拟的计算公式如下:
, (1)
式中,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,v(x)是模型空间x点处的速度,P(x,t|xs)表示时间域波场,s(x,t|xs)表示震源函数;
2)计算权重系数;在反演迭代的初期,主要恢复模型的背景速度场,所以OBS资料应该起主导作用,而在反演后期,主要提高模型的分辨率,所以拖缆数据应该占主要地位;本发明计算权重系数的公式如下:
, (2)
式中,k是当前的迭代次数,nite是反演总共的迭代次数;
3)计算正传波场;给定初始速度模型和高频雷克子波sfr,利用声波方程正演模拟得到拖缆地震记录对应的时间域正传波场Pfr;正传波场Pfr的计算公式如下:
, (3)
式中,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,v(x)是模型空间x点处的速度,sfr(x,t|xs)是高频雷克子波;
利用给定的初始速度模型和低频雷克子波sfo,利用声波方程正演模拟得到OBS地震记录对应的正传波场Pfo;正传波场Pfo的计算公式如下:
, (4)
式中,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,v(x)是模型空间x点处的速度,sfo(x,t|xs)是低频的雷克子波;
4)计算预测的地震记录;预测的拖缆地震记录的计算公式如下:
, (5)
式中,r表示水听器号,t表示时间,xs表示模型空间中炮点位置,x表示模型空间位置,xr是水听器在模型空间的位置,Pfr(x,t|xs)表示拖缆地震记录对应的时间域正传波场,dcal(r,t|xs)表示预测的拖缆地震记录,R(x|xr,xs)表示数据选择函数,当x等于xr时,R等于1,否则为0;
预测的OBS地震记录的计算公式如下:
, (6)
式中,o表示OBS号,t表示时间,xs表示模型空间中炮点位置,x表示模型空间位置,xo是OBS在模型空间的位置,Pfo(x,t|xs)表示OBS地震记录对应的时间域正传波场,dcal(o,t|xs)表示预测的OBS地震记录,R(x|xo,xs)表示数据选择函数当x等于xo时,R等于1,否则为0;
5)计算地震记录残差;将预测的拖缆地震记录和观测的拖缆地震记录作差,得到拖缆地震记录的残差,计算公式如下:
, (7)
式中,s是炮号,r是水听器号,t是时间,xs表示模型空间中炮点位置,dcal(r,t|xs)是预测的拖缆地震记录,drea(r,t|xs)是观测的拖缆地震记录,derr(r,t|xs)是拖缆地震记录的残差;
将预测的OBS地震记录和观测的OBS地震记录作差,得到OBS地震记录的残差,计算公式如下:
, (8)
式中,s是炮号,o是OBS号,t是时间,xs表示模型空间中炮点位置,dcal(o,t|xs)是预测的OBS地震记录,drea(o,t|xs)是观测的OBS地震记录,derr(o,t|xs)是OBS地震记录的残差;
6)计算目标函数值,目标函数值的计算公式如下:
, (9)
式中,E是目标函数值,α是权重系数,nr是每一炮水听器的个数,s是炮号,r是水听器号,t是时间,xs是模型空间中炮点s的位置,no是每一炮OBS的个数,o指OBS号,dcal(r,t|xs)是预测的拖缆地震记录,drea(r,t|xs)是观测的拖缆地震记录,dcal(o,t|xs)是预测的OBS地震记录,drea(o,t|xs)是观测的OBS地震记录;
7)计算反传波场;将地震记录残差作为震源函数,利用给定的初始速度模型,计算反传波场;拖缆地震记录对应的反传波场Pbr的计算公式如下:
, (10)
式中,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,v(x)是模型空间x点处的速度,sbr(x,t|xr)是第s炮拖缆地震记录的残差;
OBS地震记录对应的反传波场Pbo的计算公式如下:
, (11)
式中,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,v(x)是模型空间x点处的速度,sbo(x,t|xr)是第s炮OBS地震记录的残差;
8)将时间域拖缆和OBS数据对应的正传、反传波场变换到频率域,其计算公式如下:
, (12)
式中,T是采样时间长度,x表示模型空间位置,t表示时间,xs表示模型空间中炮点位置,ω表示角频率,P(x,t|xs)表示时间域波场,P(x,ω|xs)表示频率域波场;
9)计算目标函数的梯度;利用步骤3)计算出的正传波场和步骤7)计算出的反传波场,来计算目标函数的梯度场,计算公式如下:
, (13)
式中,α是权重系数,nr是水听器的个数,no是OBS的个数,ωr是使用拖缆数据时选择的角频率,ωo是使用OBS数据时选择的角频率,s是炮号,x是模型空间的位置,xs是炮点在模型空间的位置,Pfr(x,ωr|xs)是拖缆数据对应的正传波场,Pfo(x,ωo|xs)是OBS数据对应的正传波场,Pbr(x,ωr|xs)是拖缆数据对应的反传波场,Pbo(x,ωo|xs)是OBS数据对应的反传波场,上标*表示复共轭转置;
10)计算目标函数的共轭梯度;共轭梯度的计算公式如下:
, (14)
式中,dk-1是第k-1次的共轭梯度,βk是一个使相邻两次共轭梯度正交的系数,其计算公式如下:
, (15)
式中,gk是梯度向量,dk-1是共轭梯度向量,符号<,>表示内积运算;
11)计算步长;步长通过线性搜索方法计算,计算公式如下:
, (16)
式中,ε是试探步长,向量d’ cal是当λ=ε时,更新速度模型后,合成的拖缆地震记录,向量dcal是合成的拖缆地震记录,向量drea是观测的拖缆地震记录;试探步长的计算公式如下:
, (17)
式中,向量v表示模型空间中所有点的速度,dk表示共轭梯度向量,max表示从向量中选择最大值,符号||表示绝对值运算,即当向量中元素值为负数时,取其相反数;
12)更新速度模型;速度模型更新的计算公式如下:
, (18)
式中,vk+1是第k次迭代更新后的速度,λk是第k次迭代更新的步长;
13)判断是否满足迭代终止条件;如果满足终止条件,则输出计算结果,否则,将更新后的速度模型作为新的初始速度模型,返回步骤2)。
其中,公式1、3、4、10和11中的波场通过时间2阶、空间12阶精度的规则网格有限差分方法求解,采用完全匹配层边界条件来压制模型边界处的反射波。
其中,所述步骤3)和7)中,反演所用的初始速度模型,通过对真实速度模型平滑得到。
为了更好的说明本发明,下面列举一实施例。
为了进一步说明本方法的实现思路及实现过程并证明方法的有效性,用marmousi2模型进行测试,并和常规全波形反演的结果进行比较。
S1:建立一个真实速度模型,如图2所示。真实速度模型宽度为9216 m,深度为3048m。采用正方形网格离散,网格大小为24 × 24 m。
S2:观测系统:采用海上拖缆数据采集观测系统,从水平方向3360 m处开始以96 m的间隔布设了51个炮点。炮点从左向右依次激发。拖缆上共有140个水听器接收地震记录,道间距24 m,最小偏移距为24 m,最大偏移距为3360 m。11个海底观测站以864 m的间隔均匀布设在水平方向120 m 至 8760 m的海底。地震记录采样时间为3 s,时间采样间隔为2ms。
S3:由真实速度模型(详见图2)和震源函数为10 Hz的雷克子波,利用公式1,通过时间2阶、空间12阶精度的规则网有限差分方法,采用完全匹配层边界条件,正演模拟得到观测的拖缆地震记录。图4是第17炮观测的拖缆地震记录。
S4:由真实速度模型和震源函数为5 Hz的雷克子波,利用公式1,采用同样的正演模拟方法,得到观测的OBS地震记录。图5是第17炮观测的OBS地震记录。
S5:利用公式2,计算权重系数。
S6:由初始速度模型(详见图3)和震源函数为10 Hz的雷克子波,利用公式3和5,采用同样的正演模拟方法,得到预测的拖缆地震记录和拖缆地震记录对应的正传波场。
S7:由初始速度模型和震源函数为5 Hz的雷克子波,利用公式4和6,采用同样的正演模拟方法,得到预测的OBS地震记录和OBS地震记录对应的正传波场。
S8:由观测的拖缆地震记录和预测的拖缆地震记录,利用公式7,计算拖缆地震记录的残差。
S9:由初始速度模型和拖缆地震记录的残差,利用公式10,计算拖缆地震记录对应的反传波场。
S10:由观测的OBS地震记录和预测的OBS地震记录,利用公式8,计算OBS地震记录的残差。
S11:由步骤8和10计算的地震记录残差,利用公式9,计算目标函数值。
S12:由初始速度模型和OBS地震记录的残差,利用公式11,计算OBS地震记录对应的反传波场。
S13:由S5、S6、S8和S10计算的时间域正传和反传波场,利用公式12,计算频率域正传和反传波场,并根据计算的频率波场,利用公式13,计算目标函数的梯度。在本实施例中,共使用了12个频率,分别是拖缆地震记录的6个频率,即6、7.32、8.94、10.91、13.32、16.26Hz,OBS地震记录的6个频率,即2、2.40、2.88、3.46、4.15、4.98 Hz。将12个频率分为6组(详见表1),对于每一个频率组,都迭代15次。
S14:利用公式14计算共轭梯度场。
S15:利用公式16和17计算步长。
S16:利用公式18更新速度模型。判断更新后的速度模型是否满足迭代终止条件,如果满足,则输出结果;如果不满足,则返回步骤5。
图6是只使用拖缆数据,常规全波形反演的结果。和真实速度模型对比,可以发现反演结果只有浅层接近真实速度模型,而600 m以下的区域并没有正确的反演出来。图7是只利用OBS数据的反演结果。和真实速度模型对比,可以发展,反演结果的分辨率较低,反演结果不能满足需要。图8是本发明的反演结果,和图6对比,可以发现,模型600 m以下区域的速度得到了正确的反演,和图7对比,可以发现模型的分辨率得到了显著的提高。
以上所述仅为本发明的较佳实施例而己,并不以本发明为限制,凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等,均应包含在本发明的专利涵盖范围内。
表1 本实施例使用的频率组
Claims (5)
1.一种拖缆地震数据和少量OBS数据联合的波形反演方法,其特征在于,包括如下步骤:
1)构建拖缆地震数据和少量OBS数据联合的混合域全波形反演目标函数;
2)获取野外观测的拖缆和OBS地震记录;
3)对观测地震记录进行预处理;预处理包括带通滤波、F-K滤波等;
4)采用高阶统计方法提取地震子波;
5)计算目标函数中对拖缆地震数据和OBS数据的作用进行控制的权重系数;
6)给定初始速度模型,计算预测的拖缆地震记录及其对应的正传波场;
7)计算OBS地震记录及其对应的正传波场;
8)利用计算和观测的拖缆和OBS地震记录,计算目标函数值;
9)计算拖缆地震记录的残差和拖缆地震记录对应的反传波场;
10)计算OBS地震记录的残差和OBS地震记录对应的反传波场;
11)将步骤6)、7)、9)和10)中计算的时间域波场变换到频率域;
12)计算目标函数的梯度场及共轭梯度场;
13)计算当前迭代更新的步长及对初始模型进行更新的修正量;
14)迭代更新速度模型;
15)判断是否满足迭代终止条件;如果满足终止条件,则输出计算结果,否则,将更新后的速度模型作为新的初始速度模型,返回步骤5)。
2.根据权利要求1所述的一种拖缆地震数据和少量OBS数据联合的波形反演方法,其特征在于,所述步骤1)中,构建的目标函数如下:
, (1)
式中,E是目标函数值,α是权重系数,nr是每一炮水听器的个数,s是炮号,r是水听器号,t是时间,xs是模型空间中炮点s的位置,no是每一炮OBS的个数,o指OBS号,dcal(r,t|xs)是预测的拖缆地震记录,drea(r,t|xs)是观测的拖缆地震记录,dcal(o,t|xs)是预测的OBS地震记录,drea(o,t|xs)是观测的OBS地震记录;本发明通过最小化该目标函数来得到最优的反演结果。
3.根据权利要求1所述的一种拖缆地震数据和少量OBS数据联合的波形反演方法,其特征在于,所述步骤5)中,在反演迭代的初期,主要恢复模型的背景速度场,所以OBS资料应该起主导作用,而在反演后期,主要提高模型的分辨率,所以拖缆数据应该占主要地位;因此,本发明采用下式计算权重系数:
, (2)
式中,k是当前的迭代次数,nite是反演的总迭代次数。
4.根据权利要求1所述的一种拖缆地震数据和少量OBS数据联合的波形反演方法,其特征在于,所述步骤12)中,梯度的计算公式如下:
, (3)
式中,α是权重系数,nr是水听器的个数,no是OBS的个数,ωr是使用拖缆数据时选择的角频率,ωo是使用OBS数据时选择的角频率,s是炮号,x是模型空间的位置,xs是炮点在模型空间的位置,Pfr(x,ωr|xs)是拖缆数据对应的正传波场,Pfo(x,ωo|xs)是OBS数据对应的正传波场,Pbr(x,ωr|xs)是拖缆数据对应的反传波场,Pbo(x,ωo|xs)是OBS数据对应的反传波场,上标*表示复共轭转置。
5.根据权利要求1所述的一种拖缆地震数据和少量OBS数据联合的波形反演方法,其特征在于,从拖缆数据中提取的子波的主频高,而从OBS数据中提取的子波的主频低;因此,步骤6)和9)计算的拖缆地震记录对应的正传和反传波场的主频高,可以用来提高模型的分辨率;步骤7)和10)计算的OBS地震记录对应的正传和反传波场的主频低,可以用来恢复模型的背景速度场;步骤12)中利用拖缆和OBS地震记录对应的正传和反传波场来计算梯度场,克服了传统全波形反演只利用拖缆数据进行反演时,由于缺少低频信息,反演结果不准确的问题。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710864481.4A CN109541681B (zh) | 2017-09-22 | 2017-09-22 | 一种拖缆地震数据和少量obs数据联合的波形反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710864481.4A CN109541681B (zh) | 2017-09-22 | 2017-09-22 | 一种拖缆地震数据和少量obs数据联合的波形反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109541681A true CN109541681A (zh) | 2019-03-29 |
CN109541681B CN109541681B (zh) | 2020-07-17 |
Family
ID=65828237
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710864481.4A Expired - Fee Related CN109541681B (zh) | 2017-09-22 | 2017-09-22 | 一种拖缆地震数据和少量obs数据联合的波形反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109541681B (zh) |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110058302A (zh) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | 一种基于预条件共轭梯度加速算法的全波形反演方法 |
CN111781639A (zh) * | 2020-06-04 | 2020-10-16 | 同济大学 | 针对obs多分量数据的炮检互易弹性波全波形反演方法 |
CN112255685A (zh) * | 2020-09-28 | 2021-01-22 | 广州海洋地质调查局 | 一种obs与海面拖缆地震数据联合成像方法及处理终端 |
CN113466933A (zh) * | 2021-06-11 | 2021-10-01 | 中国海洋大学 | 基于深度加权的地震斜率层析成像方法 |
CN113777654A (zh) * | 2021-08-06 | 2021-12-10 | 同济大学 | 一种基于伴随状态法初至波走时层析的海水速度建模方法 |
CN114488286A (zh) * | 2022-01-25 | 2022-05-13 | 中国海洋大学 | 基于振幅加权的拖缆与海底地震资料联合波形反演方法 |
CN114660655A (zh) * | 2022-03-22 | 2022-06-24 | 西北工业大学青岛研究院 | 一种海底地震仪与拖缆地震联合资料采集的方法 |
CN115469357A (zh) * | 2021-06-11 | 2022-12-13 | 中国海洋大学 | 地震拖缆与obs资料联合叠前最小二乘偏移成像方法 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105445798A (zh) * | 2014-08-21 | 2016-03-30 | 中国石油化工股份有限公司 | 一种基于梯度处理的全波形反演方法和系统 |
CN105549079A (zh) * | 2016-01-12 | 2016-05-04 | 中国矿业大学(北京) | 一种地球物理参数的全波形反演模型的建立方法和装置 |
CN105676277A (zh) * | 2015-12-30 | 2016-06-15 | 中国石油大学(华东) | 一种提高高陡构造速度反演效率的全波形联合反演方法 |
US20160216389A1 (en) * | 2015-01-23 | 2016-07-28 | Advanced Geophysical Technology Inc. | Beat tone full waveform inversion |
CN106054244A (zh) * | 2016-06-16 | 2016-10-26 | 吉林大学 | 截断时窗的低通滤波多尺度全波形反演方法 |
CN106405651A (zh) * | 2016-11-14 | 2017-02-15 | 中国石油化工股份有限公司 | 一种基于测井匹配的全波形反演初始模型构建方法 |
CN106950596A (zh) * | 2017-04-11 | 2017-07-14 | 中国石油大学(华东) | 一种基于子波迭代估计的有限差分对比源全波形反演方法 |
-
2017
- 2017-09-22 CN CN201710864481.4A patent/CN109541681B/zh not_active Expired - Fee Related
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105445798A (zh) * | 2014-08-21 | 2016-03-30 | 中国石油化工股份有限公司 | 一种基于梯度处理的全波形反演方法和系统 |
US20160216389A1 (en) * | 2015-01-23 | 2016-07-28 | Advanced Geophysical Technology Inc. | Beat tone full waveform inversion |
CN105676277A (zh) * | 2015-12-30 | 2016-06-15 | 中国石油大学(华东) | 一种提高高陡构造速度反演效率的全波形联合反演方法 |
CN105549079A (zh) * | 2016-01-12 | 2016-05-04 | 中国矿业大学(北京) | 一种地球物理参数的全波形反演模型的建立方法和装置 |
CN106054244A (zh) * | 2016-06-16 | 2016-10-26 | 吉林大学 | 截断时窗的低通滤波多尺度全波形反演方法 |
CN106405651A (zh) * | 2016-11-14 | 2017-02-15 | 中国石油化工股份有限公司 | 一种基于测井匹配的全波形反演初始模型构建方法 |
CN106950596A (zh) * | 2017-04-11 | 2017-07-14 | 中国石油大学(华东) | 一种基于子波迭代估计的有限差分对比源全波形反演方法 |
Non-Patent Citations (1)
Title |
---|
林朋 等: "基于共轭梯度法的全波形反演", 《煤田地质与勘探》 * |
Cited By (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110058302A (zh) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | 一种基于预条件共轭梯度加速算法的全波形反演方法 |
CN111781639A (zh) * | 2020-06-04 | 2020-10-16 | 同济大学 | 针对obs多分量数据的炮检互易弹性波全波形反演方法 |
CN112255685A (zh) * | 2020-09-28 | 2021-01-22 | 广州海洋地质调查局 | 一种obs与海面拖缆地震数据联合成像方法及处理终端 |
CN113466933A (zh) * | 2021-06-11 | 2021-10-01 | 中国海洋大学 | 基于深度加权的地震斜率层析成像方法 |
CN115469357A (zh) * | 2021-06-11 | 2022-12-13 | 中国海洋大学 | 地震拖缆与obs资料联合叠前最小二乘偏移成像方法 |
CN115469357B (zh) * | 2021-06-11 | 2024-09-24 | 中国海洋大学 | 地震拖缆与obs资料联合叠前最小二乘偏移成像方法 |
CN113777654A (zh) * | 2021-08-06 | 2021-12-10 | 同济大学 | 一种基于伴随状态法初至波走时层析的海水速度建模方法 |
CN114488286A (zh) * | 2022-01-25 | 2022-05-13 | 中国海洋大学 | 基于振幅加权的拖缆与海底地震资料联合波形反演方法 |
CN114488286B (zh) * | 2022-01-25 | 2023-03-10 | 中国海洋大学 | 基于振幅加权的拖缆与海底地震资料联合波形反演方法 |
CN114660655A (zh) * | 2022-03-22 | 2022-06-24 | 西北工业大学青岛研究院 | 一种海底地震仪与拖缆地震联合资料采集的方法 |
CN114660655B (zh) * | 2022-03-22 | 2024-07-26 | 西北工业大学青岛研究院 | 一种海底地震仪与拖缆地震联合资料采集的方法 |
Also Published As
Publication number | Publication date |
---|---|
CN109541681B (zh) | 2020-07-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109541681A (zh) | 一种拖缆地震数据和少量obs数据联合的波形反演方法 | |
CN107203002B (zh) | 反演速度模型的建立方法和地下结构的像的获得方法 | |
US8560242B2 (en) | Pseudo-analytical method for the solution of wave equations | |
CN106405651B (zh) | 一种基于测井匹配的全波形反演初始速度模型构建方法 | |
CN108802813B (zh) | 一种多分量地震资料偏移成像方法及系统 | |
CN108241173B (zh) | 一种地震资料偏移成像方法及系统 | |
CN110579795B (zh) | 基于被动源地震波形及其逆时成像的联合速度反演方法 | |
CN113156493B (zh) | 一种使用归一化震源的时频域全波形反演方法及装置 | |
CN109507726A (zh) | 时间域弹性波多参数全波形的反演方法及系统 | |
CN110007340A (zh) | 基于角度域直接包络反演的盐丘速度密度估计方法 | |
KR20170009609A (ko) | 반복적 직접 파형 역산 및 완전 파형 역산을 이용한 탄성파 영상화 장치 및 방법 | |
KR101820850B1 (ko) | 직접 파형 역산의 반복 적용을 이용한 탄성파 영상화 장치 및 방법 | |
CN107356972B (zh) | 一种各向异性介质的成像方法 | |
CN109633749A (zh) | 基于散射积分法的非线性菲涅尔体地震走时层析成像方法 | |
Qu et al. | 3-D Least-Squares Reverse Time Migration in Curvilinear-τ Domain | |
CN113740906A (zh) | 一种水下垂直缆地震波干涉成像方法及装置 | |
CN108363097B (zh) | 一种地震资料偏移成像方法 | |
CN106855638A (zh) | 一种匹配追踪地震谱分解方法及装置 | |
CN112230274B (zh) | 面向随钻导向的声波方程频率域逆时偏移快速成像方法 | |
Yang et al. | Full waveform inversion of combined towed streamer and limited OBS seismic data: a theoretical study | |
CN113050163B (zh) | 一种振幅、相位信息可调节的全波形反演方法及装置 | |
CN107526102B (zh) | 纵波与转换波联合偏移速度建模方法和装置 | |
CN114721044A (zh) | 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统 | |
CN111208568B (zh) | 一种时间域多尺度全波形反演方法及系统 | |
CN111175822B (zh) | 改进直接包络反演与扰动分解的强散射介质反演方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200717 Termination date: 20210922 |