CN106257309A - 叠后地震数据体处理方法及装置 - Google Patents

叠后地震数据体处理方法及装置 Download PDF

Info

Publication number
CN106257309A
CN106257309A CN201610058962.1A CN201610058962A CN106257309A CN 106257309 A CN106257309 A CN 106257309A CN 201610058962 A CN201610058962 A CN 201610058962A CN 106257309 A CN106257309 A CN 106257309A
Authority
CN
China
Prior art keywords
seismic data
seismic
inverse
sample
value
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
Application number
CN201610058962.1A
Other languages
English (en)
Other versions
CN106257309B (zh
Inventor
吴吉忠
穆立华
赵宝银
石文武
邢文军
杨晓利
郭爱华
张建坤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Petrochina Co Ltd
Original Assignee
Petrochina Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Petrochina Co Ltd filed Critical Petrochina Co Ltd
Priority to CN201610058962.1A priority Critical patent/CN106257309B/zh
Publication of CN106257309A publication Critical patent/CN106257309A/zh
Application granted granted Critical
Publication of CN106257309B publication Critical patent/CN106257309B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/303Analysis for determining velocity profiles or travel times

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

本发明公开一种叠后地震数据体处理方法及装置,属于地震数据处理领域。方法包括:从叠后地震数据体对应的地震测线中选取M条样本测线,根据M条样本测线对应的校正Q值场和M条样本测线对应的加权系数场确定M条样本测线对应的最终Q值场,根据M条样本测线对应的样本速度场和最终Q值场更新李庆忠经验公式,采用更新后的李庆忠经验公式确定地震工区的Q值场;向GPU发送地震工区的Q值场和叠后地震数据体以便于GPU根据地震工区的Q值场对叠后地震数据体进行处理。本发明解决了叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。本发明用于叠后地震数据体的处理。

Description

叠后地震数据体处理方法及装置
技术领域
本发明涉及地震数据处理领域,特别涉及一种叠后地震数据体处理方法及装置。
背景技术
地震勘探是一种利用地下介质弹性和密度的差异,通过观测和分析大地对人工激发地震波的响应,推断地下岩层的性质和形态的地球物理勘探方法。其中,可以采用采集设备对地震波进行采集得到地震数据体,地震数据体可以叠加形成叠后地震数据体。在地震勘探中,为了便于观测和分析,通常需要对叠后地震数据体进行处理以提高叠后地震数据体的分辨率。
现有技术中,可以采用反Q滤波法对叠后地震数据体进行处理,以提高叠后地震数据体的分辨率。具体地,可以从已知的工区叠加速度场中获取地震数据体对应的速度场v,该速度场v中包括地震数据体对应的地震波在地层上各点的传播速度,然后根据该速度场v采用李庆忠经验公式Q=α·vβ计算得到地层上各点的Q值,进而根据地层上各点的Q值对叠后地震数据体进行反Q滤波处理。
在实现本发明的过程中,发明人发现现有技术至少存在以下问题:
现有技术中采用李庆忠经验公式确定Q值,由于李庆忠经验公式中的α、β为固定值,针对不同的地质状况,采用李庆忠经验公式确定的Q值与实际Q值相差较大,因此,现有技术对叠后地震数据体处理的灵活性较低、可靠性较差。
发明内容
为了解决现有技术中的问题,本发明提供一种叠后地震数据体处理方法及装置。所述技术方案如下:
第一方面,提供一种叠后地震数据体处理方法,用于中央处理器CPU,所述方法包括:
从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线,所述M为大于或者等于2的整数,且M=0.003N,所述N为所述叠后地震数据体的地震主测线的条数,每条所述地震测线对应一个地震剖面;
从所述地震工区的叠加速度场中确定所述M条样本测线中的每条样本测线对应的速度场,得到所述M条样本测线对应的样本速度场ν;
根据所述M条样本测线对应的样本速度场ν,采用李庆忠经验公式确定所述M条样本测线对应的初始Q值场Q0,所述李庆忠经验公式为Q=α·νβ,α=14,β=2.2;
对所述M条样本测线对应的初始Q值场Q0进行校正,得到所述M条样本测线对应的校正Q值场Q1
根据所述M条样本测线对应的校正Q值场Q1和所述M条样本测线对应的加权系数场λfinal,确定所述M条样本测线对应的最终Q值场Q2
根据所述M条样本测线对应的样本速度场ν和所述M条样本测线对应的最终Q值场Q2,采用最小二乘法更新所述李庆忠经验公式,得到更新后的李庆忠经验公式;
根据所述地震工区的叠加速度场,采用所述更新后的李庆忠经验公式确定所述地震工区的Q值场Qall
向图形处理器GPU发送所述地震工区的Q值场Qall和所述叠后地震数据体,以便于所述GPU根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
进一步地,在所述根据所述M条样本测线对应的校正Q值场Q1和所述M条样本测线对应的加权系数场λfinal,确定所述M条样本测线对应的最终Q值场Q2之前,所述方法还包括:
根据预设的加权系数组对所述M条样本测线对应的校正Q值场Q1进行处理,得到所述M条样本测线对应的加权Q值场组,所述预设的加权系数组中包括P个加权系数λ,所述加权Q值场组中包括P套加权Q值场,所述P为大于或者等于1的整数;
根据所述加权Q值场组对所述M条样本测线组成的样本地震数据体进行反Q滤波处理,得到反Q滤波数据体组,所述反Q滤波数据体组中包括P套反Q滤波数据体,每套反Q滤波数据体包括M条反Q样本测线,每条所述反Q样本测线对应一个地震剖面;
采用傅里叶变换对所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析;
根据频谱分析结果进行处理,确定所述M条反Q样本测线对应的加权系数场;
将所述M条反Q样本测线对应的加权系数场确定为所述M条样本测线对应的加权系数场λfinal
进一步地,所述采用傅里叶变换对所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析,包括:
在所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面上选择三个矩形时窗,所述三个矩形时窗包括:浅层矩形时窗、中层矩形时窗和深层矩形时窗;
采用傅里叶变换对所述每条反Q样本测线对应的地震剖面上位于所述三个矩形时窗中的每个矩形时窗内的部分进行处理,得到所述每条反Q样本测线对应的3个反Q样本频谱;
根据所述每套反Q滤波数据体中的M条反Q样本测线对应的反Q样本频谱,得到所述每套反Q滤波数据体对应的3M个反Q样本频谱;
根据所述反Q滤波数据体组中的P套反Q滤波数据体对应的反Q样本频谱得到3×M×P个反Q样本频谱,所述3×M×P个反Q样本频谱中包括:M×P个浅层矩形时窗对应的反Q样本频谱、M×P个中层矩形时窗对应的反Q样本频谱和M×P个深层矩形时窗对应的反Q样本频谱;
对所述3×M×P个反Q样本频谱中的每个反Q样本频谱进行频谱分析。
进一步地,所述根据频谱分析结果进行处理,确定所述M条反Q样本测线对应的加权系数场,包括:
根据频谱分析结果确定所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱,所述第一目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述浅层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,所述第二目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述中层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,所述第三目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述深层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱;
分别确定所述第一目标频谱对应的矩形时窗对应的第一目标加权系数、所述第二目标频谱对应的矩形时窗对应的第二目标加权系数和所述第三目标频谱对应的矩形时窗对应的第三目标加权系数,得到所述每条反Q样本测线对应的3个目标加权系数,其中,所述第一目标加权系数、所述第二目标加权系数和所述第三目标加权系数都为所述预设的加权系数组中的加权系数;
根据所述M条反Q样本测线中的每条反Q样本测线对应的3个目标加权系数,得到所述M条反Q样本测线对应的3M个目标加权系数;
对所述3M个目标加权系数进行空间插值,得到所述地震工区对应的加权系数场;
对所述地震工区对应的加权系数场进行平滑处理,得到处理后的加权系数场;
从所述处理后的加权系数场中抽取出所述M条反Q样本测线对应的加权系数场。
进一步地,在所述向图形处理器GPU发送所述地震工区的Q值场Qall和所述叠后地震数据体之前,所述方法还包括:
按照预设采样间隔对所述地震工区的Q值场Qall进行采样,得到采样Q值场Qs
对所述叠后地震数据体进行分块处理,得到D个地震数据块,所述D为大于或者等于2的整数;
所述向图形处理器GPU发送所述地震工区的Q值场Qall和叠后地震数据体,包括:
依次向所述GPU发送所述采样Q值场Qs和所述D个地震数据块中的每个地震数据块。
进一步地,所述对所述叠后地震数据体进行分块处理,得到D个地震数据块,包括:
计算所述叠后地震数据体中的地震数据的总道数E;
根据所述GPU在一个周期内处理的地震数据的道数Tracemax和所述叠后地震数据体中的地震数据的总道数E,对所述叠后地震数据体进行分块处理,得到所述D个地震数据块;
其中,E=Tracemax×(D-1)+L,所述L为所述D个地震数据块中的最后一个地震数据块中的地震数据的道数。
第二方面,提供一种叠后地震数据体处理方法,用于图形处理器GPU,所述方法包括:
接收中央处理器CPU发送的地震工区的Q值场Qall和叠后地震数据体;
根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
进一步地,所述接收中央处理器CPU发送的地震工区的Q值场Qall和叠后地震数据体,包括:
接收所述CPU发送的采样Q值场Qs和D个地震数据块中的每个地震数据块,所述采样Q值场Qs是所述CPU按照预设采样间隔对所述地震工区的Q值场Qall进行采样得到的,所述D个地震数据块是所述CPU对所述叠后地震数据体进行分块处理得到的;
所述根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体,包括:
根据所述采样Q值场Qs对所述D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块;
根据所述D个处理后的地震数据块,得到所述处理后的叠后地震数据体。
进一步地,所述根据所述采样Q值场Qs对所述D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块,包括:
对所述采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ
根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到所述D个处理后的地震数据块。
进一步地,所述对所述采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ,包括:
根据所述采样Q值场Qs和Q值插值计算公式确定每个时间深度τ处的Q值Qτ
其中,所述Q值插值计算公式为:
Q τ j = Q τ ( j - 1 ) ( τ ( j + 1 ) - τ j ) + Q τ ( j + 1 ) ( τ j - τ ( j - 1 ) ) τ ( j + 1 ) - τ ( j - 1 ) ;
τ ( j - 1 ) = ( int ( τ j 0.02 ) - 1 ) × 0.02 , τ ( j + 1 ) = ( int ( τ j 0.02 ) + 1 ) × 0.02 , 所述τj、所述τ(j-1)和所述τ(j+1)都为时间深度,所述为时间深度τj处的Q值插值,所述为时间深度τ(j-1)处的Q值插值,所述为时间深度τ(j+1)处的Q值插值,表示对所述取整。
进一步地,所述根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到所述D个处理后的地震数据块,包括:
根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道;
根据所述每个地震数据块中的所有处理后的地震数据道,得到处理后的地震数据块;
根据所述D个地震数据块中的所有处理后的地震数据块,得到所述D个处理后的地震数据块。
进一步地,所述根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道,包括:
分别确定所述GPU中的线程的线程格尺寸和线程块尺寸;
对所述D个地震数据块中的每个地震数据块中的每个地震数据道进行快速傅里叶变换;
根据快速傅里叶变换结果,确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
根据所述角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz
根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
根据所述GPU中的所有线程处理得到的所述每个地震数据道对应的所有时间深度处的所有频率点对应的地震幅值,处理得到所述处理后的地震数据道。
进一步地,所述GPU中的线程的线程格尺寸为:所述Tracemax为所述GPU在一个周期内处理的地震数据的道数,b_size=2m,所述m为正整数,且4≤m≤6,所述S为地震数据道的采样点数;
所述GPU中的线程的线程块尺寸为:(b_size,b_size),b_size=2m,所述m为正整数,且4≤m≤6;
所述根据快速傅里叶变换结果,确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω,包括:
根据快速傅里叶变换结果,采用角频率采样间隔计算公式确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
其中,所述角频率采样间隔计算公式为:△ω=2π/(Gdt),所述G为满足快速傅里叶变换算法的时间采样间隔,所述dt为地震数据的时间采样间隔;
所述根据所述角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,包括:
根据所述目标低角频率值ωmin和所述角频率采样间隔△ω,采用起算角频率计算公式计算所述每个地震数据道对应的起算角频率值ωp
根据所述目标高角频率值ωmax和所述角频率采样间隔△ω,采用终止角频率计算公式计算所述每个地震数据道对应的终止角频率值ωz
其中,所述起算角频率计算公式为:所述终止角频率计算公式为:所述△f为频率增量,所述ωmax+2π△f为通过反Q滤波处理后期望达到的角频率的上限;
所述根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,包括:
根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,采用地震幅值计算公式计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
其中,所述地震幅值计算公式为: 所述A为所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且所述A为复数,所述k为所述每个频率点的序号,所述△ω为所述角频率采样间隔,所述τj为时间深度,所述为时间深度τj处的Q值;
所述根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A,包括:
根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,采用幅值和计算公式计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
其中,所述幅值和计算公式为:
Σ A = Σ k = int ( ω min Δ ω ) int ( ω max + 2 π Δ f Δ ω ) Re { F ( k Δ ω ) exp ( kΔωτ j 2 Q τ j ) } , 所述A为所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且所述A为复数,所述k为所述每个频率点的序号,所述△ω为所述角频率采样间隔,所述τj为时间深度,所述为时间深度τj处的Q值,所述Re表示取实部。
第三方面,提供一种叠后地震数据体处理装置,用于中央处理器CPU,所述装置包括:
选取模块,用于从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线,所述M为大于或者等于2的整数,且M=0.003N,所述N为所述叠后地震数据体的地震主测线的条数,每条所述地震测线对应一个地震剖面;
第一确定模块,用于从所述地震工区的叠加速度场中确定所述M条样本测线中的每条样本测线对应的速度场,得到所述M条样本测线对应的样本速度场ν;
第二确定模块,用于根据所述M条样本测线对应的样本速度场ν,采用李庆忠经验公式确定所述M条样本测线对应的初始Q值场Q0,所述李庆忠经验公式为Q=α·νβ,α=14,β=2.2;
校正模块,用于对所述M条样本测线对应的初始Q值场Q0进行校正,得到所述M条样本测线对应的校正Q值场Q1
第三确定模块,用于根据所述M条样本测线对应的校正Q值场Q1和所述M条样本测线对应的加权系数场λfinal,确定所述M条样本测线对应的最终Q值场Q2
更新模块,用于根据所述M条样本测线对应的样本速度场ν和所述M条样本测线对应的最终Q值场Q2,采用最小二乘法更新所述李庆忠经验公式,得到更新后的李庆忠经验公式;
第四确定模块,用于根据所述地震工区的叠加速度场,采用所述更新后的李庆忠经验公式确定所述地震工区的Q值场Qall
发送模块,用于向图形处理器GPU发送所述地震工区的Q值场Qall和所述叠后地震数据体,以便于所述GPU根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
进一步地,所述装置还包括:
第一处理模块,用于根据预设的加权系数组对所述M条样本测线对应的校正Q值场Q1进行处理,得到所述M条样本测线对应的加权Q值场组,所述预设的加权系数组中包括P个加权系数λ,所述加权Q值场组中包括P套加权Q值场,所述P为大于或者等于1的整数;
第二处理模块,用于根据所述加权Q值场组对所述M条样本测线组成的样本地震数据体进行反Q滤波处理,得到反Q滤波数据体组,所述反Q滤波数据体组中包括P套反Q滤波数据体,每套反Q滤波数据体包括M条反Q样本测线,每条所述反Q样本测线对应一个地震剖面;
频谱分析模块,用于采用傅里叶变换对所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析;
第五确定模块,用于根据频谱分析结果进行处理,确定所述M条反Q样本测线对应的加权系数场;
第六确定模块,用于将所述M条反Q样本测线对应的加权系数场确定为所述M条样本测线对应的加权系数场λfinal
进一步地,所述频谱分析模块,具体用于:
在所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面上选择三个矩形时窗,所述三个矩形时窗包括:浅层矩形时窗、中层矩形时窗和深层矩形时窗;
采用傅里叶变换对所述每条反Q样本测线对应的地震剖面上位于所述三个矩形时窗中的每个矩形时窗内的部分进行处理,得到所述每条反Q样本测线对应的3个反Q样本频谱;
根据所述每套反Q滤波数据体中的M条反Q样本测线对应的反Q样本频谱,得到所述每套反Q滤波数据体对应的3M个反Q样本频谱;
根据所述反Q滤波数据体组中的P套反Q滤波数据体对应的反Q样本频谱得到3×M×P个反Q样本频谱,所述3×M×P个反Q样本频谱中包括:M×P个浅层矩形时窗对应的反Q样本频谱、M×P个中层矩形时窗对应的反Q样本频谱和M×P个深层矩形时窗对应的反Q样本频谱;
对所述3×M×P个反Q样本频谱中的每个反Q样本频谱进行频谱分析。
进一步地,所述第五确定模块,具体用于:
根据频谱分析结果确定所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱,所述第一目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述浅层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,所述第二目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述中层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,所述第三目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述深层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱;
分别确定所述第一目标频谱对应的矩形时窗对应的第一目标加权系数、所述第二目标频谱对应的矩形时窗对应的第二目标加权系数和所述第三目标频谱对应的矩形时窗对应的第三目标加权系数,得到所述每条反Q样本测线对应的3个目标加权系数,其中,所述第一目标加权系数、所述第二目标加权系数和所述第三目标加权系数都为所述预设的加权系数组中的加权系数;
根据所述M条反Q样本测线中的每条反Q样本测线对应的3个目标加权系数,得到所述M条反Q样本测线对应的3M个目标加权系数;
对所述3M个目标加权系数进行空间插值,得到所述地震工区对应的加权系数场;
对所述地震工区对应的加权系数场进行平滑处理,得到处理后的加权系数场;
从所述处理后的加权系数场中抽取出所述M条反Q样本测线对应的加权系数场。
进一步地,所述装置还包括:
采样模块,用于按照预设采样间隔对所述地震工区的Q值场Qall进行采样,得到采样Q值场Qs
分块模块,用于对所述叠后地震数据体进行分块处理,得到D个地震数据块,所述D为大于或者等于2的整数;
所述发送模块,具体用于:
依次向所述GPU发送所述采样Q值场Qs和所述D个地震数据块中的每个地震数据块。
进一步地,所述分块模块,具体用于:
计算所述叠后地震数据体中的地震数据的总道数E;
根据所述GPU在一个周期内处理的地震数据的道数Tracemax和所述叠后地震数据体中的地震数据的总道数E,对所述叠后地震数据体进行分块处理,得到所述D个地震数据块;
其中,E=Tracemax×(D-1)+L,所述L为所述D个地震数据块中的最后一个地震数据块中的地震数据的道数。
第四方面,提供一种叠后地震数据体处理装置,用于图形处理器GPU,所述装置包括:
接收模块,用于接收中央处理器CPU发送的地震工区的Q值场Qall和叠后地震数据体;
处理模块,用于根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
进一步地,所述接收模块,具体用于:
接收所述CPU发送的采样Q值场Qs和D个地震数据块中的每个地震数据块,所述采样Q值场Qs是所述CPU按照预设采样间隔对所述地震工区的Q值场Qall进行采样得到的,所述D个地震数据块是所述CPU对所述叠后地震数据体进行分块处理得到的;
所述处理模块,包括:
处理子模块,用于根据所述采样Q值场Qs对所述D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块;
得到子模块,用于根据所述D个处理后的地震数据块,得到所述处理后的叠后地震数据体。
进一步地,所述处理子模块,包括:
第一处理单元,用于对所述采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ
第二处理单元,用于根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到所述D个处理后的地震数据块。
进一步地,所述第一处理单元,具体用于:
根据所述采样Q值场Qs和Q值插值计算公式确定每个时间深度τ处的Q值Qτ
其中,所述Q值插值计算公式为:
Q τ j = Q τ ( j - 1 ) ( τ ( j + 1 ) - τ j ) + Q τ ( j + 1 ) ( τ j - τ ( j - 1 ) ) τ ( j + 1 ) - τ ( j - 1 ) ;
τ ( j - 1 ) = ( int ( τ j 0.02 ) - 1 ) × 0.02 , τ ( j + 1 ) = ( int ( τ j 0.02 ) + 1 ) × 0.02 , 所述τj、所述τ(j-1)和所述τ(j+1)都为时间深度,所述为时间深度τj处的Q值插值,所述为时间深度τ(j-1)处的Q值插值,所述为时间深度τ(j+1)处的Q值插值,表示对所述取整。
进一步地,所述第二处理单元,包括:
处理子单元,用于根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道;
第一得到子单元,用于根据所述每个地震数据块中的所有处理后的地震数据道,得到处理后的地震数据块;
第二得到子单元,用于根据所述D个地震数据块中的所有处理后的地震数据块,得到所述D个处理后的地震数据块。
进一步地,所述处理子单元,包括:
第一确定子单元,用于分别确定所述GPU中的线程的线程格尺寸和线程块尺寸;
变换子单元,用于对所述D个地震数据块中的每个地震数据块中的每个地震数据道进行快速傅里叶变换;
第二确定子单元,用于根据快速傅里叶变换结果,确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
第三确定子单元,用于根据所述角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz
第一计算子单元,用于根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
第二计算子单元,用于根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
第一处理子单元,用于根据所述GPU中的所有线程处理得到的所述每个地震数据道对应的所有时间深度处的所有频率点对应的地震幅值,处理得到所述处理后的地震数据道。
进一步地,所述GPU中的线程的线程格尺寸为:所述Tracemax为所述GPU在一个周期内处理的地震数据的道数,b_size=2m,所述m为正整数,且4≤m≤6,所述S为地震数据道的采样点数;
所述GPU中的线程的线程块尺寸为:(b_size,b_size),b_size=2m,所述m为正整数,且4≤m≤6;
所述第二确定子单元,具体用于:
根据快速傅里叶变换结果,采用角频率采样间隔计算公式确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
其中,所述角频率采样间隔计算公式为:△ω=2π/(Gdt),所述G为满足快速傅里叶变换算法的时间采样间隔,所述dt为地震数据的时间采样间隔;
所述第三确定子单元,具体用于:
根据所述目标低角频率值ωmin和所述角频率采样间隔△ω,采用起算角频率计算公式计算所述每个地震数据道对应的起算角频率值ωp
根据所述目标高角频率值ωmax和所述角频率采样间隔△ω,采用终止角频率计算公式计算所述每个地震数据道对应的终止角频率值ωz
其中,所述起算角频率计算公式为:所述终止角频率计算公式为:所述△f为频率增量,所述ωmax+2π△f为通过反Q滤波处理后期望达到的角频率的上限;
所述第一计算子单元,具体用于:
根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,采用地震幅值计算公式计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
其中,所述地震幅值计算公式为: 所述A为所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且所述A为复数,所述k为所述每个频率点的序号,所述△ω为所述角频率采样间隔,所述τj为时间深度,所述为时间深度τj处的Q值;
所述第二计算子单元,具体用于:
根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,采用幅值和计算公式计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
其中,所述幅值和计算公式为:
Σ A = Σ k = int ( ω min Δ ω ) int ( ω max + 2 π Δ f Δ ω ) Re { F ( k Δ ω ) exp ( kΔωτ j 2 Q τ j ) } , 所述A为所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且所述A为复数,所述k为所述每个频率点的序号,所述△ω为所述角频率采样间隔,所述τj为时间深度,所述为时间深度τj处的Q值,所述Re表示取实部。
本发明提供的技术方案带来的有益效果是:
本发明提供的叠后地震数据体处理方法及装置,方法包括:从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线;从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场;根据M条样本测线对应的样本速度场,采用李庆忠经验公式确定M条样本测线对应的初始Q值场;对M条样本测线对应的初始Q值场进行校正,得到M条样本测线对应的校正Q值场;根据M条样本测线对应的校正Q值场和M条样本测线对应的加权系数场,确定M条样本测线对应的最终Q值场;根据M条样本测线对应的样本速度场和M条样本测线对应的最终Q值场,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式;根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场;向GPU发送地震工区的Q值场和叠后地震数据体,以便于GPU根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。由于针对叠后地震数据体对李庆忠经验公式进行了更新,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题。
进一步地,通过CPU对叠后地震数据体进行分块处理得到地震数据块,使得GPU根据地震数据块对叠后地震数据体进行处理,由于GPU的处理能力较强,本发明通过采用GPU对叠后地震数据体进行反Q滤波处理,提高了对叠后地震数据体的处理速度和处理效率。
应当理解的是,以上的一般描述和后文的细节描述仅是示例性的,并不能限制本发明。
附图说明
为了更清楚地说明本发明实施例中的技术方案,下面将对实施例描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1是本发明实施例提供的一种叠后地震数据体处理方法的方法流程图;
图2是本发明实施例提供的另一种叠后地震数据体处理方法的方法流程图;
图3-1是本发明实施例提供的再一种叠后地震数据体处理方法的方法流程图;
图3-2是图3-1所示实施例提供的一种采用傅里叶变换对每条反Q样本测线对应的地震剖面进行频谱分析的方法流程图;
图3-3是图3-1所示实施例提供的一种根据频谱分析结果进行处理确定M条反Q样本测线对应的加权系数场的方法流程图;
图3-4是图3-1所示实施例提供的一种对叠后地震数据体进行分块处理的方法流程图;
图3-5是图3-1所示实施例提供的一种根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理的方法流程图;
图3-6是图3-1所示实施例提供的一种根据采样Q值场对D个地震数据块中的每个地震数据块进行反Q滤波处理的方法流程图;
图3-7是图3-1所示实施例提供的一种根据所有时间深度处的Q值对每个地震数据块中的所有时间深度处的地震幅值进行反Q滤波处理的方法流程图;
图3-8是图3-1所示实施例提供的一种根据所有时间深度处的Q值对每个地震数据道的所有时间深度处的地震幅值进行反Q滤波处理的方法流程图;
图3-9是本发明实施例提供的一种处理前的地震剖面图;
图3-10是图3-9所示的地震剖面图上的部分区域的放大图;
图3-11是本发明实施例提供的对图3-9所示的地震剖面图进行处理得到的地震剖面图;
图3-12是图3-11所示的地震剖面图上的部分区域的放大图;
图3-13是本发明实施例提供的一种处理前的水平时间切片图;
图3-14是是本发明实施例提供的对图3-13所示的水平时间切片图进行处理得到的水平时间切片图;
图4是本发明实施例提供的一种叠后地震数据体处理装置的框图;
图5是本发明实施例提供的另一种叠后地震数据体处理装置的框图;
图6是本发明实施例提供的再一种叠后地震数据体处理装置的框图;
图7-1是本发明实施例提供的又一种叠后地震数据体处理装置的框图;
图7-2是图7-1所示实施例提供的一种处理模块的框图;
图7-3是图7-1所示实施例提供的一种处理子模块的框图;
图7-4是图7-1所示实施例提供的一种第二处理单元的框图;
图7-5是图7-1所示实施例提供的一种处理子单元的框图。
此处的附图被并入说明书中并构成本说明书的一部分,示出了符合本发明的实施例,并与说明书一起用于解释本发明的原理。
具体实施方式
为了使本发明的目的、技术方案和优点更加清楚,下面将结合附图对本发明作进一步地详细描述,显然,所描述的实施例仅仅是本发明一部份实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例,都属于本发明保护的范围。
请参考图1,其示出了本发明实施例提供的一种叠后地震数据体处理方法的方法流程图,该叠后地震数据体处理方法可以用于中央处理器(英文:Central ProcessingUnit;简称:CPU),参见图1,该方法流程具体包括:
在步骤101中,从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线。
其中,M为大于或者等于2的整数,且M=0.003N,N为叠后地震数据体的地震主测线的条数,每条地震测线对应一个地震剖面。
在步骤102中,从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场ν。
在步骤103中,根据M条样本测线对应的样本速度场ν,采用李庆忠经验公式确定M条样本测线对应的初始Q值场Q0
其中,李庆忠经验公式为Q=α·νβ,α=14,β=2.2。
在步骤104中,对M条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1
在步骤105中,根据M条样本测线对应的校正Q值场Q1和M条样本测线对应的加权系数场λfinal,确定M条样本测线对应的最终Q值场Q2
在步骤106中,根据M条样本测线对应的样本速度场ν和M条样本测线对应的最终Q值场Q2,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式。
在步骤107中,根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场Qall
在步骤108中,向图形处理器GPU发送地震工区的Q值场Qall和叠后地震数据体,以便于GPU根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
综上所述,本发明实施例提供的叠后地震数据体处理方法,通过从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线;从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场;根据M条样本测线对应的样本速度场,采用李庆忠经验公式确定M条样本测线对应的初始Q值场;对M条样本测线对应的初始Q值场进行校正,得到M条样本测线对应的校正Q值场;根据M条样本测线对应的校正Q值场和M条样本测线对应的加权系数场,确定M条样本测线对应的最终Q值场;根据M条样本测线对应的样本速度场和M条样本测线对应的最终Q值场,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式;根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场;向GPU发送地震工区的Q值场和叠后地震数据体,以便于GPU根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。由于针对叠后地震数据体对李庆忠经验公式进行了更新,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
进一步地,在步骤105之前,该方法还包括:
根据预设的加权系数组对M条样本测线对应的校正Q值场Q1进行处理,得到M条样本测线对应的加权Q值场组,预设的加权系数组中包括P个加权系数λ,加权Q值场组中包括P套加权Q值场,P为大于或者等于1的整数;
根据加权Q值场组对M条样本测线组成的样本地震数据体进行反Q滤波处理,得到反Q滤波数据体组,反Q滤波数据体组中包括P套反Q滤波数据体,每套反Q滤波数据体包括M条反Q样本测线,每条反Q样本测线对应一个地震剖面;
采用傅里叶变换对反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析;
根据频谱分析结果进行处理,确定M条反Q样本测线对应的加权系数场;
将M条反Q样本测线对应的加权系数场确定为M条样本测线对应的加权系数场λfinal
进一步地,采用傅里叶变换对反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析,包括:
在反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面上选择三个矩形时窗,三个矩形时窗包括:浅层矩形时窗、中层矩形时窗和深层矩形时窗;
采用傅里叶变换对每条反Q样本测线对应的地震剖面上位于三个矩形时窗中的每个矩形时窗内的部分进行处理,得到每条反Q样本测线对应的3个反Q样本频谱;
根据每套反Q滤波数据体中的M条反Q样本测线对应的反Q样本频谱,得到每套反Q滤波数据体对应的3M个反Q样本频谱;
根据反Q滤波数据体组中的P套反Q滤波数据体对应的反Q样本频谱得到3×M×P个反Q样本频谱,3×M×P个反Q样本频谱中包括:M×P个浅层矩形时窗对应的反Q样本频谱、M×P个中层矩形时窗对应的反Q样本频谱和M×P个深层矩形时窗对应的反Q样本频谱;
对3×M×P个反Q样本频谱中的每个反Q样本频谱进行频谱分析。
进一步地,根据频谱分析结果进行处理,确定M条反Q样本测线对应的加权系数场,包括:
根据频谱分析结果确定反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱,第一目标频谱为每条反Q样本测线对应的地震剖面上位于浅层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,第二目标频谱为每条反Q样本测线对应的地震剖面上位于中层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,第三目标频谱为每条反Q样本测线对应的地震剖面上位于深层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱;
分别确定第一目标频谱对应的矩形时窗对应的第一目标加权系数、第二目标频谱对应的矩形时窗对应的第二目标加权系数和第三目标频谱对应的矩形时窗对应的第三目标加权系数,得到每条反Q样本测线对应的3个目标加权系数,其中,第一目标加权系数、第二目标加权系数和第三目标加权系数都为预设的加权系数组中的加权系数;
根据M条反Q样本测线中的每条反Q样本测线对应的3个目标加权系数,得到M条反Q样本测线对应的3M个目标加权系数;
对3M个目标加权系数进行空间插值,得到地震工区对应的加权系数场;
对地震工区对应的加权系数场进行平滑处理,得到处理后的加权系数场;
从处理后的加权系数场中抽取出M条反Q样本测线对应的加权系数场。
进一步地,在步骤108之前,该方法还包括:
按照预设采样间隔对地震工区的Q值场Qall进行采样,得到采样Q值场Qs
对叠后地震数据体进行分块处理,得到D个地震数据块,D为大于或者等于2的整数;
向图形处理器GPU发送地震工区的Q值场Qall和叠后地震数据体,包括:
依次向GPU发送采样Q值场Qs和D个地震数据块中的每个地震数据块。
进一步地,对叠后地震数据体进行分块处理,得到D个地震数据块,包括:
计算叠后地震数据体中的地震数据的总道数E;
根据GPU在一个周期内处理的地震数据的道数Tracemax和叠后地震数据体中的地震数据的总道数E,对叠后地震数据体进行分块处理,得到D个地震数据块;
其中,E=Tracemax×(D-1)+L,L为D个地震数据块中的最后一个地震数据块中的地震数据的道数。
进一步地,步骤104可以包括:
采用谱比法确定M条样本测线对应的层Q值场Qi,i为垂直地震剖面VSP测井中的地层序号;
根据M条样本测线对应的层Q值场Qi,确定M条样本测线对应的时间深度的叠加Q值场Qvsp
根据M条样本测线对应的时间深度的叠加Q值场Qvsp,对M条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1
进一步地,步骤106可以包括:
根据M条样本测线对应的样本速度场ν和M条样本测线对应的最终Q值场Q2,采用最小二乘法确定李庆忠经验公式Q=α·νβ中的α值与β值,得到更新α值和更新β值;
将更新α值和更新β值代入李庆忠经验公式,得到更新后的李庆忠经验公式。
进一步地,根据M条样本测线对应的层Q值场Qi,确定M条样本测线对应的时间深度的叠加Q值场Qvsp,包括:
根据M条样本测线对应的层Q值场Qi,利用叠加Q值场公式确定M条样本测线对应的时间深度的叠加Q值场Qvsp
叠加Q值场公式为:
其中,i为垂直地震剖面VSP测井中的地层序号,△Ti为M条样本测线对应的地震波在第i层地层中的单程旅行时间,n为地层总数,T为M条样本测线对应的地震波从地表垂直传播至第i层地层的时间。
进一步地,在反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面上选择三个矩形时窗,包括:
在反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面上选择三个矩形时窗,使每个矩形时窗的大小相等,且使每个矩形时窗内的同相轴的清晰度大于预设清晰度,同相轴的连续度大于预设连续度。
进一步地,在步骤108之前,该方法还包括:
对叠后地震数据体进行去噪处理。
进一步地,步骤101可以包括:
从地震工区的叠后地震数据体对应的地震测线中选取M条存在垂直地震剖面VSP测井的地震测线作为M条样本测线;
其中,M条样本测线包括:地震主测线和联络测线中的至少一种。
进一步地,从地震工区的叠后地震数据体对应的地震测线中选取M条存在垂直地震剖面VSP测井的地震测线作为M条样本测线,包括:
根据大地坐标,将地震工区的所有VSP测井投射到工区平面图上;
根据工区平面图,从地震工区的叠后地震数据体对应的地震测线中选取M条存在垂直地震剖面VSP测井的地震测线作为M条样本测线。
进一步地,采用谱比法确定M条样本测线对应的层Q值场Qi,i为垂直地震剖面VSP测井中的地层序号,包括:
根据层Q值场公式确定M条样本测线对应的层Q值场Qi
层Q值场为: l n | a i ( f ) a i - 1 ( f ) | = c o n s t + π δ t Q i f ;
其中,f为M条样本测线中地层的频谱中的频率值,ai(f)为M条样本测线中第i层地层的频谱中频率f对应的地震幅值,ai-1(f)为M条样本测线中第i-1层地层的频谱中频率f对应的地震幅值,const表示常数,δt为地震波到达第i层地层的时间与地震波到达第i-1层地层的时间的差值。
进一步地,在步骤108之前,该方法还可以包括:
采用傅里叶变换对M条样本测线中的每条样本测线对应的地震剖面进行频谱分析,得到M条样本测线对应的目标低角频率值ωmin和目标高角频率值ωmax
该方法还可以包括:
向GPU发送目标低角频率值ωmin和目标高角频率值ωmax
进一步地,采用傅里叶变换对M条样本测线中的每条样本测线对应的地震剖面进行频谱分析,得到M条样本测线对应的目标低角频率值ωmin和目标高角频率值ωmax,包括:
在M条样本测线中的每条样本测线对应的地震剖面上选择三个矩形时窗,三个矩形时窗包括:浅层矩形时窗、中层矩形时窗和深层矩形时窗;
采用傅里叶变换对每条样本测线对应的地震剖面上位于三个矩形时窗中的每个矩形时窗内的部分进行处理,得到每条样本测线对应的3个样本频谱;
根据M条样本测线中的每条样本测线对应的3个样本频谱,得到M条样本测线对应的3M个样本频谱;
确定3M个样本频谱中的每个样本频谱对应的最高角频率值和最低角频率值,得到3M个最高角频率值和3M个最低角频率值;
将3M个最低角频率值中最小的角频率值作为M条样本测线对应的目标低角频率值ωmin
将3M个最高角频率值中最大的角频率值作为M条样本测线对应的目标高角频率值ωmax
综上所述,本发明实施例提供的叠后地震数据体处理方法,通过从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线;从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场;根据M条样本测线对应的样本速度场,采用李庆忠经验公式确定M条样本测线对应的初始Q值场;对M条样本测线对应的初始Q值场进行校正,得到M条样本测线对应的校正Q值场;根据M条样本测线对应的校正Q值场和M条样本测线对应的加权系数场,确定M条样本测线对应的最终Q值场;根据M条样本测线对应的样本速度场和M条样本测线对应的最终Q值场,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式;根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场;向GPU发送地震工区的Q值场和叠后地震数据体,以便于GPU根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。由于针对叠后地震数据体对李庆忠经验公式进行了更新,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
进一步地,通过CPU对叠后地震数据体进行分块处理得到地震数据块,使得GPU根据地震数据块对叠后地震数据体进行处理,由于GPU的处理能力较强,本发明实施例通过采用GPU对叠后地震数据体进行反Q滤波处理,提高了对叠后地震数据体的处理速度和处理效率。
请参考图2,其示出了本发明实施例提供的另一种叠后地震数据体处理方法的方法流程图,该叠后地震数据体处理方法可以用于图形处理器(英文:Graphic ProcessingUnit;简称:GPU),参见图2,该方法流程具体包括:
在步骤201中,接收中央处理器CPU发送的地震工区的Q值场Qall和叠后地震数据体。
在步骤202中,根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
综上所述,本发明实施例提供的叠后地震数据体处理方法,通过接收CPU发送的地震工区的Q值场和叠后地震数据体,根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体,由于地震工区的Q值场是CPU根据更新后的李庆忠经验公式计算得到的,且更新后的李庆忠经验公式是CPU针对叠后地震数据体对李庆忠经验公式进行更新得到的,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
进一步地,步骤201可以包括:
接收CPU发送的采样Q值场Qs和D个地震数据块中的每个地震数据块,采样Q值场Qs是CPU按照预设采样间隔对地震工区的Q值场Qall进行采样得到的,D个地震数据块是CPU对叠后地震数据体进行分块处理得到的;
步骤202可以包括:
根据采样Q值场Qs对D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块;
根据D个处理后的地震数据块,得到处理后的叠后地震数据体。
进一步地,根据采样Q值场Qs对D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块,包括:
对采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ
根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到D个处理后的地震数据块。
进一步地,对采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ,包括:
根据采样Q值场Qs和Q值插值计算公式确定每个时间深度τ处的Q值Qτ
其中,Q值插值计算公式为:
Q τ j = Q τ ( j - 1 ) ( τ ( j + 1 ) - τ j ) + Q τ ( j + 1 ) ( τ j - τ ( j - 1 ) ) τ ( j + 1 ) - τ ( j - 1 ) ;
τ ( j - 1 ) = ( int ( τ j 0.02 ) - 1 ) × 0.02 , τ ( j + 1 ) = ( int ( τ j 0.02 ) + 1 ) × 0.02 , τj、τ(j-1)和τ(j+1)都为时间深度,为时间深度τj处的Q值插值,为时间深度τ(j-1)处的Q值插值,为时间深度τ(j+1)处的Q值插值,表示对取整。
进一步地,根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到D个处理后的地震数据块,包括:
根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道;
根据每个地震数据块中的所有处理后的地震数据道,得到处理后的地震数据块;
根据D个地震数据块中的所有处理后的地震数据块,得到D个处理后的地震数据块。
进一步地,根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道,包括:
分别确定GPU中的线程的线程格尺寸和线程块尺寸;
对D个地震数据块中的每个地震数据块中的每个地震数据道进行快速傅里叶变换;
根据快速傅里叶变换结果,确定D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
根据角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定每个地震数据道对应的起算角频率值ωp和终止角频率值ωz
根据GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
根据每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
根据GPU中的所有线程处理得到的每个地震数据道对应的所有时间深度处的所有频率点对应的地震幅值,处理得到处理后的地震数据道。
进一步地,GPU中的线程的线程格尺寸为:Tracemax为GPU在一个周期内处理的地震数据的道数,b_size=2m,m为正整数,且4≤m≤6,S为地震数据道的采样点数;
GPU中的线程的线程块尺寸为:(b_size,b_size),b_size=2m,m为正整数,且4≤m≤6;
根据快速傅里叶变换结果,确定D个地震数据块中所有地震数据对应的角频率采样间隔△ω,包括:
根据快速傅里叶变换结果,采用角频率采样间隔计算公式确定D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
其中,角频率采样间隔计算公式为:△ω=2π/(Gdt),G为满足快速傅里叶变换算法的时间采样间隔,dt为地震数据的时间采样间隔;
根据角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,包括:
根据目标低角频率值ωmin和角频率采样间隔△ω,采用起算角频率计算公式计算每个地震数据道对应的起算角频率值ωp
根据目标高角频率值ωmax和角频率采样间隔△ω,采用终止角频率计算公式计算每个地震数据道对应的终止角频率值ωz
其中,起算角频率计算公式为:终止角频率计算公式为:△f为频率增量,ωmax+2π△f为通过反Q滤波处理后期望达到的角频率的上限;
根据GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,包括:
根据GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,采用地震幅值计算公式计算每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
其中,地震幅值计算公式为: A为每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且A为复数,k为每个频率点的序号,△ω为角频率采样间隔,τj为时间深度,为时间深度τj处的Q值;
根据每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A,包括:
根据每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,采用幅值和计算公式计算每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
其中,幅值和计算公式为: Σ A = Σ k = int ( ω min Δ ω ) int ( ω max + 2 π Δ f Δ ω ) Re { F ( k Δ ω ) exp ( kΔωτ j 2 Q τ j ) } , A为每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且A为复数,k为每个频率点的序号,△ω为角频率采样间隔,τj为时间深度,为时间深度τj处的Q值,Re表示取实部。
进一步地,在所述根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到所述D个处理后的地震数据块之前,所述方法还包括:
接收所述CPU发送的目标低角频率值ωmin和目标高角频率值ωmax
进一步地,在步骤202之后,该方法还可以包括:
对处理后的叠后地震数据体进行图形显示。
综上所述,本发明实施例提供的叠后地震数据体处理方法,通过接收CPU发送的地震工区的Q值场和叠后地震数据体,根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体,由于地震工区的Q值场是CPU根据更新后的李庆忠经验公式计算得到的,且更新后的李庆忠经验公式是CPU针对叠后地震数据体对李庆忠经验公式进行更新得到的,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
进一步地,通过CPU对叠后地震数据体进行分块处理得到地震数据块,GPU根据地震数据块对叠后地震数据体进行处理,由于GPU的处理能力较强,本发明实施例通过采用GPU对叠后地震数据体进行反Q滤波处理,提高了对叠后地震数据体的处理速度和处理效率。
请参考图3-1,其示出了本发明实施例提供的再一种叠后地震数据体处理方法的方法流程图,该叠后地震数据体处理方法可以用于CPU-GPU平台,参见图3-1,该方法流程具体包括:
在步骤301中,CPU从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线。
其中,M为大于或者等于2的整数,且M=0.003N,N为叠后地震数据体的地震主测线的条数,每条地震测线对应一个地震剖面。在本发明实施例中,当N大于667时,M=0.003N,当N小于667时,M等于2。
示例地,在选取样本测线时,可以从地震工区的叠后地震数据体对应的地震测线中选取存在垂直地震剖面(英文:Vertical Seismic Profiling;简称:VSP)测井的地震测线作为样本测线,该样本测线可以包括:地震主测线和联络测线中的至少一种,需要说明的是,该样本测线也可以是叠后地震数据体对应的地震测线中除地震主测线和联络测线中之外的其他测线,且M条样本测线中任意相邻的两条样本测线之间的距离可以相等,也可以不相等,本发明实施例对此不做限定。
可选地,在了便于样本测线的选择,可以先根据大地坐标,将地震工区的所有VSP测井投射到工区平面图上,然后根据工区平面图,从地震工区的叠后地震数据体对应的地震测线中选取M条存在垂直地震剖面VSP测井的地震测线作为M条样本测线。其中,将将地震工区的所有VSP测井投射到工区平面图上可以便于从宏观上了解地震工区,便于样本测线的选择。
在步骤302中,CPU从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场ν。
在本发明实施例中,每条地震测线可以对应一个速度场,该速度场中包括该速度场对应的地震测线对应的地震剖面上的各个点对应的传播速度,地震工区的叠加速度场中包括所有地震测线对应的速度场,CPU可以直接从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场ν,本发明实施例对此不做限定。
在步骤303中,CPU根据M条样本测线对应的样本速度场ν,采用李庆忠经验公式确定M条样本测线对应的初始Q值场Q0
其中,李庆忠经验公式为Q=α·νβ,α=14,β=2.2。
得到M条样本测线对应的样本速度场ν后,可以将该M条样本测线对应的样本速度场ν代入李庆忠经验公式Q=α·νβ,求出M条样本测线对应的初始Q值场Q0
在步骤304中,CPU对M条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1
得到M条样本测线对应的初始Q值场Q0后,CPU对M条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1
可选地,在对M条样本测线对应的初始Q值场Q0进行校正时,CPU可以先采用谱比法确定M条样本测线对应的层Q值场Qi,i为垂直地震剖面VSP测井中的地层序号,然后根据M条样本测线对应的层Q值场Qi,确定M条样本测线对应的时间深度的叠加Q值场Qvsp,进而根据M条样本测线对应的时间深度的叠加Q值场Qvsp,对M条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1
在本发明实施例中,CPU采用谱比法确定M条样本测线对应的层Q值场Qi可以包括:CPU根据层Q值场公式确定M条样本测线对应的层Q值场Qi,其中,层Q值场为:f为M条样本测线中地层的频谱中的频率值,ai(f)为M条样本测线中第i层地层的频谱中频率f对应的地震幅值,ai-1(f)为M条样本测线中第i-1层地层的频谱中频率f对应的地震幅值,const表示常数,δt为地震波到达第i层地层的时间与地震波到达第i-1层地层的时间的差值。其中,f、ai(f)、ai-1(f)和δt都可以预先获知,且CPU获知f、ai(f)、ai-1(f)和δt的过程可以参考现有技术,本发明实施例在此不再赘述。
在本发明实施例中,CPU然后根据M条样本测线对应的层Q值场Qi,确定M条样本测线对应的时间深度的叠加Q值场Qvsp可以包括:CPU根据M条样本测线对应的层Q值场Qi,利用叠加Q值场公式确定M条样本测线对应的时间深度的叠加Q值场Qvsp,其中,叠加Q值场公式为: i为垂直地震剖面VSP测井中的地层序号,△Ti为M条样本测线对应的地震波在第i层地层中的单程旅行时间,n为地层总数,T为M条样本测线对应的地震波从地表垂直传播至第i层地层的时间。在本发明实施例中,i、△Ti、T和n的获取过程均可以参考现有技术,本发明实施例在此不再赘述。
在本发明实施例中,CPU根据M条样本测线对应的时间深度的叠加Q值场Qvsp,对M条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1的过程可以包括:CPU将M条样本测线中的每条样本测线上位于VSP井位置(井口大地坐标为(x,y))的各时间深度点的初始Q值与该时间深度点对应的Qvsp逐一作相除运算获得校正系数η(x,y,τ)=Qvsp(τ)/Q0(τ),其中,Q0(τ)是时间深度τ处的初始Q值,Qvsp(τ)是时间深度τ处由VSP井获得的Q值,η(x,y,τ)是大地坐标为(x,y)且时间深度为τ时的时间深度点对应的校正系数,利用校正系数η(x,y,τ)在其大地坐标所在的样本测线剖面内进行线性插值并做二维平滑处理,得到每条样本测线对应的校正系数场,进而根据M条样本测线中的所有样本测线对应的校正系数场得到M条样本测线对应的校正系数场,在M条样本测线对应的地震剖面上逐时间深度点将校正系数场中的校正系数与初始Q值场Q0中的Q值进行相乘运算,得到校正后的Q值场Q1
在步骤305中,CPU根据预设的加权系数组对M条样本测线对应的校正Q值场Q1进行处理,得到M条样本测线对应的加权Q值场组。
其中,预设的加权系数组中包括P个加权系数λ,加权Q值场组中包括P套加权Q值场,P为大于或者等于1的整数。示例地,在本发明实施例中,P等于5,加权系数组中的加权系数可以包括:0.8、0.9、1.0、1.1、1.2,本发明实施例对此不做限定。
在本发明实施例中,CPU根据预设的加权系数组对M条样本测线对应的校正Q值场Q1进行处理,得到M条样本测线对应的加权Q值场组可以包括:CPU将加权系数组中的每个加权系数都与M条样本测线对应的校正Q值场Q1相乘,得到M条样本测线对应的加权Q值场组。示例地,当加权系数可以包括:0.8、0.9、1.0、1.1、1.2时,加权Q值场组中的加权Q值场包括:0.8Q1、0.9Q1、1.0Q1、1.1Q1和1.2Q1
在步骤306中,CPU根据加权Q值场组对M条样本测线组成的样本地震数据体进行反Q滤波处理,得到反Q滤波数据体组。
得到M条样本测线对应的加权Q值场组后,CPU可以根据M条样本测线对应的加权Q值场组对M条样本测线组成的样本地震数据体进行反Q滤波处理,得到反Q滤波数据体组。其中,反Q滤波数据体组中包括P套反Q滤波数据体,每套反Q滤波数据体包括M条反Q样本测线,每条反Q样本测线对应一个地震剖面。
在本发明实施例中,CPU根据加权Q值场组对M条样本测线组成的样本地震数据体进行反Q滤波处理的过程可以包括:CPU采用公式对M条样本测线组成的样本地震数据体进行处理,其中,在该公式中,f(t)表示地震幅值,F(ω)表示地震数据的傅里叶变换结果,ω表示角频率,τ表示时间深度,Qq表示加权Q值场组中的Q值,表示傅里叶逆变换,本发明实施例在此不再赘述。
在步骤307中,CPU采用傅里叶变换对反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析。
CPU得到反Q滤波数据体组后,可以对反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析。示例地,请参考图3-2,其示例出的是图3-1所示实施例提供的一种采用傅里叶变换对每条反Q样本测线对应的地震剖面进行频谱分析的方法流程图,参见图3-2,该方法流程可以包括:
在子步骤3071中,在反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面上选择三个矩形时窗。
在本发明实施例中,反Q滤波数据体组中包括P套反Q滤波数据体,每套反Q滤波数据体包括M条反Q样本测线,每条反Q样本测线对应一个地震剖面。CPU可以在每条饭Q样本测线对应的地震剖面上选择三个矩形时窗,其中,三个矩形时窗包括:浅层矩形时窗、中层矩形时窗和深层矩形时窗。
需要说明的是,在本发明实施例中,CPU在选择矩形时窗时,选择的三个矩形时窗的大小相等,且三个矩形时窗无重叠区域,每个矩形时窗内的同相轴的清晰度大于预设清晰度,同相轴的连续度大于预设连续度。预设清晰度用于区分矩形时窗内的同相轴是否清晰,当矩形时窗内的同相轴的清晰度大于或者等于预设清晰度时,表示矩形时窗内的同相轴较为清晰,当矩形时窗内的同相轴的清晰度小于预设清晰度时,表示矩形时窗内的同相轴的不清晰;预设连续度用于分矩形时窗内的同相轴是否连续,当矩形时窗内的同相轴的连续度大于或者等于预设连续度时,表示矩形时窗内的同相轴连续,当矩形时窗内的同相轴的连续度小于预设连续度时,表示矩形时窗内的同相轴的不连续,其中,预设清晰度和预设连续度的具体数值均可以根据实际需要进行设置,本发明实施例对此不做限定。
在子步骤3072中,采用傅里叶变换对每条反Q样本测线对应的地震剖面上位于三个矩形时窗中的每个矩形时窗内的部分进行处理,得到每条反Q样本测线对应的3个反Q样本频谱。
在本发明实施例中,每条反Q样本测线对应的地震剖面上的三个矩形时窗中的每个矩形时窗中可以包括每条反Q样本测线对应的地震剖面的一部分,为了简化分析过程,CPU可以采用傅里叶变换对每条反Q样本测线对应的地震剖面上位于每个矩形时窗内的部分进行处理,得到每条反Q样本测线对应的3个反Q样本频谱。其中,CPU采用傅里叶变换对每条反Q样本测线对应的地震剖面上位于三个矩形时窗中的每个矩形时窗内的部分进行处理,得到每条反Q样本测线对应的3个反Q样本频谱的过程可以参考现有技术,本发明实施例在此不再赘述。
在子步骤3073中,根据每套反Q滤波数据体中的M条反Q样本测线对应的反Q样本频谱,得到每套反Q滤波数据体对应的3M个反Q样本频谱。
CPU采用傅里叶变换对每条反Q样本测线对应的地震剖面上位于三个矩形时窗中的每个矩形时窗内的部分进行处理,得到每条反Q样本测线对应的3个反Q样本频谱后,由于每套反Q滤波数据体中包括M条反Q样本测线,因此,CPU可以根据每套反Q滤波数据体中的M条反Q样本测线对应的反Q样本频谱,得到每套反Q滤波数据体对应的3M个反Q样本频谱。示例地,当M=2时,3M=6,也即是,每套反Q滤波数据体对应的6个反Q样本频谱。
在子步骤3074中,根据反Q滤波数据体组中的P套反Q滤波数据体对应的反Q样本频谱得到3×M×P个反Q样本频谱。
CPU得到每套反Q滤波数据体对应的3M个反Q样本频谱后,由于反Q滤波数据体组中包括P套反Q滤波数据体,因此,CPU可以根据反Q滤波数据体组中的P套反Q滤波数据体对应的反Q样本频谱,得到3×M×P个反Q样本频谱。其中,3×M×P个反Q样本频谱中包括:M×P个浅层矩形时窗对应的反Q样本频谱、M×P个中层矩形时窗对应的反Q样本频谱和M×P个深层矩形时窗对应的反Q样本频谱。示例地,当M=2,P=5时,3×M×P=30,也即是,CPU根据反Q滤波数据体组中的5套反Q滤波数据体对应的反Q样本频谱得到30个反Q样本频谱,该30个反Q样本频谱中包括:10个浅层矩形时窗对应的反Q样本频谱、10个中层矩形时窗对应的反Q样本频谱和10个深层矩形时窗对应的反Q样本频谱,本发明实施例在此不再赘述。
在子步骤3075中,对3×M×P个反Q样本频谱中的每个反Q样本频谱进行频谱分析。
CPU得到3×M×P个反Q样本频谱后,可以对该3×M×P个反Q样本频谱中的每个反Q样本频谱进行频谱分析,示例地,参见子步骤3074可知,3×M×P=30,因此,CPU对该30个反Q样本频谱中的每个反Q样本频谱进行频谱分析。其中,CPU对每个反Q样本频谱进行频谱分析的过程可以参考现有技术,本发明实施例在此不再赘述。
在步骤308中,CPU根据频谱分析结果进行处理,确定M条反Q样本测线对应的加权系数场。
CPU采用傅里叶变换对反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析后,可以根据频谱分析结果进行处理,确定M条反Q样本测线对应的加权系数场。示例地,请参考图3-3,其示例出的是图3-1所示实施例提供的一种根据频谱分析结果进行处理确定M条反Q样本测线对应的加权系数场的方法流程图,参见图3-3,该方法流程可以包括:
在子步骤3081中,根据频谱分析结果确定反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱。
在本发明实施例中,参见子步骤3074可知,每条反Q样本测线对应的地震剖面上位于浅层矩形时窗内的部分对应P个反Q样本频谱,位于中层矩形时窗内的部分对应P个反Q样本频谱,位于深层矩形时窗内的部分对应P个反Q样本频谱。CPU可以根据频谱分析结果确定反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱。其中,第一目标频谱为每条反Q样本测线对应的地震剖面上位于浅层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,第二目标频谱为每条反Q样本测线对应的地震剖面上位于中层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,第三目标频谱为每条反Q样本测线对应的地震剖面上位于深层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,该第一目标频谱、第二目标频谱和第三目标频谱都可以称为恰当的补偿频谱,本发明实施例对此不做限定。
在子步骤3082中,分别确定第一目标频谱对应的矩形时窗对应的第一目标加权系数、第二目标频谱对应的矩形时窗对应的第二目标加权系数和第三目标频谱对应的矩形时窗对应的第三目标加权系数,得到每条反Q样本测线对应的3个目标加权系数。
CPU确定每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱后,可以分别确定第一目标频谱对应的矩形时窗对应的第一目标加权系数、第二目标频谱对应的矩形时窗对应的第二目标加权系数和第三目标频谱对应的矩形时窗对应的第三目标加权系数,得到每条反Q样本测线对应的3个目标加权系数。其中,第一目标加权系数、第二目标加权系数和第三目标加权系数都为预设的加权系数组中的加权系数。
需要说明的是,由于第一目标频谱、第二目标频谱和第三目标频谱都是反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的反Q样本频谱,而反Q滤波数据体组是CPU根据M条样本测线对应的加权Q值场组对M条样本测线组成的样本地震数据体进行反Q滤波处理得到的,因此,CPU可以反推出第一目标频谱对应的矩形时窗对应的第一目标加权系数、第二目标频谱对应的矩形时窗对应的第二目标加权系数和第三目标频谱对应的矩形时窗对应的第三目标加权系数,本发明实施例在此不再赘述。
在子步骤3083中,根据M条反Q样本测线中的每条反Q样本测线对应的3个目标加权系数,得到M条反Q样本测线对应的3M个目标加权系数。
由于每条反Q样本测线对应的3个目标加权系数,因此,M条样本测线可以对应3M个目标加权系数,CPU可以直接根据M条反Q样本测线中的每条烦Q样本测线对应的3个目标加权系数得到3M个目标加权系数,本发明实施在此不再赘述。
在子步骤3084中,对3M个目标加权系数进行空间插值,得到地震工区对应的加权系数场。
得到M条反Q样本测线对应的3M个目标加权系数后,CPU可以对M条反Q样本测线对应的3M个目标加权系数进行空间插值,得到地震工区对应的加权系数场。
在本发明实施例中,CPU对3M个目标加权系数进行空间插值的过程可以包括:对于每条反Q样本测线而言,假设该每条反Q样本测线对应的地震剖面上浅层矩形时窗的中心点对应的大地坐标是(x1,y1),在时间深度τ处,该浅层矩形时窗的中心点对应的目标加权系数为λ1,假设该每条反Q样本测线对应的地震剖面上中层矩形时窗的中心点对应的大地坐标是(x2,y2),在时间深度τ处,该浅层矩形时窗的中心点对应的目标加权系数为λ2,假设该每条反Q样本测线对应的地震剖面上深层矩形时窗的中心点对应的大地坐标是(x3,y3),在时间深度τ处,该浅层矩形时窗的中心点对应的目标加权系数为λ3,可以得到3个空间加权系数点,也即,每条反Q样本测线对应3个空间加权系数点,因此,M条反Q样本测线可以对应3M个空间加权系数点,该3M个空间加权系数点都为已知的空间加权系数点,可以采用空间三维插值方法对该3M个已知的空间加权系数点进行空间插值,得到地震工区对应的加权系数场。其中,采用空间三维插值方法对该3M个已知的空间加权系数点进行空间插值的过程可以参考现有技术,本发明实施例在此不再赘述。
在子步骤3085中,对地震工区对应的加权系数场进行平滑处理,得到处理后的加权系数场。
由于在子步骤3084中,CPU对3M个目标加权系数进行空间插值得到地震工区对应的加权系数场是离散的加权系数场,该离散的加权系数场中包括的是多个离散的空间加权系数点,因此,CPU可以对该离散的加权系数场进行平滑处理,得到处理后的加权系数场。示例地,CPU对地震工区对应的加权系数场做空间三维平滑处理,得到处理后的加权系数场,本发明实施例在此不再赘述。
在子步骤3086中,从处理后的加权系数场中抽取出M条反Q样本测线对应的加权系数场。
CPU对地震工区对应的加权系数场进行平滑处理得到处理后的加权系数场后,可以从处理后的加权系数场中抽取出M条反Q样本测线对应的加权系数场,示例地,CPU根据M条样本测线对应的空间坐标值,从处理后的加权系数场中抽取出M条反Q样本测线对应的加权系数场。其中,CPU根据M条样本测线对应的空间坐标值从处理后的加权系数场中抽取出M条反Q样本测线对应的加权系数场的过程可以参考现有技术,本发明实施例在此不再赘述。
在步骤309中,CPU将M条反Q样本测线对应的加权系数场确定为M条样本测线对应的加权系数场λfinal
抽取出M条反Q样本测线对应的加权系数场后,CPU可以将M条反Q样本测线对应的加权系数场确定为M条样本测线对应的加权系数场λfinal,本发明实施例在此不再赘述。
在步骤310中,CPU根据M条样本测线对应的校正Q值场Q1和M条样本测线对应的加权系数场λfinal,确定M条样本测线对应的最终Q值场Q2
由于在步骤304中,CPU得到了M条样本测线对应的校正Q值场Q1,在步骤309中CPU得到了M条样本测线对应的加权系数场λfinal,因此,CPU可以根据步骤304中得到的M条样本测线对应的校正Q值场Q1和步骤309中得到的M条样本测线对应的加权系数场λfinal确定M条样本测线对应的最终Q值场Q2
示例地,在本发明实施例中,CPU可以将M条样本测线对应的校正Q值场Q1和M条样本测线对应的加权系数场λfinal相乘,将M条样本测线对应的校正Q值场Q1与M条样本测线对应的加权系数场λfinal的乘积作为M条样本测线对应的最终Q值场Q2,本发明实施例对此不做限定。
在步骤311中,CPU根据M条样本测线对应的样本速度场ν和M条样本测线对应的最终Q值场Q2,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式。
由于在步骤302中CPU得到了M条样本测线对应的样本速度场ν,在步骤310中CPU得到了M条样本测线对应的最终Q值场Q2,因此,CPU可以根据M条样本测线对应的样本速度场ν和M条样本测线对应的最终Q值场Q2采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式。
其中,CPU根据M条样本测线对应的样本速度场ν和M条样本测线对应的最终Q值场Q2采用最小二乘法更新李庆忠经验公式得到更新后的李庆忠经验公式可以包括:CPU根据M条样本测线对应的样本速度场ν和M条样本测线对应的最终Q值场Q2,采用最小二乘法确定李庆忠经验公式Q=α·νβ中的α值与β值,得到更新α值和更新β值,然后将更新α值和更新β值代入李庆忠经验公式,得到更新后的李庆忠经验公式。具体地,CPU将李庆忠经验公式Q=α·νβ中的α与β都看做未知数,将M条样本测线对应的样本速度场ν中的速度值和M条样本测线对应的最终Q值场Q2中的Q值代入李庆忠经验公式Q=α·νβ,计算得到更新α值和更新β值,然后将Q和ν看做未知数,将更新α值和更新β值代入李庆忠经验公式,得到更新后的李庆忠经验公式,本发明实施例在此不再赘述。
在步骤312中,CPU根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场Qall
得到更新后的李庆忠经验公式后,CPU可以根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场Qall,示例地,CPU将地震工区的叠加速度场代入更新后的李庆忠经验公式,得到地震工区的Q值场Qall,本发明实施例在此不再赘述。
在步骤313中,按照预设采样间隔对地震工区的Q值场Qall进行采样,得到采样Q值场Qs
CPU得到地震工区的Q值场Qall后,可以按照预设采样间隔对地震工区的Q值场Qall进行采样,得到采样Q值场Qs,其中,预设采样间隔可以根据实际需要进行设置,本发明实施例对此不做限定,示例地,预设采样间隔为20ms(中文:毫秒),CPU对地震工区的Q值场Qall进行采样的过程可以参考现有技术,本发明实施例在此不再赘述。
在步骤314中,对叠后地震数据体进行分块处理,得到D个地震数据块,D为大于或者等于2的整数。
由于叠后地震数据体中的数据量较大,因此,CPU可以对叠后地震数据体进行分块处理,得到D个地震数据块,其中,D为大于或者等于2的整数,且D的具体数值可以根据实际情况确定,本发明实施例对此不做限定。
示例地,请参考图3-4,其示例出的是图3-1所示实施例提供的一种对叠后地震数据体进行分块处理的方法流程图,参见图3-4,该方法流程可以包括:
在子步骤3141中,计算叠后地震数据体中的地震数据的总道数E。
其中,CPU计算叠后地震数据体中的地震数据的总道数E具体实现过程可以参考现有技术,本发明实施例在此不再赘述。
在子步骤3142中,根据GPU在一个周期内处理的地震数据的道数Tracemax和叠后地震数据体中的地震数据的总道数E,对叠后地震数据体进行分块处理,得到D个地震数据块。
得到叠后地震数据体中的地震数据的总道数E后,CPU可以根据GPU在一个周期内处理的地震数据的道数Tracemax和叠后地震数据体中的地震数据的总道数E,对叠后地震数据体进行分块处理,得到D个地震数据块,其中,E=Tracemax×(D-1)+L,L为D个地震数据块中的最后一个地震数据块中的地震数据的道数,其中,Tracemax的取值为2u,256≤2u≤0.001E,且u为正整数,本发明实施例对此不做限定。
在步骤315中,CPU向GPU发送地震工区的Q值场Qall和叠后地震数据体。
CPU得到地震工区的Q值场Qall后,可以向GPU发送地震工区的Q值场Qall和叠后地震数据体,以便于GPU根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
需要说明的是,由于在步骤314中CPU对叠后地震数据体进行分块处理得到了D个地震数据块,且在步骤313中CPU按照预设采样间隔对地震工区的Q值场Qall进行采样得到采样Q值场Qs,因此,CPU向GPU发送地震工区的Q值场Qall和叠后地震数据体也即是CPU依次向GPU发送采样Q值场Qs和D个地震数据块中的每个地震数据块,本发明实施例对此不作限定。
还需要说明的是,为了便于对叠后地震数据体进行处理,CPU在向GPU发送叠后地震数据体之前,可以对叠后地震数据体进行去噪处理,示例地,CPU对叠后地震数据体进行保幅去噪处理,本发明实施例对此不做限定。
还需要说明的是,在向GPU发送地震工区的Q值场Qall和叠后地震数据体之前,CPU可以采用傅里叶变换对M条样本测线中的每条样本测线对应的地震剖面进行频谱分析,得到M条样本测线对应的目标低角频率值ωmin和目标高角频率值ωmax,CPU还可以向GPU发送该目标低角频率值ωmin和目标高角频率值ωmax。其中,CPU采用傅里叶变换对M条样本测线中的每条样本测线对应的地震剖面进行频谱分析,得到M条样本测线对应的目标低角频率值ωmin和目标高角频率值ωmax的过程可以包括:CPU在叠后地震数据体对应的M条样本测线中的每条样本测线对应的地震剖面上选择三个矩形时窗,然后采用傅里叶变换对每条样本测线对应的地震剖面上位于三个矩形时窗中的每个矩形时窗内的部分进行处理,得到每条样本测线对应的3个样本频谱,根据M条样本测线中的每条样本测线对应的3个样本频谱,得到M条样本测线对应的3M个样本频谱,进而确定3M个样本频谱中的每个样本频谱对应的最高角频率值和最低角频率值,得到3M个最高角频率值和3M个最低角频率值,然后将3M个最低角频率值中最小的角频率值作为M条样本测线对应的目标低角频率值ωmin,将3M个最高角频率值中最大的角频率值作为M条样本测线对应的目标高角频率值ωmax。其中,该三个矩形时窗包括:浅层矩形时窗、中层矩形时窗和深层矩形时窗,在本发明实施例中,CPU在选择矩形时窗时,选择的三个矩形时窗的大小相等,且三个矩形时窗无重叠区域,每个矩形时窗内的同相轴的清晰度大于预设清晰度,同相轴的连续度大于预设连续度。预设清晰度用于区分矩形时窗内的同相轴是否清晰,当矩形时窗内的同相轴的清晰度大于或者等于预设清晰度时,表示矩形时窗内的同相轴较为清晰,当矩形时窗内的同相轴的清晰度小于预设清晰度时,表示矩形时窗内的同相轴的不清晰;预设连续度用于分矩形时窗内的同相轴是否连续,当矩形时窗内的同相轴的连续度大于或者等于预设连续度时,表示矩形时窗内的同相轴连续,当矩形时窗内的同相轴的连续度小于预设连续度时,表示矩形时窗内的同相轴的不连续,其中,预设清晰度和预设连续度的具体数值均可以根据实际需要进行设置,本发明实施例对此不做限定。
在步骤316中,GPU接收CPU发送的地震工区的Q值场Qall和叠后地震数据体。
CPU向GPU发送地震工区的Q值场Qall和叠后地震数据体时,GPU可以接收CPU发送的地震工区的Q值场Qall和叠后地震数据体。在本发明实施例中,由于CPU向GPU发送地震工区的Q值场Qall和叠后地震数据体可以包括CPU依次向GPU发送采样Q值场Qs和D个地震数据块中的每个地震数据块。因此,GPU接收CPU发送的地震工区的Q值场Qall和叠后地震数据体可以包括GPU接收CPU发送的采样Q值场Qs和D个地震数据块中的每个地震数据块,其中,采样Q值场Qs是CPU按照预设采样间隔对地震工区的Q值场Qall进行采样得到的,D个地震数据块是CPU对叠后地震数据体进行分块处理得到的。
在步骤317中,GPU根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
接收到CPU发送的地震工区的Q值场Qall和叠后地震数据体后,GPU可以根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。示例地,请参考图3-5,其示例出的是图3-1所示实施例提供的一种根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理的方法流程图,参见图3-5,该方法流程可以包括:
在子步骤3171中,根据采样Q值场Qs对D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块。
由于GPU接收到的实际上是地震数据块和采样Q值场Qs,因此,GPU根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体也即是GPU根据采样Q值场Qs对D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块,进而根据D个处理后的地震数据块处理得到处理后的叠后地震数据体。
示例地,请参考图3-6,其示例出的是图3-1所示实施例提供的一种根据采样Q值场对D个地震数据块中的每个地震数据块进行反Q滤波处理的方法流程图,参见图3-6,该方法流程可以包括:
在子步骤31711中,对采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ
由于采样Q值场Qs是CPU按照预设采样间隔对地震工区的Q值场Qall进行采样得到的,因此,采样Q值场Qs中的Q值较少,因此,GPU可以对采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ
其中,GPU对采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ可以包括:GPU根据采样Q值场Qs和Q值插值计算公式确定每个时间深度τ处的Q值Qτ。其中,Q值插值计算公式为:
Q τ j = Q τ ( j - 1 ) ( τ ( j + 1 ) - τ j ) + Q τ ( j + 1 ) ( τ j - τ ( j - 1 ) ) τ ( j + 1 ) - τ ( j - 1 ) ;
τ ( j - 1 ) = ( int ( τ j 0.02 ) - 1 ) × 0.02 , τ ( j + 1 ) = ( int ( τ j 0.02 ) + 1 ) × 0.02 , τj、τ(j-1)和τ(j+1)都为时间深度,为时间深度τj处的Q值插值,为时间深度τ(j-1)处的Q值插值,为时间深度τ(j+1)处的Q值插值,表示对取整。
在子步骤31712中,根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到D个处理后的地震数据块。
得到每个时间深度τ处的Q值Qτ后,GPU可以根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到D个处理后的地震数据块。示例地,请参考图3-7,其示例出的是图3-1所示实施例提供的一种根据所有时间深度处的Q值对D个地震数据块中的每个地震数据块中的所有时间深度处的地震幅值进行反Q滤波处理的方法流程图,参见图3-7,该方法流程可以包括:
在子步骤317121中,根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道。
在本发明实施例中,每个地震数据道包括多个时间深度τ,每个时间深度τ对应一个Q值Qτ和一个地震幅值,因此,GPU可以根据D个地震数据块中的每个地震数据块中的每个地震数据道的每个时间深度τ处的Q值Qτ对该每个时间深度τ处的地震幅值进行反Q滤波处理,进而根据每个地震数据道的所有时间深度τ处的Q值Qτ对该每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理。
示例地,请参考图3-8,其示例出的是图3-1所示实施例提供的一种根据所有时间深度处的Q值对D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度处的地震幅值进行反Q滤波处理的方法流程图,参见图3-8,该方法流程可以包括:
在子步骤3171211中,分别确定GPU中的线程的线程格尺寸和线程块尺寸。
GPU中可以包括多个线程,GPU可以确定GPU中的线程的线程格尺寸和线程块尺寸。其中,GPU中的线程的线程格尺寸为:Tracemax为GPU在一个周期内处理的地震数据的道数,b_size=2m,m为正整数,且4≤m≤6,S为地震数据道的采样点数,GPU中的线程的线程块尺寸为:(b_size,b_size),b_size=2m,m为正整数,且4≤m≤6,本发明实施例对此不做限定。
在子步骤3171212中,对D个地震数据块中的每个地震数据块中的每个地震数据道进行快速傅里叶变换。
其中,GPU对D个地震数据块中的每个地震数据块中的每个地震数据道进行快速傅里叶变换的过程可以参考现有技术,本发明实施例在此不再赘述。
在子步骤3171213中,根据快速傅里叶变换结果,确定D个地震数据块中所有地震数据对应的角频率采样间隔△ω。
对D个地震数据块中的每个地震数据块中的每个地震数据道进行快速傅里叶变换后,GPU可以根据快速傅里叶变换结果,确定D个地震数据块中所有地震数据对应的角频率采样间隔△ω。
其中,GPU根据快速傅里叶变换结果,确定D个地震数据块中所有地震数据对应的角频率采样间隔△ω可以包括:GPU根据快速傅里叶变换结果,采用角频率采样间隔计算公式确定D个地震数据块中所有地震数据对应的角频率采样间隔△ω。其中,角频率采样间隔计算公式为:△ω=2π/(Gdt),G为满足快速傅里叶变换算法的时间采样间隔,dt为地震数据的时间采样间隔,G和dt的求取过程可以参考现有技术,本发明实施例在此不再赘述。
在子步骤3171214中,根据角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定每个地震数据道对应的起算角频率值ωp和终止角频率值ωz
得到D个地震数据块中所有地震数据对应的角频率采样间隔△ω后,GPU可以根据角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,其中,GPU根据角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定每个地震数据道对应的起算角频率值ωp和终止角频率值ωz的过程可以包括:GPU根据目标低角频率值ωmin和角频率采样间隔△ω,采用起算角频率计算公式计算每个地震数据道对应的起算角频率值ωp,根据目标高角频率值ωmax和角频率采样间隔△ω,采用终止角频率计算公式计算每个地震数据道对应的终止角频率值ωz,其中,起算角频率计算公式为:终止角频率计算公式为:△f为频率增量,ωmax+2π△f为通过反Q滤波处理后期望达到的角频率的上限,△f和ωmax+2π△f的求取过程可以参考现有技术,本发明实施例在此不再赘述。
需要说明的是,GPU还可以计算xindex=blockDim.x*blockIdx.x+threadIdx.x,yindex=blockDim.y*blockIdx.y+threadIdx.y,其中xindex与yindex分别为线程在线程格x方向与y方向的索引值,blockDim.x与blockDim.y分别是线程在线程块x方向与y方向的尺寸值,blockIdx.x与blockIdx.y分别是线程块在线程格中x方向与y方向的索引值,threadIdx.x与threadIdx.y分别是线程在线程块中x方向与y方向的索引值,本发明实施例在此不再赘述。
在子步骤3171215中,根据GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A。
得到GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ、每个地震数据道对应的起算角频率值ωp和终止角频率值ωz后,GPU可以根据GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A。示例地,GPU可以根据GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,采用地震幅值计算公式计算每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;其中,地震幅值计算公式为: A为每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且A为复数,k为每个频率点的序号,△ω为角频率采样间隔,τj为时间深度,为时间深度τj处的Q值。
在子步骤3171216中,根据每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A。
得到每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A后,GPU可以根据每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A。示例地,GPU可以根据每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,采用幅值和计算公式计算每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
其中,幅值和计算公式为: Σ A = Σ k = int ( ω min Δ ω ) int ( ω max + 2 π Δ f Δ ω ) Re { F ( k Δ ω ) exp ( kΔωτ j 2 Q τ j ) } , A为每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且A为复数,k为每个频率点的序号,△ω为角频率采样间隔,τj为时间深度,为时间深度τj处的Q值,Re表示取实部。
在子步骤3171217中,根据GPU中的所有线程处理得到的每个地震数据道对应的所有时间深度处的所有频率点对应的地震幅值,处理得到处理后的地震数据道。
得到每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A后,GPU可以根据GPU中的所有线程处理得到的每个地震数据道对应的所有时间深度处的所有频率点对应的地震幅值,处理得到处理后的地震数据道,本发明实施例在此不再赘述。
在子步骤317122中,根据每个地震数据块中的所有处理后的地震数据道,得到处理后的地震数据块。
得到每个地震数据块中的处理后的地震数据道后,GPU可以根据每个地震数据块中的所有处理后的地震数据道,得到处理后的地震数据块。
在子步骤317123中,根据D个地震数据块中的所有处理后的地震数据块,得到D个处理后的地震数据块。
得到D个地震数据块中的每个处理后的地震数据块后,GPU可以根据D个地震数据块中的所有处理后的地震数据块,得到D个处理后的地震数据块。
在子步骤3172中,根据D个处理后的地震数据块,得到处理后的叠后地震数据体。
得到D个处理后的地震数据块后,GPU可以根据D个处理后的地震数据块,得到处理后的叠后地震数据体。
需要说明的是,在本发明实施例中,GPU在根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到D个处理后的地震数据块之前,可以先接收CPU发送的目标低角频率值ωmin和目标高角频率值ωmax;且GPU在根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体之后,还可以对处理后的叠后地震数据体进行图形显示,其中,GPU显示的图形可以指示地下构造的地层沉积样式与微小断裂,可以用于构造解释与砂体追踪,本发明实施例在此不再赘述。
综上所述,本发明实施例提供的叠后地震数据体处理方法,通过从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线;从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场;根据M条样本测线对应的样本速度场,采用李庆忠经验公式确定M条样本测线对应的初始Q值场;对M条样本测线对应的初始Q值场进行校正,得到M条样本测线对应的校正Q值场;根据M条样本测线对应的校正Q值场和M条样本测线对应的加权系数场,确定M条样本测线对应的最终Q值场;根据M条样本测线对应的样本速度场和M条样本测线对应的最终Q值场,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式;根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场;向GPU发送地震工区的Q值场和叠后地震数据体,以便于GPU根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。由于针对叠后地震数据体对李庆忠经验公式进行了更新,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
进一步地,现有技术中,通常在CPU中对叠后地震数据体进行处理,由于CPU在一个计算周期内只能对一个地震数据道进行处理,且其为单道输入输出,读写频繁,处理效率较低。本发明实施例提供的叠后地震数据体处理方法,通过CPU对叠后地震数据体进行分块处理得到地震数据块,GPU根据地震数据块对叠后地震数据体进行处理,由于GPU在一个计算周期内可以对多个地震数据道进行处理,GPU的处理能力较强,本发明实施例通过采用GPU对叠后地震数据体进行反Q滤波处理,提高了对叠后地震数据体的处理速度和处理效率。
本发明实施例提供的叠后地震数据体处理方法可以提高叠后地震数据的分辨率,特别适用于处理大规模的三维偏移叠加地震数据体,对提高叠后地震数据体的分辨率具有重要应用价值。
下面以一个具体的例子对本发明实施例提供的叠后地震数据体处理方法进行简单说明。具体地,以冀东油田柏各庄区块为例,CPU型号为Xeon E5-2687,GPU型号为QuadroK5000,具体的处理过程可以包括以下步骤:
(1)CPU从地震工区的叠后地震数据体对应的地震测线中选取2条样本测线。
具体包括:CPU统计地震工区内VSP测井数为3,根据大地坐标,将地震工区的所有VSP测井投射到工区平面图上,统计地震工区内地震主测线条数N=500,从500条地震主测线条中选取M=2条地震主测线作为样本测线,且使每条样本测线上分布一口VSP测井。
在选取样本测线之前,CPU对叠后地震数据体进行去噪处理,其中,叠后地震数据体的采样间隔为2ms,记录时长为2500ms,三维偏移叠后地震数据体中共包含300000个成像道(地震数据道),水平地震测线之间的道间距为15m(中文:米),垂直地震测线的道间距为40m,叠后地震数据体的覆盖面积是180平方千米,示例地,如图3-9所示,其示出的是垂直坐标为8000m处的沿地震测线方向切片的三维偏移叠加剖面。其中,在该图3-9中,横轴表示垂直距离LS,单位为km(中文:千米),纵轴表示时间深度τ,单位为s(中文:秒)。如图3-10所示,其示出的是图3-9所示的区域390部分的放大图。
(2)CPU从地震工区的叠加速度场中抽取出2条样本测线中的每条样本测线对应的速度场,得到该2条样本测线对应的样本速度场ν,该样本速度场ν内的速度值按照样本测线的测线号的顺序依次排列。
(3)CPU将样本速度场ν内的速度值代入李庆忠经验公式Q=α·νβ(α=14,β=2.2)得到该2条样本测线对应的初始的Q值场Q0
(4)CPU对该2条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1
具体包括:采用谱比法确定该2条样本测线对应的层Q值场Qi,i为垂直地震剖面VSP测井中的地层序号,然后将层Q值场Qi转换为该2条样本测线对应的时间深度的叠加Q值场Qvsp,进而根据该2条样本测线对应的时间深度的叠加Q值场Qvsp对该2条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1
(5)CPU根据预设的加权系数组对2条样本测线对应的校正Q值场Q1进行处理,得到2条样本测线对应的加权Q值场组。
具体包括:CPU将校正Q值场Q1分别与加权系数组中的加权系数λ进行相乘(λ取值为0.8,0.9,1.0,1.1,1.2),得到加权Q值场组,该加权Q值场组中包括5套加权Q值场。
(6)CPU根据加权Q值场组对2条样本测线组成的样本地震数据体进行反Q滤波处理,得到反Q滤波数据体组,该反Q滤波数据体组中包括5套反Q滤波数据体,每套反Q滤波数据体中包括2条反Q样本测线。
(7)CPU采用傅里叶变换对反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析。
(8)CPU根据频谱分析结果进行处理,确定2条反Q样本测线对应的加权系数场。
具体包括:,根据频谱分析结果确定反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱,分别确定第一目标频谱对应的矩形时窗对应的第一目标加权系数、第二目标频谱对应的矩形时窗对应的第二目标加权系数和第三目标频谱对应的矩形时窗对应的第三目标加权系数,得到每条反Q样本测线对应的3个目标加权系数,根据2条反Q样本测线中的每条反Q样本测线对应的3个目标加权系数,得到M条反Q样本测线对应的6个目标加权系数,对6个目标加权系数进行空间插值,得到地震工区对应的加权系数场,对地震工区对应的加权系数场进行平滑处理,得到处理后的加权系数场,从处理后的加权系数场中抽取出2条反Q样本测线对应的加权系数场。
(9)CPU将2条反Q样本测线对应的加权系数场确定为2条样本测线对应的加权系数场λfinal
(10)CPU将步骤(9)得到的加权系数场λfinal与步骤(4)得到的校正Q值场Q1进行相乘运算,得到2条样本测线对应的最终Q值场Q2
(11)CPU根据2条样本测线对应的样本速度场ν和2条样本测线对应的最终Q值场Q2,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式。
(12)CPU根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场Qall
(13)CPU按照预设采样间隔对地震工区的Q值场Qall进行采样,得到采样Q值场Qs
(14)CPU对叠后地震数据体进行分块处理,得到D个地震数据块。
具体包括:CPU计算叠后数据体的总道数为E=300000,设GPU在一个周期内处理的地震数据的道数Trace_max=256,将总道数为E=300000的叠后数据体按照每块含有Trace_max=256道进行分块操作,共产生D=1171块内含256道地震数据的数据块,会剩余不足256道的尾道部分,该尾道部分的道数为224,将该部分数据作为最后一块数据体。
(15)CPU将得到的地震工区的Q值场Qall和地震数据块输入到GPU。
具体地,CPU将分块得到的D个地震数据块逐块输入GPU内。
(16)GPU根据地震工区的Q值场Qall对D个地震数据块中的每个地震数据愉块进行反Q滤波处理,得到处理后的叠后地震数据体。
GPU处理完后,可以将处理结果输出到CPU端,CPU与GPU相互协同,共同完成所有分块数据的输入、反Q滤波处理与输出操作。本发明应用CPU-GPU平台对地震工区的所有地震数据处理完用时618秒,而现有技术中仅在CPU上执行相同数据的反Q滤波处理用时13600秒,通过对比可以看出,应用CPU-GPU平台后,计算效率显著提升,加速比达到22倍。
(17)GPU通过显示软件对处理后的叠后地震数据体进行图像显示。
GPU显示的处理后的叠后地震数据体的图像可以如图3-11所示,其中,该图3-11是对图3-9所示的地震剖面处理后的图,在该图3-11中,横轴表示垂直距离LS,单位为km,纵轴表示时间深度τ,单位为s。图3-12是图3-11中的区域3110的放大图。对比图3-9与图3-11,以及图3-10与图3-12可见,应用本发明后,地震数据的分辨率明显提高,原先叠合的同相轴380被更好的分离开来,同相轴380横向连续性变好。
图3-13是处理前的叠后地震数据体在1200ms时刻的水平时间切片,图3-14是应用本发明实施例提供的叠后地震数据体处理方法对图3-13所示的水平时间切片处理后的时间切片图,参见图3-13和图3-14,其中,横坐标Line表示地震主测线的序号,纵坐标表示共深度点(英文:Common Depth Point;简称:CDP),通过对比图3-13和图3-14可以看出,应用本发明实施例提供的叠后地震数据体处理方法对图3-13所示的叠后地震数据体在1200ms时刻的水平时间切片处理后,河道空间展布更加清楚,细节信息明显增多。
综上所述,本发明实施例提供的叠后地震数据体处理方法,通过从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线;从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场;根据M条样本测线对应的样本速度场,采用李庆忠经验公式确定M条样本测线对应的初始Q值场;对M条样本测线对应的初始Q值场进行校正,得到M条样本测线对应的校正Q值场;根据M条样本测线对应的校正Q值场和M条样本测线对应的加权系数场,确定M条样本测线对应的最终Q值场;根据M条样本测线对应的样本速度场和M条样本测线对应的最终Q值场,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式;根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场;向GPU发送地震工区的Q值场和叠后地震数据体,以便于GPU根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。由于针对叠后地震数据体对李庆忠经验公式进行了更新,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
进一步地,现有技术中,通常在CPU中对叠后地震数据体进行处理,由于CPU在一个计算周期内只能对一个地震数据道进行处理,且其为单道输入输出,读写频繁,处理效率较低。本发明实施例提供的叠后地震数据体处理方法,通过CPU对叠后地震数据体进行分块处理得到地震数据块,GPU根据地震数据块对叠后地震数据体进行处理,由于GPU在一个计算周期内可以对多个地震数据道进行处理,GPU的处理能力较强,本发明实施例通过采用GPU对叠后地震数据体进行反Q滤波处理,提高了对叠后地震数据体的处理速度和处理效率。
本发明实施例提供的叠后地震数据体处理方法可以提高叠后地震数据的分辨率,特别适用于处理大规模的三维偏移叠加地震数据体,对提高叠后地震数据体的分辨率具有重要应用价值。
下述为本发明装置实施例,可以用于执行本发明方法实施例。对于本发明装置实施例中未披露的细节,请参照本发明方法实施例。
请参考图4,其示出了本发明实施例提供的一种叠后地震数据体处理装置400的框图,该叠后地震数据体处理装置400可以用于执行图1所示实施例提供的叠后地震数据体处理方法或者图3-1所示实施例提供的叠后地震数据体处理方法中的部分方法,参见图4,该叠后地震数据体处理装置400可以包括:
选取模块401,用于从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线,M为大于或者等于2的整数,且M=0.003N,N为叠后地震数据体的地震主测线的条数,每条地震测线对应一个地震剖面;
第一确定模块402,用于从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场ν;
第二确定模块403,用于根据M条样本测线对应的样本速度场ν,采用李庆忠经验公式确定M条样本测线对应的初始Q值场Q0,李庆忠经验公式为Q=α·νβ,α=14,β=2.2;
校正模块404,用于对M条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1
第三确定模块405,用于根据M条样本测线对应的校正Q值场Q1和M条样本测线对应的加权系数场λfinal,确定M条样本测线对应的最终Q值场Q2
更新模块406,用于根据M条样本测线对应的样本速度场ν和M条样本测线对应的最终Q值场Q2,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式;
第四确定模块407,用于根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场Qall
发送模块408,用于向图形处理器GPU发送地震工区的Q值场Qall和叠后地震数据体,以便于GPU根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
综上所述,本发明实施例提供的叠后地震数据体处理装置,通过从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线;从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场;根据M条样本测线对应的样本速度场,采用李庆忠经验公式确定M条样本测线对应的初始Q值场;对M条样本测线对应的初始Q值场进行校正,得到M条样本测线对应的校正Q值场;根据M条样本测线对应的校正Q值场和M条样本测线对应的加权系数场,确定M条样本测线对应的最终Q值场;根据M条样本测线对应的样本速度场和M条样本测线对应的最终Q值场,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式;根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场;向GPU发送地震工区的Q值场和叠后地震数据体,以便于GPU根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。由于针对叠后地震数据体对李庆忠经验公式进行了更新,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
请参考图5,其示出了本发明实施例提供的另一种叠后地震数据体处理装置400的框图,该叠后地震数据体处理装置500可以用于执行图1所示实施例提供的叠后地震数据体处理方法或者图3-1所示实施例提供的叠后地震数据体处理方法中的部分方法,参见图5,该叠后地震数据体处理装置500可以包括但不限于:
选取模块501,用于从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线,M为大于或者等于2的整数,且M=0.003N,N为叠后地震数据体的地震主测线的条数,每条地震测线对应一个地震剖面;
第一确定模块502,用于从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场ν;
第二确定模块503,用于根据M条样本测线对应的样本速度场ν,采用李庆忠经验公式确定M条样本测线对应的初始Q值场Q0,李庆忠经验公式为Q=α·νβ,α=14,β=2.2;
校正模块504,用于对M条样本测线对应的初始Q值场Q0进行校正,得到M条样本测线对应的校正Q值场Q1
第三确定模块505,用于根据M条样本测线对应的校正Q值场Q1和M条样本测线对应的加权系数场λfinal,确定M条样本测线对应的最终Q值场Q2
更新模块506,用于根据M条样本测线对应的样本速度场ν和M条样本测线对应的最终Q值场Q2,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式;
第四确定模块507,用于根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场Qall
发送模块508,用于向图形处理器GPU发送地震工区的Q值场Qall和叠后地震数据体,以便于GPU根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
进一步地,请继续参考图5,该叠后地震数据体处理装置500还可以包括:
第一处理模块509,用于根据预设的加权系数组对M条样本测线对应的校正Q值场Q1进行处理,得到M条样本测线对应的加权Q值场组,预设的加权系数组中包括P个加权系数λ,加权Q值场组中包括P套加权Q值场,P为大于或者等于1的整数;
第二处理模块510,用于根据加权Q值场组对M条样本测线组成的样本地震数据体进行反Q滤波处理,得到反Q滤波数据体组,反Q滤波数据体组中包括P套反Q滤波数据体,每套反Q滤波数据体包括M条反Q样本测线,每条反Q样本测线对应一个地震剖面;
频谱分析模块511,用于采用傅里叶变换对反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析;
第五确定模块512,用于根据频谱分析结果进行处理,确定M条反Q样本测线对应的加权系数场;
第六确定模块513,用于将M条反Q样本测线对应的加权系数场确定为M条样本测线对应的加权系数场λfinal
进一步地,频谱分析模块511,具体用于:
在反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面上选择三个矩形时窗,三个矩形时窗包括:浅层矩形时窗、中层矩形时窗和深层矩形时窗;
采用傅里叶变换对每条反Q样本测线对应的地震剖面上位于三个矩形时窗中的每个矩形时窗内的部分进行处理,得到每条反Q样本测线对应的3个反Q样本频谱;
根据每套反Q滤波数据体中的M条反Q样本测线对应的反Q样本频谱,得到每套反Q滤波数据体对应的3M个反Q样本频谱;
根据反Q滤波数据体组中的P套反Q滤波数据体对应的反Q样本频谱得到3×M×P个反Q样本频谱,3×M×P个反Q样本频谱中包括:M×P个浅层矩形时窗对应的反Q样本频谱、M×P个中层矩形时窗对应的反Q样本频谱和M×P个深层矩形时窗对应的反Q样本频谱;
对3×M×P个反Q样本频谱中的每个反Q样本频谱进行频谱分析。
进一步地,第五确定模块512,具体用于:
根据频谱分析结果确定反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱,第一目标频谱为每条反Q样本测线对应的地震剖面上位于浅层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,第二目标频谱为每条反Q样本测线对应的地震剖面上位于中层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,第三目标频谱为每条反Q样本测线对应的地震剖面上位于深层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱;
分别确定第一目标频谱对应的矩形时窗对应的第一目标加权系数、第二目标频谱对应的矩形时窗对应的第二目标加权系数和第三目标频谱对应的矩形时窗对应的第三目标加权系数,得到每条反Q样本测线对应的3个目标加权系数,其中,第一目标加权系数、第二目标加权系数和第三目标加权系数都为预设的加权系数组中的加权系数;
根据M条反Q样本测线中的每条反Q样本测线对应的3个目标加权系数,得到M条反Q样本测线对应的3M个目标加权系数;
对3M个目标加权系数进行空间插值,得到地震工区对应的加权系数场;
对地震工区对应的加权系数场进行平滑处理,得到处理后的加权系数场;
从处理后的加权系数场中抽取出M条反Q样本测线对应的加权系数场。
进一步地,请继续参考图5,该叠后地震数据体处理装置500还可以包括:
采样模块514,用于按照预设采样间隔对地震工区的Q值场Qall进行采样,得到采样Q值场Qs
分块模块515,用于对叠后地震数据体进行分块处理,得到D个地震数据块,D为大于或者等于2的整数;
发送模块508,具体用于:
依次向GPU发送采样Q值场Qs和D个地震数据块中的每个地震数据块。
进一步地,分块模块515,具体用于:
计算叠后地震数据体中的地震数据的总道数E;
根据GPU在一个周期内处理的地震数据的道数Tracemax和叠后地震数据体中的地震数据的总道数E,对叠后地震数据体进行分块处理,得到D个地震数据块;
其中,E=Tracemax×(D-1)+L,L为D个地震数据块中的最后一个地震数据块中的地震数据的道数。
综上所述,本发明实施例提供的叠后地震数据体处理装置,通过从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线;从地震工区的叠加速度场中确定M条样本测线中的每条样本测线对应的速度场,得到M条样本测线对应的样本速度场;根据M条样本测线对应的样本速度场,采用李庆忠经验公式确定M条样本测线对应的初始Q值场;对M条样本测线对应的初始Q值场进行校正,得到M条样本测线对应的校正Q值场;根据M条样本测线对应的校正Q值场和M条样本测线对应的加权系数场,确定M条样本测线对应的最终Q值场;根据M条样本测线对应的样本速度场和M条样本测线对应的最终Q值场,采用最小二乘法更新李庆忠经验公式,得到更新后的李庆忠经验公式;根据地震工区的叠加速度场,采用更新后的李庆忠经验公式确定地震工区的Q值场;向GPU发送地震工区的Q值场和叠后地震数据体,以便于GPU根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。由于针对叠后地震数据体对李庆忠经验公式进行了更新,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
进一步地,现有技术中,通常在CPU中对叠后地震数据体进行处理,由于CPU在一个计算周期内只能对一个地震数据道进行处理,且其为单道输入输出,读写频繁,处理效率较低。本发明实施例提供的叠后地震数据体处理装置,通过CPU对叠后地震数据体进行分块处理得到地震数据块,GPU根据地震数据块对叠后地震数据体进行处理,由于GPU在一个计算周期内可以对多个地震数据道进行处理,GPU的处理能力较强,本发明实施例通过采用GPU对叠后地震数据体进行反Q滤波处理,提高了对叠后地震数据体的处理速度和处理效率。
本发明实施例提供的叠后地震数据体处理装置可以提高叠后地震数据的分辨率,特别适用于处理大规模的三维偏移叠加地震数据体,对提高叠后地震数据体的分辨率具有重要应用价值。
请参考图6,其示出了本发明实施例提供的再一种叠后地震数据体处理装置400的框图,该叠后地震数据体处理装置600可以用于执行图2所示实施例提供的叠后地震数据体处理方法或者图3-1所示实施例提供的叠后地震数据体处理方法中的部分方法,参见图6,该叠后地震数据体处理装置600可以包括:
接收模块610,用于接收中央处理器CPU发送的地震工区的Q值场Qall和叠后地震数据体;
处理模块620,用于根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
综上所述,本发明实施例提供的叠后地震数据体处理装置,通过接收CPU发送的地震工区的Q值场和叠后地震数据体,根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体,由于地震工区的Q值场是CPU根据更新后的李庆忠经验公式计算得到的,且更新后的李庆忠经验公式是CPU针对叠后地震数据体对李庆忠经验公式进行更新得到的,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
请参考图7-1,其示出了本发明实施例提供的又一种叠后地震数据体处理装置400的框图,该叠后地震数据体处理装置700可以用于执行图2所示实施例提供的叠后地震数据体处理方法或者图3-1所示实施例提供的叠后地震数据体处理方法中的部分方法,参见图7-1,该叠后地震数据体处理装置700可以包括:
接收模块710,用于接收中央处理器CPU发送的地震工区的Q值场Qall和叠后地震数据体;
处理模块720,用于根据地震工区的Q值场Qall对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
进一步地,接收模块710,具体用于:
接收CPU发送的采样Q值场Qs和D个地震数据块中的每个地震数据块,采样Q值场Qs是CPU按照预设采样间隔对地震工区的Q值场Qall进行采样得到的,D个地震数据块是CPU对叠后地震数据体进行分块处理得到的;
请参考图7-2,其示出的是图7-1所示实施例提供的一种处理模块720的框图,参见图7-2,该处理模块720可以包括:
处理子模块721,用于根据采样Q值场Qs对D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块;
得到子模块722,用于根据D个处理后的地震数据块,得到处理后的叠后地震数据体。
进一步地,请参考图7-3,其示出的是图7-1所示实施例提供的一种处理子模块721的框图,参见图7-3,该处理子模块721可以包括:
第一处理单元7211,用于对采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ
第二处理单元7212,用于根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到D个处理后的地震数据块。
进一步地,第一处理单元7211,具体用于:
根据采样Q值场Qs和Q值插值计算公式确定每个时间深度τ处的Q值Qτ
其中,Q值插值计算公式为:
Q τ j = Q τ ( j - 1 ) ( τ ( j + 1 ) - τ j ) + Q τ ( j + 1 ) ( τ j - τ ( j - 1 ) ) τ ( j + 1 ) - τ ( j - 1 ) ;
τ ( j - 1 ) = ( int ( τ j 0.02 ) - 1 ) × 0.02 , τ ( j + 1 ) = ( int ( τ j 0.02 ) + 1 ) × 0.02 , τj、τ(j-1)和τ(j+1)都为时间深度,为时间深度τj处的Q值插值,为时间深度τ(j-1)处的Q值插值,为时间深度τ(j+1)处的Q值插值,表示对取整。
进一步地,请参考图7-4,其示出的是图7-1所示实施例提供的一种第二处理单元7212的框图,参见图7-4,该第二处理单元7212可以包括:
处理子单元72121,用于根据所有时间深度τ处的Q值Qτ对D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道;
第一得到子单元72122,用于根据每个地震数据块中的所有处理后的地震数据道,得到处理后的地震数据块;
第二得到子单元72123,用于根据D个地震数据块中的所有处理后的地震数据块,得到D个处理后的地震数据块。
进一步地,请参考图7-5,其示出的是图7-1所示实施例提供的一种处理子单元72121的框图,参见图7-5,该处理子单元72121可以包括:
第一确定子单元721211,用于分别确定GPU中的线程的线程格尺寸和线程块尺寸;
变换子单元721212,用于对D个地震数据块中的每个地震数据块中的每个地震数据道进行快速傅里叶变换;
第二确定子单元721213,用于根据快速傅里叶变换结果,确定D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
第三确定子单元721214,用于根据角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定每个地震数据道对应的起算角频率值ωp和终止角频率值ωz
第一计算子单元721215,用于根据GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
第二计算子单元721216,用于根据每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
第一处理子单元72127,用于根据GPU中的所有线程处理得到的每个地震数据道对应的所有时间深度处的所有频率点对应的地震幅值,处理得到处理后的地震数据道。
进一步地,GPU中的线程的线程格尺寸为:Tracemax为GPU在一个周期内处理的地震数据的道数,b_size=2m,m为正整数,且4≤m≤6,S为地震数据道的采样点数;
GPU中的线程的线程块尺寸为:(b_size,b_size),b_size=2m,m为正整数,且4≤m≤6;
第二确定子单元721213,具体用于:
根据快速傅里叶变换结果,采用角频率采样间隔计算公式确定D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
其中,角频率采样间隔计算公式为:△ω=2π/(Gdt),G为满足快速傅里叶变换算法的时间采样间隔,dt为地震数据的时间采样间隔;
第三确定子单元721214,具体用于:
根据目标低角频率值ωmin和角频率采样间隔△ω,采用起算角频率计算公式计算每个地震数据道对应的起算角频率值ωp
根据目标高角频率值ωmax和角频率采样间隔△ω,采用终止角频率计算公式计算每个地震数据道对应的终止角频率值ωz
其中,起算角频率计算公式为:终止角频率计算公式为:△f为频率增量,ωmax+2π△f为通过反Q滤波处理后期望达到的角频率的上限;
第一计算子单元721215,具体用于:
根据GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,采用地震幅值计算公式计算每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
其中,地震幅值计算公式为: A为每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且A为复数,k为每个频率点的序号,△ω为角频率采样间隔,τj为时间深度,为时间深度τj处的Q值;
第二计算子单元721216,具体用于:
根据每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,采用幅值和计算公式计算每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
其中,幅值和计算公式为:
Σ A = Σ k = int ( ω min Δ ω ) int ( ω max + 2 π Δ f Δ ω ) Re { F ( k Δ ω ) exp ( kΔωτ j 2 Q τ j ) } ,
A为每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且A为复数,k为每个频率点的序号,△ω为角频率采样间隔,τj为时间深度,为时间深度τj处的Q值,Re表示取实部。
综上所述,本发明实施例提供的叠后地震数据体处理装置,通过接收CPU发送的地震工区的Q值场和叠后地震数据体,根据地震工区的Q值场对叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体,由于地震工区的Q值场是CPU根据更新后的李庆忠经验公式计算得到的,且更新后的李庆忠经验公式是CPU针对叠后地震数据体对李庆忠经验公式进行更新得到的,因此,解决了对叠后地震数据体处理的灵活性较低、可靠性较差的问题,达到了提高叠后地震数据体处理的灵活性和可靠性的效果。
进一步地,现有技术中,通常在CPU中对叠后地震数据体进行处理,由于CPU在一个计算周期内只能对一个地震数据道进行处理,且其为单道输入输出,读写频繁,处理效率较低。本发明实施例提供的叠后地震数据体处理装置,通过CPU对叠后地震数据体进行分块处理得到地震数据块,GPU根据地震数据块对叠后地震数据体进行处理,由于GPU在一个计算周期内可以对多个地震数据道进行处理,GPU的处理能力较强,本发明实施例通过采用GPU对叠后地震数据体进行反Q滤波处理,提高了对叠后地震数据体的处理速度和处理效率。
本发明实施例提供的叠后地震数据体处理装置可以提高叠后地震数据的分辨率,特别适用于处理大规模的三维偏移叠加地震数据体,对提高叠后地震数据体的分辨率具有重要应用价值。
需要说明的是:上述实施例提供的叠后地震数据体处理装置在处理叠后地震数据体时,仅以上述各功能模块的划分进行举例说明,实际应用中,可以根据需要而将上述功能分配由不同的功能模块完成,即将设备的内部结构划分成不同的功能模块,以完成以上描述的全部或者部分功能。另外,上述实施例提供的叠后地震数据体处理装置与叠后地震数据体处理方法实施例属于同一构思,其具体实现过程详见方法实施例,这里不再赘述。
本领域普通技术人员可以理解实现上述实施例的全部或部分步骤可以通过硬件来完成,也可以通过程序来指令相关的硬件完成,所述的程序可以存储于一种计算机可读存储介质中,上述提到的存储介质可以是只读存储器,磁盘或光盘等。
以上所述仅为本发明的较佳实施例,并不用以限制本发明,凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (26)

1.一种叠后地震数据体处理方法,其特征在于,用于中央处理器CPU,所述方法包括:
从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线,所述M为大于或者等于2的整数,且M=0.003N,所述N为所述叠后地震数据体的地震主测线的条数,每条所述地震测线对应一个地震剖面;
从所述地震工区的叠加速度场中确定所述M条样本测线中的每条样本测线对应的速度场,得到所述M条样本测线对应的样本速度场ν;
根据所述M条样本测线对应的样本速度场ν,采用李庆忠经验公式确定所述M条样本测线对应的初始Q值场Q0,所述李庆忠经验公式为Q=α·νβ,α=14,β=2.2;
对所述M条样本测线对应的初始Q值场Q0进行校正,得到所述M条样本测线对应的校正Q值场Q1
根据所述M条样本测线对应的校正Q值场Q1和所述M条样本测线对应的加权系数场λfinal,确定所述M条样本测线对应的最终Q值场Q2
根据所述M条样本测线对应的样本速度场ν和所述M条样本测线对应的最终Q值场Q2,采用最小二乘法更新所述李庆忠经验公式,得到更新后的李庆忠经验公式;
根据所述地震工区的叠加速度场,采用所述更新后的李庆忠经验公式确定所述地震工区的Q值场Qall
向图形处理器GPU发送所述地震工区的Q值场Qall和所述叠后地震数据体,以便于所述GPU根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
2.根据权利要求1所述的方法,其特征在于,在所述根据所述M条样本测线对应的校正Q值场Q1和所述M条样本测线对应的加权系数场λfinal,确定所述M条样本测线对应的最终Q值场Q2之前,所述方法还包括:
根据预设的加权系数组对所述M条样本测线对应的校正Q值场Q1进行处理,得到所述M条样本测线对应的加权Q值场组,所述预设的加权系数组中包括P个加权系数λ,所述加权Q值场组中包括P套加权Q值场,所述P为大于或者等于1的整数;
根据所述加权Q值场组对所述M条样本测线组成的样本地震数据体进行反Q滤波处理,得到反Q滤波数据体组,所述反Q滤波数据体组中包括P套反Q滤波数据体,每套反Q滤波数据体包括M条反Q样本测线,每条所述反Q样本测线对应一个地震剖面;
采用傅里叶变换对所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析;
根据频谱分析结果进行处理,确定所述M条反Q样本测线对应的加权系数场;
将所述M条反Q样本测线对应的加权系数场确定为所述M条样本测线对应的加权系数场λfinal
3.根据权利要求2所述的方法,其特征在于,
所述采用傅里叶变换对所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析,包括:
在所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面上选择三个矩形时窗,所述三个矩形时窗包括:浅层矩形时窗、中层矩形时窗和深层矩形时窗;
采用傅里叶变换对所述每条反Q样本测线对应的地震剖面上位于所述三个矩形时窗中的每个矩形时窗内的部分进行处理,得到所述每条反Q样本测线对应的3个反Q样本频谱;
根据所述每套反Q滤波数据体中的M条反Q样本测线对应的反Q样本频谱,得到所述每套反Q滤波数据体对应的3M个反Q样本频谱;
根据所述反Q滤波数据体组中的P套反Q滤波数据体对应的反Q样本频谱得到3×M×P个反Q样本频谱,所述3×M×P个反Q样本频谱中包括:M×P个浅层矩形时窗对应的反Q样本频谱、M×P个中层矩形时窗对应的反Q样本频谱和M×P个深层矩形时窗对应的反Q样本频谱;
对所述3×M×P个反Q样本频谱中的每个反Q样本频谱进行频谱分析。
4.根据权利要求3所述的方法,其特征在于,所述根据频谱分析结果进行处理,确定所述M条反Q样本测线对应的加权系数场,包括:
根据频谱分析结果确定所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱,所述第一目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述浅层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,所述第二目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述中层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,所述第三目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述深层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱;
分别确定所述第一目标频谱对应的矩形时窗对应的第一目标加权系数、所述第二目标频谱对应的矩形时窗对应的第二目标加权系数和所述第三目标频谱对应的矩形时窗对应的第三目标加权系数,得到所述每条反Q样本测线对应的3个目标加权系数,其中,所述第一目标加权系数、所述第二目标加权系数和所述第三目标加权系数都为所述预设的加权系数组中的加权系数;
根据所述M条反Q样本测线中的每条反Q样本测线对应的3个目标加权系数,得到所述M条反Q样本测线对应的3M个目标加权系数;
对所述3M个目标加权系数进行空间插值,得到所述地震工区对应的加权系数场;
对所述地震工区对应的加权系数场进行平滑处理,得到处理后的加权系数场;
从所述处理后的加权系数场中抽取出所述M条反Q样本测线对应的加权系数场。
5.根据权利要求1至4任一所述的方法,其特征在于,
在所述向图形处理器GPU发送所述地震工区的Q值场Qall和所述叠后地震数据体之前,所述方法还包括:
按照预设采样间隔对所述地震工区的Q值场Qall进行采样,得到采样Q值场Qs
对所述叠后地震数据体进行分块处理,得到D个地震数据块,所述D为大于或者等于2的整数;
所述向图形处理器GPU发送所述地震工区的Q值场Qall和叠后地震数据体,包括:
依次向所述GPU发送所述采样Q值场Qs和所述D个地震数据块中的每个地震数据块。
6.根据权利要求5所述的方法,其特征在于,所述对所述叠后地震数据体进行分块处理,得到D个地震数据块,包括:
计算所述叠后地震数据体中的地震数据的总道数E;
根据所述GPU在一个周期内处理的地震数据的道数Tracemax和所述叠后地震数据体中的地震数据的总道数E,对所述叠后地震数据体进行分块处理,得到所述D个地震数据块;
其中,E=Tracemax×(D-1)+L,所述L为所述D个地震数据块中的最后一个地震数据块中的地震数据的道数。
7.一种叠后地震数据体处理方法,其特征在于,用于图形处理器GPU,所述方法包括:
接收中央处理器CPU发送的地震工区的Q值场Qall和叠后地震数据体;
根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
8.根据权利要求7所述的方法,其特征在于,
所述接收中央处理器CPU发送的地震工区的Q值场Qall和叠后地震数据体,包括:
接收所述CPU发送的采样Q值场Qs和D个地震数据块中的每个地震数据块,所述采样Q值场Qs是所述CPU按照预设采样间隔对所述地震工区的Q值场Qall进行采样得到的,所述D个地震数据块是所述CPU对所述叠后地震数据体进行分块处理得到的;
所述根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体,包括:
根据所述采样Q值场Qs对所述D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块;
根据所述D个处理后的地震数据块,得到所述处理后的叠后地震数据体。
9.根据权利要求8所述的方法,其特征在于,
所述根据所述采样Q值场Qs对所述D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块,包括:
对所述采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ
根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到所述D个处理后的地震数据块。
10.根据权利要求9所述的方法,其特征在于,所述对所述采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ,包括:
根据所述采样Q值场Qs和Q值插值计算公式确定每个时间深度τ处的Q值Qτ
其中,所述Q值插值计算公式为:
Q τ j = Q τ ( j - 1 ) ( τ ( j + 1 ) - τ j ) + Q τ ( j + 1 ) ( τ j - τ ( j - 1 ) ) τ ( j + 1 ) - τ ( j - 1 ) ;
τ ( j - 1 ) = ( int ( τ j 0.02 ) - 1 ) × 0.02 , τ ( j + 1 ) = ( int ( τ j 0.02 ) + 1 ) × 0.02 , 所述τj、所述τ(j-1)和所述τ(j+1)都为时间深度,所述为时间深度τj处的Q值插值,所述为时间深度τ(j-1)处的Q值插值,所述为时间深度τ(j+1)处的Q值插值,表示对所述取整。
11.根据权利要求9或10所述的方法,其特征在于,所述根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到所述D个处理后的地震数据块,包括:
根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道;
根据所述每个地震数据块中的所有处理后的地震数据道,得到处理后的地震数据块;
根据所述D个地震数据块中的所有处理后的地震数据块,得到所述D个处理后的地震数据块。
12.根据权利要求11所述的方法,其特征在于,
所述根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道,包括:
分别确定所述GPU中的线程的线程格尺寸和线程块尺寸;
对所述D个地震数据块中的每个地震数据块中的每个地震数据道进行快速傅里叶变换;
根据快速傅里叶变换结果,确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
根据所述角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz
根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
根据所述GPU中的所有线程处理得到的所述每个地震数据道对应的所有时间深度处的所有频率点对应的地震幅值,处理得到所述处理后的地震数据道。
13.根据权利要求12所述的方法,其特征在于,
所述GPU中的线程的线程格尺寸为:所述Tracemax为所述GPU在一个周期内处理的地震数据的道数,b_size=2m,所述m为正整数,且4≤m≤6,所述S为地震数据道的采样点数;
所述GPU中的线程的线程块尺寸为:(b_size,b_size),b_size=2m,所述m为正整数,且4≤m≤6;
所述根据快速傅里叶变换结果,确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω,包括:
根据快速傅里叶变换结果,采用角频率采样间隔计算公式确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
其中,所述角频率采样间隔计算公式为:△ω=2π/(Gdt),所述G为满足快速傅里叶变换算法的时间采样间隔,所述dt为地震数据的时间采样间隔;
所述根据所述角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,包括:
根据所述目标低角频率值ωmin和所述角频率采样间隔△ω,采用起算角频率计算公式计算所述每个地震数据道对应的起算角频率值ωp
根据所述目标高角频率值ωmax和所述角频率采样间隔△ω,采用终止角频率计算公式计算所述每个地震数据道对应的终止角频率值ωz
其中,所述起算角频率计算公式为:所述终止角频率计算公式为:所述△f为频率增量,所述ωmax+2π△f为通过反Q滤波处理后期望达到的角频率的上限;
所述根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,包括:
根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,采用地震幅值计算公式计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
其中,所述地震幅值计算公式为: 所述A为所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且所述A为复数,所述k为所述每个频率点的序号,所述△ω为所述角频率采样间隔,所述τj为时间深度,所述为时间深度τj处的Q值;
所述根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A,包括:
根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,采用幅值和计算公式计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
其中,所述幅值和计算公式为:
Σ A = Σ k = int ( ω min Δ ω ) int ( ω max + 2 π Δ f Δ ω ) Re { F ( k Δ ω ) exp ( kΔωτ j 2 Q τ j ) } ,
所述A为所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且所述A为复数,所述k为所述每个频率点的序号,所述△ω为所述角频率采样间隔,所述τj为时间深度,所述为时间深度τj处的Q值,所述Re表示取实部。
14.一种叠后地震数据体处理装置,其特征在于,用于中央处理器CPU,所述装置包括:
选取模块,用于从地震工区的叠后地震数据体对应的地震测线中选取M条样本测线,所述M为大于或者等于2的整数,且M=0.003N,所述N为所述叠后地震数据体的地震主测线的条数,每条所述地震测线对应一个地震剖面;
第一确定模块,用于从所述地震工区的叠加速度场中确定所述M条样本测线中的每条样本测线对应的速度场,得到所述M条样本测线对应的样本速度场ν;
第二确定模块,用于根据所述M条样本测线对应的样本速度场ν,采用李庆忠经验公式确定所述M条样本测线对应的初始Q值场Q0,所述李庆忠经验公式为Q=α·νβ,α=14,β=2.2;
校正模块,用于对所述M条样本测线对应的初始Q值场Q0进行校正,得到所述M条样本测线对应的校正Q值场Q1
第三确定模块,用于根据所述M条样本测线对应的校正Q值场Q1和所述M条样本测线对应的加权系数场λfinal,确定所述M条样本测线对应的最终Q值场Q2
更新模块,用于根据所述M条样本测线对应的样本速度场ν和所述M条样本测线对应的最终Q值场Q2,采用最小二乘法更新所述李庆忠经验公式,得到更新后的李庆忠经验公式;
第四确定模块,用于根据所述地震工区的叠加速度场,采用所述更新后的李庆忠经验公式确定所述地震工区的Q值场Qall
发送模块,用于向图形处理器GPU发送所述地震工区的Q值场Qall和所述叠后地震数据体,以便于所述GPU根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
15.根据权利要求14所述的装置,其特征在于,所述装置还包括:
第一处理模块,用于根据预设的加权系数组对所述M条样本测线对应的校正Q值场Q1进行处理,得到所述M条样本测线对应的加权Q值场组,所述预设的加权系数组中包括P个加权系数λ,所述加权Q值场组中包括P套加权Q值场,所述P为大于或者等于1的整数;
第二处理模块,用于根据所述加权Q值场组对所述M条样本测线组成的样本地震数据体进行反Q滤波处理,得到反Q滤波数据体组,所述反Q滤波数据体组中包括P套反Q滤波数据体,每套反Q滤波数据体包括M条反Q样本测线,每条所述反Q样本测线对应一个地震剖面;
频谱分析模块,用于采用傅里叶变换对所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面进行频谱分析;
第五确定模块,用于根据频谱分析结果进行处理,确定所述M条反Q样本测线对应的加权系数场;
第六确定模块,用于将所述M条反Q样本测线对应的加权系数场确定为所述M条样本测线对应的加权系数场λfinal
16.根据权利要求15所述的装置,其特征在于,
所述频谱分析模块,具体用于:
在所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的地震剖面上选择三个矩形时窗,所述三个矩形时窗包括:浅层矩形时窗、中层矩形时窗和深层矩形时窗;
采用傅里叶变换对所述每条反Q样本测线对应的地震剖面上位于所述三个矩形时窗中的每个矩形时窗内的部分进行处理,得到所述每条反Q样本测线对应的3个反Q样本频谱;
根据所述每套反Q滤波数据体中的M条反Q样本测线对应的反Q样本频谱,得到所述每套反Q滤波数据体对应的3M个反Q样本频谱;
根据所述反Q滤波数据体组中的P套反Q滤波数据体对应的反Q样本频谱得到3×M×P个反Q样本频谱,所述3×M×P个反Q样本频谱中包括:M×P个浅层矩形时窗对应的反Q样本频谱、M×P个中层矩形时窗对应的反Q样本频谱和M×P个深层矩形时窗对应的反Q样本频谱;
对所述3×M×P个反Q样本频谱中的每个反Q样本频谱进行频谱分析。
17.根据权利要求16所述的装置,其特征在于,
所述第五确定模块,具体用于:
根据频谱分析结果确定所述反Q滤波数据体组中的每套反Q滤波数据体中的每条反Q样本测线对应的第一目标频谱、第二目标频谱和第三目标频谱,所述第一目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述浅层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,所述第二目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述中层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱,所述第三目标频谱为所述每条反Q样本测线对应的地震剖面上位于所述深层矩形时窗内的部分对应的P个反Q样本频谱中与横轴围成的封闭区域的面积最大的反Q样本频谱;
分别确定所述第一目标频谱对应的矩形时窗对应的第一目标加权系数、所述第二目标频谱对应的矩形时窗对应的第二目标加权系数和所述第三目标频谱对应的矩形时窗对应的第三目标加权系数,得到所述每条反Q样本测线对应的3个目标加权系数,其中,所述第一目标加权系数、所述第二目标加权系数和所述第三目标加权系数都为所述预设的加权系数组中的加权系数;
根据所述M条反Q样本测线中的每条反Q样本测线对应的3个目标加权系数,得到所述M条反Q样本测线对应的3M个目标加权系数;
对所述3M个目标加权系数进行空间插值,得到所述地震工区对应的加权系数场;
对所述地震工区对应的加权系数场进行平滑处理,得到处理后的加权系数场;
从所述处理后的加权系数场中抽取出所述M条反Q样本测线对应的加权系数场。
18.根据权利要求14至17任一所述的装置,其特征在于,所述装置还包括:
采样模块,用于按照预设采样间隔对所述地震工区的Q值场Qall进行采样,得到采样Q值场Qs
分块模块,用于对所述叠后地震数据体进行分块处理,得到D个地震数据块,所述D为大于或者等于2的整数;
所述发送模块,具体用于:
依次向所述GPU发送所述采样Q值场Qs和所述D个地震数据块中的每个地震数据块。
19.根据权利要求18所述的装置,其特征在于,
所述分块模块,具体用于:
计算所述叠后地震数据体中的地震数据的总道数E;
根据所述GPU在一个周期内处理的地震数据的道数Tracemax和所述叠后地震数据体中的地震数据的总道数E,对所述叠后地震数据体进行分块处理,得到所述D个地震数据块;
其中,E=Tracemax×(D-1)+L,所述L为所述D个地震数据块中的最后一个地震数据块中的地震数据的道数。
20.一种叠后地震数据体处理装置,其特征在于,用于图形处理器GPU,所述装置包括:
接收模块,用于接收中央处理器CPU发送的地震工区的Q值场Qall和叠后地震数据体;
处理模块,用于根据所述地震工区的Q值场Qall对所述叠后地震数据体进行反Q滤波处理,得到处理后的叠后地震数据体。
21.根据权利要求20所述的装置,其特征在于,
所述接收模块,具体用于:
接收所述CPU发送的采样Q值场Qs和D个地震数据块中的每个地震数据块,所述采样Q值场Qs是所述CPU按照预设采样间隔对所述地震工区的Q值场Qall进行采样得到的,所述D个地震数据块是所述CPU对所述叠后地震数据体进行分块处理得到的;
所述处理模块,包括:
处理子模块,用于根据所述采样Q值场Qs对所述D个地震数据块中的每个地震数据块进行反Q滤波处理,得到D个处理后的地震数据块;
得到子模块,用于根据所述D个处理后的地震数据块,得到所述处理后的叠后地震数据体。
22.根据权利要求21所述的装置,其特征在于,
所述处理子模块,包括:
第一处理单元,用于对所述采样Q值场Qs进行处理得到每个时间深度τ处的Q值Qτ
第二处理单元,用于根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的所有时间深度τ处的地震幅值进行反Q滤波处理,得到所述D个处理后的地震数据块。
23.根据权利要求22所述的装置,其特征在于,
所述第一处理单元,具体用于:
根据所述采样Q值场Qs和Q值插值计算公式确定每个时间深度τ处的Q值Qτ
其中,所述Q值插值计算公式为:
Q τ j = Q τ ( j - 1 ) ( τ ( j + 1 ) - τ j ) + Q τ ( j + 1 ) ( τ j - τ ( j - 1 ) ) τ ( j + 1 ) - τ ( j - 1 ) ;
τ ( j - 1 ) = ( int ( τ j 0.02 ) - 1 ) × 0.02 , τ ( j + 1 ) = ( int ( τ j 0.02 ) + 1 ) × 0.02 , 所述τj、所述τ(j-1)和所述τ(j+1)都为时间深度,所述为时间深度τj处的Q值插值,所述为时间深度τ(j-1)处的Q值插值,所述为时间深度τ(j+1)处的Q值插值,表示对所述取整。
24.根据权利要求22或23所述的装置,其特征在于,所述第二处理单元,包括:
处理子单元,用于根据所有时间深度τ处的Q值Qτ对所述D个地震数据块中的每个地震数据块中的每个地震数据道的所有时间深度τ处的地震幅值进行反Q滤波处理,得到处理后的地震数据道;
第一得到子单元,用于根据所述每个地震数据块中的所有处理后的地震数据道,得到处理后的地震数据块;
第二得到子单元,用于根据所述D个地震数据块中的所有处理后的地震数据块,得到所述D个处理后的地震数据块。
25.根据权利要求24所述的装置,其特征在于,
所述处理子单元,包括:
第一确定子单元,用于分别确定所述GPU中的线程的线程格尺寸和线程块尺寸;
变换子单元,用于对所述D个地震数据块中的每个地震数据块中的每个地震数据道进行快速傅里叶变换;
第二确定子单元,用于根据快速傅里叶变换结果,确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
第三确定子单元,用于根据所述角频率采样间隔△ω和预先接收的目标低角频率值ωmin、目标高角频率值ωmax,确定所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz
第一计算子单元,用于根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
第二计算子单元,用于根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
第一处理子单元,用于根据所述GPU中的所有线程处理得到的所述每个地震数据道对应的所有时间深度处的所有频率点对应的地震幅值,处理得到所述处理后的地震数据道。
26.根据权利要求25所述的装置,其特征在于,
所述GPU中的线程的线程格尺寸为:所述Tracemax为所述GPU在一个周期内处理的地震数据的道数,b_size=2m,所述m为正整数,且4≤m≤6,所述S为地震数据道的采样点数;
所述GPU中的线程的线程块尺寸为:(b_size,b_size),b_size=2m,所述m为正整数,且4≤m≤6;
所述第二确定子单元,具体用于:
根据快速傅里叶变换结果,采用角频率采样间隔计算公式确定所述D个地震数据块中所有地震数据对应的角频率采样间隔△ω;
其中,所述角频率采样间隔计算公式为:△ω=2π/(Gdt),所述G为满足快速傅里叶变换算法的时间采样间隔,所述dt为地震数据的时间采样间隔;
所述第三确定子单元,具体用于:
根据所述目标低角频率值ωmin和所述角频率采样间隔△ω,采用起算角频率计算公式计算所述每个地震数据道对应的起算角频率值ωp
根据所述目标高角频率值ωmax和所述角频率采样间隔△ω,采用终止角频率计算公式计算所述每个地震数据道对应的终止角频率值ωz
其中,所述起算角频率计算公式为:所述终止角频率计算公式为:所述△f为频率增量,所述ωmax+2π△f为通过反Q滤波处理后期望达到的角频率的上限;
所述第一计算子单元,具体用于:
根据所述GPU中的线程的线程格尺寸、线程块尺寸、所有时间深度τ处的Q值Qτ中的时间深度τj处的Q值所述每个地震数据道对应的起算角频率值ωp和终止角频率值ωz,采用地震幅值计算公式计算所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A;
其中,所述地震幅值计算公式为: 所述A为所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且所述A为复数,所述k为所述每个频率点的序号,所述△ω为所述角频率采样间隔,所述τj为时间深度,所述为时间深度τj处的Q值;
所述第二计算子单元,具体用于:
根据所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值A,采用幅值和计算公式计算所述每个地震数据道对应的时间深度τj处的所有频率点对应的地震幅值∑A;
其中,所述幅值和计算公式为:
Σ A = Σ k = int ( ω min Δ ω ) int ( ω max + 2 π Δ f Δ ω ) Re { F ( k Δ ω ) exp ( kΔωτ j 2 Q τ j ) } ,
所述A为所述每个地震数据道对应的时间深度τj处的每个频率点对应的地震幅值且所述A为复数,所述k为所述每个频率点的序号,所述△ω为所述角频率采样间隔,所述τj为时间深度,所述为时间深度τj处的Q值,所述Re表示取实部。
CN201610058962.1A 2016-01-28 2016-01-28 叠后地震数据体处理方法及装置 Active CN106257309B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610058962.1A CN106257309B (zh) 2016-01-28 2016-01-28 叠后地震数据体处理方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610058962.1A CN106257309B (zh) 2016-01-28 2016-01-28 叠后地震数据体处理方法及装置

Publications (2)

Publication Number Publication Date
CN106257309A true CN106257309A (zh) 2016-12-28
CN106257309B CN106257309B (zh) 2018-11-16

Family

ID=57713760

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610058962.1A Active CN106257309B (zh) 2016-01-28 2016-01-28 叠后地震数据体处理方法及装置

Country Status (1)

Country Link
CN (1) CN106257309B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108415076A (zh) * 2018-02-08 2018-08-17 中国地质调查局油气资源调查中心 一种基于线性反演的保幅保边界信噪增强方法
CN109100794A (zh) * 2017-06-20 2018-12-28 中国石油化工股份有限公司 一种时窗加权的相干速度反演方法及系统
CN111060969A (zh) * 2019-12-25 2020-04-24 恒泰艾普(北京)能源科技研究院有限公司 一种井控q补偿方法
CN111880218A (zh) * 2020-07-13 2020-11-03 西南石油大学 基于品质因子的反演子波字典构建方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100299071A1 (en) * 2007-12-14 2010-11-25 Kiyashchenko Denis Method of processing data obtained from seismic prospecting
CN102053273A (zh) * 2009-10-29 2011-05-11 中国石油化工股份有限公司 一种对地震波信号进行反q滤波的方法
CN103424777A (zh) * 2013-07-01 2013-12-04 中国科学院地质与地球物理研究所 一种提高地震成像分辨率的方法
CN104635268A (zh) * 2015-03-09 2015-05-20 成都晶石石油科技有限公司 地震资料约束下的品质因子计算方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100299071A1 (en) * 2007-12-14 2010-11-25 Kiyashchenko Denis Method of processing data obtained from seismic prospecting
CN102053273A (zh) * 2009-10-29 2011-05-11 中国石油化工股份有限公司 一种对地震波信号进行反q滤波的方法
CN103424777A (zh) * 2013-07-01 2013-12-04 中国科学院地质与地球物理研究所 一种提高地震成像分辨率的方法
CN104635268A (zh) * 2015-03-09 2015-05-20 成都晶石石油科技有限公司 地震资料约束下的品质因子计算方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
余连勇 等: "稳定的反Q滤波统一算法及其在地震资料高分辨率处理中的应用", 《中国海上油气》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109100794A (zh) * 2017-06-20 2018-12-28 中国石油化工股份有限公司 一种时窗加权的相干速度反演方法及系统
CN109100794B (zh) * 2017-06-20 2020-12-01 中国石油化工股份有限公司 一种时窗加权的相干速度反演方法及系统
CN108415076A (zh) * 2018-02-08 2018-08-17 中国地质调查局油气资源调查中心 一种基于线性反演的保幅保边界信噪增强方法
CN108415076B (zh) * 2018-02-08 2019-08-30 中国地质调查局油气资源调查中心 一种基于线性反演的保幅保边界信噪增强方法
CN111060969A (zh) * 2019-12-25 2020-04-24 恒泰艾普(北京)能源科技研究院有限公司 一种井控q补偿方法
CN111880218A (zh) * 2020-07-13 2020-11-03 西南石油大学 基于品质因子的反演子波字典构建方法

Also Published As

Publication number Publication date
CN106257309B (zh) 2018-11-16

Similar Documents

Publication Publication Date Title
RU2694621C1 (ru) Способ и устройство для обработки сейсмических данных
EP2715405B1 (en) Method of processing seismic data by providing surface offset common image gathers
CN105277978B (zh) 一种确定近地表速度模型的方法及装置
CN102890290B (zh) 一种起伏地表条件下的叠前深度偏移方法
CN102841379B (zh) 一种基于共散射点道集的叠前时间偏移与速度分析方法
KR20180067650A (ko) 진폭 보존을 갖는 fwi 모델 도메인 각도 스택들
CN106257309B (zh) 叠后地震数据体处理方法及装置
CN105182408A (zh) 一种合成地震记录的制作方法和装置
Mordret et al. Helmholtz tomography of ambient noise surface wave data to estimate Scholte wave phase velocity at Valhall Life of the Field
CN103954992B (zh) 一种反褶积方法及装置
CN107065013B (zh) 一种地震尺度下的层速度确定方法及装置
EP2877880B1 (en) System and method for migration velocity modeling
CN105093301B (zh) 共成像点反射角角道集的生成方法及装置
CN101545986A (zh) 基于最大能量旅行时计算的三维积分叠前深度偏移方法
CN104459794A (zh) 共反射点道集时变时间差值的校正方法及装置
US20160161619A1 (en) Systems and methods for aligning a monitor seismic survey with a baseline seismic survey
CN106597545B (zh) 一种水平裂缝地震叠前反演方法和装置
CN109884710A (zh) 针对激发井深设计的微测井层析成像方法
CN109188520A (zh) 薄储层厚度预测方法及装置
CN107942379A (zh) 一种提高复杂断块速度模型精度的方法
CN104570116A (zh) 基于地质标志层的时差分析校正方法
CN104199088B (zh) 一种提取入射角道集的方法及系统
CN104597485B (zh) 一种微小断层检测方法及断层检测装置
EA030770B1 (ru) Система и способ адаптивной сейсмической оптики
CN107340537A (zh) 一种p-sv转换波叠前逆时深度偏移的方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant