地震数据保幅重建方法和系统
技术领域
本发明涉及商用厨房灭火技术领域,尤其是涉及地震数据保幅重建方法和系统。
背景技术
在地震数据重建方面,一般有如下方法:基于预测滤波重建方法,基于波动方程重建方法和基于函数变换重建方法等。基于函数变换这类方法的特点是利用数学变换理论和信号分析原理,如Fourier变换、Wavelet变换、Curvelet变换、Radon变换和压缩感知等对数据进行重建。
Radon变换又包括线性Radon、双曲Radon和抛物Radon这三种。在地震勘探中,利用水平层状介质或小倾角地质结构的反射同相轴走时特征在CMP道集上是双曲线的性质,双曲Radon变换在地震数据多次波压制、插值重建理论研究方面得到广泛关注,但是其计算效率的低下限制了双曲Radon变换在工业上的应用。
在基于抛物Radon变换地震数据重建方面,传统的抛物Radon变换是同相轴振幅沿着不同的曲率进行叠加,没有考虑振幅随偏移距的变化特征,因此在处理实际地震数据时不具备保幅性;同时,以上方法中的迭代加权系数在不同位置的缺失道都相同,而实际上不同位置缺失的重建效果是不一样的,用相同的权系数会导致较大误差出现。
发明内容
有鉴于此,本发明的目的在于提供地震数据保幅重建方法和系统,有效改善了远偏移距缺失数据的重建效果,具有保幅性好、计算效率高的优点。
第一方面,本发明实施例提供了地震数据保幅重建方法包括:
获取原始缺失地震数据,并对所述原始缺失地震数据进行预处理得到动校正数据;
对所述动校正数据进行频域变换得到第一频域数据;
利用快速高阶抛物Radon变换算法对所述第一频域数据进行计算得到第二频域数据;
计算权系数,并对所述第二频域数据和所述权系数进行加权计算得到频域重建数据;
判断所述频域重建数据是否满足预设条件;
如果不满足所述预设条件,则将所述频域重建数据进行迭代计算得到时域重建数据。
结合第一方面,本发明实施例提供了第一方面的第一种可能的实施方式,其中,所述频域变换包括傅里叶变换,所述利用快速高阶抛物Radon变换算法对所述第一频域数据进行计算得到第二频域数据包括:
根据快速高阶抛物Radon反变换公式对所述第一频域数据进行反演求解得到变换域系数,其中,所述反演求解包括利用最小二乘框架下的奇异值分解SVD方法进行矩阵求逆计算;
将所述变换域系数代入所述快速高阶抛物Radon反变换公式得到所述第二频域数据。
结合第一方面,本发明实施例提供了第一方面的第二种可能的实施方式,其中,所述计算权系数包括:
根据下式计算权系数:
其中,W(Tn)为所述权系数,Tn为第n道数据,xn,yn为缺失道的坐标;xmax,ymax为最大偏移距;γ1和γ2取值范围为[0,1];an表示所述缺失道所处位置的缺失程度。
结合第一方面,本发明实施例提供了第一方面的第三种可能的实施方式,其中,所述预设条件为所述频域重建数据满足所述迭代计算的次数要求,所述将所述频域重建数据进行迭代计算包括:
将所述频域重建数据填补入原始缺失道;
重复进行所述快速高阶抛物Radon变换算法和所述加权计算,直至满足所述次数要求,其中,所述次数要求为三次;
对满足所述次数的所述频域重建数据进行傅里叶反变换得到所述时域重建数据。
结合第一方面的第一种可能的实施方式,本发明实施例提供了第一方面的第四种可能的实施方式,其中,根据快速高阶抛物Radon反变换公式对所述第一频域数据进行反演求解得到变换域系数包括:
根据下式计算所述变换域系数:
其中,M(f,λx,λy)为所述变换域系数,分别为X,Y方向的所述快速高阶抛物Radon变换算子,H为矩阵的共轭转置,D1(f,xn,yk)为所述第一频域数据。
第二方面,本发明实施例提供了地震数据保幅重建系统,包括:
预处理单元,用于获取原始缺失地震数据,并对所述原始缺失地震数据进行预处理得到动校正数据;
变换单元,用于对所述动校正数据进行频域变换得到第一频域数据;
第一计算单元,用于利用快速高阶抛物Radon变换算法对所述第一频域数据进行计算得到第二频域数据;
第二计算单元,计算权系数,并对所述第二频域数据和所述权系数进行加权计算得到频域重建数据;
判断单元,用于判断所述频域重建数据是否满足预设条件;
迭代单元,用于在不满足所述预设条件的情况下,则将所述频域重建数据进行迭代计算得到时域重建数据。
结合第二方面,本发明实施例提供了第二方面的第一种可能的实施方式,其中,所述频域变换包括傅里叶变换,所述第一计算单元包括:
根据快速高阶抛物Radon反变换公式对所述第一频域数据进行反演求解得到变换域系数,其中,所述反演求解包括利用最小二乘框架下的奇异值分解SVD方法进行矩阵求逆计算;
将所述变换域系数代入所述快速高阶抛物Radon反变换公式得到所述第二频域数据。
结合第二方面,本发明实施例提供了第二方面的第二种可能的实施方式,其中,所述第二计算单元包括:
根据下式计算权系数:
其中,W(Tn)为所述权系数,Tn为第n道数据,xn,yn为缺失道的坐标;xmax,ymax为最大偏移距;γ1和γ2取值范围为[0,1];an表示所述缺失道所处位置的缺失程度。
结合第二方面,本发明实施例提供了第二方面的第三种可能的实施方式,其中,所述满足预设条件为满足所述迭代计算的次数要求,所述迭代单元包括:
将所述频域重建数据填补入原始缺失道;
重复进行所述快速算法和所述加权计算,直至满足所述次数要求,其中,所述次数要求为三次;
对满足所述次数的所述频域重建数据进行傅里叶反变换得到所述时域重建数据。
结合第二方面的第一种可能的实施方式,本发明实施例提供了第二方面的第四种可能的实施方式,其中,所述第一计算单元还包括:
根据下式计算所述变换域系数:
其中,M(f,λx,λy)为所述变换域系数,分别为X,Y方向的所述快速高阶抛物Radon变换算子,H为矩阵的共轭转置,D1(f,xn,yk)为所述第一频域数据。
本发明提供地震数据保幅重建方法和系统,所述方法包括:获取原始缺失地震数据,并对原始缺失地震数据进行预处理得到动校正数据;对动校正数据进行频域变换得到第一频域数据;利用快速高阶抛物Radon变换算法对第一频域数据进行计算得到第二频域数据;计算权系数,并对第二频域数据和权系数进行加权计算得到频域重建数据;判断频域重建数据是否满足预设条件;如果不满足预设条件,则将频域重建数据进行迭代计算得到时域重建数据。本发明可以改进远偏移距数据的重建效果,并提高方法的计算效率。
本发明的其他特征和优点将在随后的说明书中阐述,并且,部分地从说明书中变得显而易见,或者通过实施本发明而了解。本发明的目的和其他优点在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。
为使本发明的上述目的、特征和优点能更明显易懂,下文特举较佳实施例,并配合所附附图,作详细说明如下。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案,下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图是本发明的一些实施方式,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明实施例一提供的地震数据保幅重建方法流程图;
图2为本发明实施例一提供的步骤S103的方法流程图;
图3为本发明实施例一提供的步骤S1061的方法流程图;
图4为本发明实施例二提供的地震数据保幅重建系统结构示意图;
图5为本发明实施例三提供的地震数据保幅重建方法流程图;
图6(a)为本发明实施例三提供的传统方法重建效果图;
图6(b)为本发明实施例三提供的本发明方法重建效果图;
图7(a)为本发明实施例三提供的近偏移距缺失第54道重建结果比较图;
图7(b)为本发明实施例三提供的中偏移距缺失第79道重建结果比较图;
图7(c)为本发明实施例三提供的远偏移距缺失第106道重建结果比较图。
图标:
10-预处理单元;20-变换单元;30-第一计算单元;40-第二计算单元;50-判断单元;60-迭代单元。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合附图对本发明的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
目前,在基于抛物Radon变换地震数据重建方面,传统的抛物Radon变换是同相轴振幅沿着不同的曲率进行叠加,没有考虑振幅随偏移距的变化特征,因此在处理实际地震数据时不具备保幅性;同时,以上方法中的迭代加权系数在不同位置的缺失道都相同,而实际上不同位置缺失的重建效果是不一样的,用相同的权系数会导致较大误差出现。基于此,本发明实施例提供的地震数据保幅重建方法和系统,可以改进远偏移距数据的重建效果,并提高方法的计算效率。
实施例一:
图1为本发明实施例一提供的地震数据保幅重建方法流程图。
参照图1,地震数据保幅重建方法包括:
步骤S101,获取原始缺失地震数据,并对原始缺失地震数据进行预处理得到动校正数据;
步骤S102,对动校正数据进行频域变换得到第一频域数据;
这里,第一频域数据是将时空域数据变换到频率-空间域数据。
步骤S103,利用快速高阶抛物Radon变换算法对第一频域数据进行计算得到第二频域数据;
步骤S104,计算权系数,并对第二频域数据和权系数进行加权计算得到频域重建数据;
步骤S105,判断频域重建数据是否满足预设条件;
步骤S1061,如果不满足预设条件,则将频域重建数据进行迭代计算得到时域重建数据。
具体地,如果满足预设条件,则执行步骤S1062,直接将频域重建数据进行傅里叶反变换得到时域重建数据。
根据本发明的示例性实施例,频域变换包括傅里叶变换,利用快速高阶抛物Radon变换算法对第一频域数据进行计算得到第二频域数据包括:
如图2所示,步骤S201,根据快速高阶抛物Radon反变换公式对第一频域数据进行反演求解得到变换域系数,其中,反演求解包括利用最小二乘框架下的奇异值分解SVD方法进行矩阵求逆计算;
步骤S202,将变换域系数代入快速高阶抛物Radon反变换公式得到第二频域数据。
根据本发明的示例性实施例,计算权系数包括:
根据公式(1)计算权系数:
其中,W(Tn)为权系数,Tn为第n道数据,xn,yn为缺失道的坐标;xmax,ymax为最大偏移距;γ1和γ2取值范围为[0,1];an表示缺失道所处位置的缺失程度。
根据本发明的示例性实施例,预设条件为频域重建数据满足迭代计算的次数要求,将频域重建数据进行迭代计算包括:
如图3所示,步骤S301,将频域重建数据填补入原始缺失道;
步骤S302,重复进行快速算法和加权计算,直至满足次数要求,其中,次数要求为三次;
步骤S303,对满足次数的频域重建数据进行傅里叶反变换得到时域重建数据。
具体地,判断频域重建数据是否满足迭代次数要求,在不满足的情况下执行步骤将其带入原始数据缺失位置,并重复快速算法和加权计算,直到迭代次数达到3次为止,得到满足条件的频域重建数据,并对其进行傅里叶反变换得到时域重建数据。
根据本发明的示例性实施例,根据快速高阶抛物Radon反变换公式对第一频域数据进行反演求解得到变换域系数包括:
根据公式(2)计算变换域系数:
其中,M(f,λx,λy)为所述变换域系数,分别为X,Y方向的所述快速高阶抛物Radon变换算子,H为矩阵的共轭转置,D1(f,xn,yk)为所述第一频域数据。
本发明提供地震数据保幅重建方法和系统,所述方法包括:获取原始缺失地震数据,并对原始缺失地震数据进行预处理得到动校正数据;对动校正数据进行频域变换得到第一频域数据;利用快速高阶抛物Radon变换算法对第一频域数据进行计算得到第二频域数据;计算权系数,并对第二频域数据和权系数进行加权计算得到频域重建数据;判断频域重建数据是否满足预设条件;如果不满足预设条件,则将频域重建数据进行迭代计算得到时域重建数据。本发明可以改进远偏移距数据的重建效果,并提高方法的计算效率。
实施例二:
图4为本发明实施例二提供的地震数据保幅重建系统结构示意图。
参照图4,地震数据保幅重建系统包括:
预处理单元10,用于获取原始缺失地震数据,并对原始缺失地震数据进行预处理得到动校正数据;
变换单元20,用于对动校正数据进行频域变换得到第一频域数据;
第一计算单元30,用于利用快速高阶抛物Radon变换算法对第一频域数据进行计算得到第二频域数据;
第二计算单元40,计算权系数,并对第二频域数据和权系数进行加权计算得到频域重建数据;
判断单元50,用于判断频域重建数据是否满足预设条件;
迭代单元60,用于在不满足预设条件的情况下,则将频域重建数据进行迭代计算得到时域重建数据。
根据本发明的示例性实施例,频域变换包括傅里叶变换,第一计算单元包括:
根据快速高阶抛物Radon反变换公式对第一频域数据进行反演求解得到变换域系数,其中,反演求解包括利用最小二乘框架下的奇异值分解SVD方法进行矩阵求逆计算;
将变换域系数代入快速高阶抛物Radon反变换公式得到第二频域数据。
根据本发明的示例性实施例,第二计算单元包括:
根据公式(1)计算权系数。
根据本发明的示例性实施例,满足要求为满足迭代计算的次数要求,迭代单元包括:
将频域重建数据填补入原始缺失道;
重复进行快速算法和加权计算,直至满足次数要求,其中,次数要求为三次;
对满足次数的频域重建数据进行傅里叶反变换得到时域重建数据。
根据本发明的示例性实施例,第一计算单元还包括:
根据公式(2)计算变换域系数。
实施例三:
图5为本发明实施例三提供的地震数据保幅重建方法流程图。
步骤S401,获取时空域三维缺失地震数据,即原始缺失地震数据,对原始缺失地震数据进行预处理,得到部分动校正数据记为D0(t,xn,yk),n=1,2,3,…,N;k=1,2,3,…,K;N表示inline方向每条测线的道数,K表示crossline方向每条测线的道数。
步骤S402,对于部分动校正后数据通过频域变换,即傅里叶变换,将时空域数据变换到频率-空间域数据D1(f,xn,yk),即第一频域数据,f表示频率。
步骤S403,在数据频带范围内,由快速高阶抛物Radon反变换公式,如公式(3)所示来反演求解地震数据的高阶抛物Radon正变换域系数,即变换域系数,在最小二乘框架下用SVD方法求其解如公式(2)所示,每个频率数据独立计算,直到频带范围内的所有频率都完成计算;
其中,M(f,λx,λy)是待求的高阶抛物Radon变换域系数,分别为X,Y方向的快速高阶抛物Radon变换算子,表达式如公式(4)和公式(5)所示:
其中,p(x),p(y)是正交多项式基函数,它们的求取方法相同,一般取2阶正交多项式,以pj(x)为例,pj(x)的求取如公式(6)所示,系数α的求取如公式(7)至公式(10)所示:
P0=1/α00(10)
需要说明的是,系数αjk和多项式Pj应按照以下顺序计算:α10,α11,P1,α20,α21,α22,P2,…,αL0,αL1,…,αLL,PL.。另外,λx和λy的取值范围如式(11)和式(12)所示:
非快速高阶抛物Radon变换算子如公式(13)(14)所示:
其中,ω=2πf,由公式(4)(5)(13)(14)可以看出,快速高阶抛物Radon变换算子参数只与偏移距有关,与频率f无关,而非快速高阶抛物Radon变换算子与频率相关,不同频率数据的变换算子都不同.因此,求解不同频率的正变换系数矩阵M(f)时,快速高阶抛物Radon变换算子L(λx)、L(λy)均相同,变换算子L(λx)、L(λy)求逆过程仅需要计算一次,而非快速方法中的求逆次数是数据的频带宽度,地震数据的有效频带宽度一般是5-60Hz,所以求逆次数不少于50次。很明显,快速算法能够显著地提高计算效率。
公式(2)中的H表示矩阵的共轭转置,所有的矩阵求逆采用SVD的方法,具体方法如下:
将待求逆矩阵统一用A表示,矩阵A的奇异值分解为:A=USV T则矩阵A的逆可以用公式(15)计算:
A-1=V[(S 2+Λ)-1S]U T (15)
其中,Λ是阻尼因子,一般取值为最大奇异值的1%,但是当最小奇异值与最大奇异值的比小于1%时,阻尼因子Λ的值取为Smax*0.01-Smin,这里的Smax,Smin分别为最大奇异值和最小奇异值.采用奇异值分解法对高阶抛物Radon变换进行求解,可以提高Radon域系数的精度,得到高分辨率的解。
步骤S404,由步骤S403得到所有频率数据计算完毕得到的高阶Radon域系数M,即变换域系数,将其带入快速高阶抛物Radon反变换公式,如公式(16)计算得到频率空间域数据,即第二频域系数,参与计算的同样是频带范围内的所有频率数据。
步骤S405,将步骤S404中缺失位置重建出的数据振幅经过加权,加权系数为W,得到新的重建数据,即频域重建数据,根据步骤S406判断频域重建数据是否满足迭代次数要求,在不满足的情况下执行步骤S4071,即将其带入原始数据缺失位置,并重复步骤S403至步骤S406,直到迭代次数达到3次为止。每一缺失道的权系数与缺失道所处位置以及缺失程度相关,计算如公式(1)所示:
其中,Tn表示第n道数据,xn,yn是缺失道的坐标;xmax,ymax为最大偏移距,如果数据是双边接收则为缺失道所在那一边的最大偏移距;γ1和γ2取值范围为[0,1];an表示缺失道所处位置的缺失程度,即缺失道所处位置在inline和crossline这两个方向连续缺失的最多道数;由于任何方法都会受到数据缺失空间大小的限制,所以an的取值是有一定范围的,在本发明中最大的缺失程度由公式(17)确定:
其中,λ为待重建数据主频对应的波长;Δs为3D数据的面元大小,如果是2D数据则为道间距;在本发明中k可以取到10;可以看出,数据采集得越密集Δs越小,可有效重建的连续缺失道数越多,这与实际情况相符合.与缺失程度an相关的比例系数0.05,是经过试验和测试得到的经验值。
步骤S4072,将步骤S406中得到满足要求的f-x-y域数据经过傅里叶反变换得到t-x-y域重建后数据。
现用实际数据来测试本发明实施例提出的方法在重建中的有效性,该数据有45条测线,线间距为200m,每条测线113道,道间距为100m,时间采样间隔为8ms,最远偏移距为5700m。
(1)重建效果的改进:
将该数据随机缺失30%,用快速加权3D高阶抛物Radon变换对其重建,任意抽取其中一条测线(inline10)检测该方法的有效性。如图6所示为随机缺失重建结果,这里图6包括图6(a)和图6(b),重建前如图6(a)所示,重建后如图6(b)所示,数据同相轴连续,相位及振幅和原始数据一致。同时比较了本发明实施例的重建效果和传统方法(不考虑缺失位置和缺失程度对重建结果的影响)重建效果,如图7所示,这里图7包括图7(a)、图7(b)、图7(c),图7(a)为inline10近偏移距缺失道trace54重建效果,图7(b)为中偏移距缺失道trace79重建效果,图7(c)为远偏移距缺失道trace106重建效果,从图中可以发现在近偏移距二者重建效果一致,但是在远偏移距我们的方法效果更好,重建道与原始道振幅差异更小。
(2)重建效率的提高:
表1给出了该地震数据在不同单频重建次数和不同变换算子大小的情况下,利用3D快速高阶抛物Radon变换(FHOPRT)和3D非快速高阶抛物Radon变换(HOPRT)的计算时间,单位是“秒”。由表1可以看出该地震数据在参数采样个数为100、单频重建次数为22的情况下,利用快速高阶算法计算效率提高了4倍以上;参数采样个数为300、单频重建次数为83的情况下,利用快速高阶算法计算效率提高了19倍以上.比较表1中不同参数采样个数和单频重建次数情况下3D HOPRT与3D FHOPRT计算时间的比值可以发现,参数采样个数和单频重建次数越多,快速算法相对于非快速算法提高的计算效率越高。
表1单频重建次数和变换算子大小对计算效率的联合影响
本发明实施例提供的地震数据保幅重建方法,与上述实施例提供的地震数据保幅重建方法具有相同的技术特征,所以也能解决相同的技术问题,达到相同的技术效果。
本发明实施例所提供的地震数据保幅重建方法和系统的计算机程序产品,包括存储了程序代码的计算机可读存储介质,所述程序代码包括的指令可用于执行前面方法实施例中所述的方法,具体实现可参见方法实施例,在此不再赘述。
所属领域的技术人员可以清楚地了解到,为描述的方便和简洁,上述描述的系统和装置的具体工作过程,可以参考前述方法实施例中的对应过程,在此不再赘述。
本发明实施例还提供一种电子设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,处理器执行计算机程序时实现上述实施例提供的地震数据保幅重建方法的步骤。本发明实施例还提供一种计算机可读存储介质,计算机可读存储介质上存储有计算机程序,计算机程序被处理器运行时执行上述实施例的地震数据保幅重建方法的步骤
所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时,可以存储在一个计算机可读取存储介质中。基于这样的理解,本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来,该计算机软件产品存储在一个存储介质中,包括若干指令用以使得一台计算机设备(可以是个人计算机,服务器,或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括:U盘、移动硬盘、只读存储器(ROM,Read-Only Memory)、随机存取存储器(RAM,Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
此外,术语“第一”、“第二”、“第三”仅用于描述目的,而不能理解为指示或暗示相对重要性。
最后应说明的是:以上所述实施例,仅为本发明的具体实施方式,用以说明本发明的技术方案,而非对其限制,本发明的保护范围并不局限于此,尽管参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化,或者对其中部分技术特征进行等同替换;而这些修改、变化或者替换,并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应所述以权利要求的保护范围为准。