WO2023124922A1 - 角度域广义Radon变换的离散方法、系统、终端及介质 - Google Patents

角度域广义Radon变换的离散方法、系统、终端及介质 Download PDF

Info

Publication number
WO2023124922A1
WO2023124922A1 PCT/CN2022/138195 CN2022138195W WO2023124922A1 WO 2023124922 A1 WO2023124922 A1 WO 2023124922A1 CN 2022138195 W CN2022138195 W CN 2022138195W WO 2023124922 A1 WO2023124922 A1 WO 2023124922A1
Authority
WO
WIPO (PCT)
Prior art keywords
discrete
angle
scattering angle
interval
domain
Prior art date
Application number
PCT/CN2022/138195
Other languages
English (en)
French (fr)
Inventor
栗学磊
魏彦杰
冯圣中
Original Assignee
深圳先进技术研究院
中国科学院深圳理工大学(筹)
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 深圳先进技术研究院, 中国科学院深圳理工大学(筹) filed Critical 深圳先进技术研究院
Publication of WO2023124922A1 publication Critical patent/WO2023124922A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the present application belongs to the technical field of angle domain information extraction, and in particular relates to a discrete method, system, terminal and medium of angle domain generalized Radon transform.
  • the inventor proposed and established the angle-domain generalized Radon transform (AD-GRT) inversion method (Li Xuelei et al., 2020), which can realize multi-parameter inversion of acoustic waves and elastic media without iteration.
  • AD-GRT angle-domain generalized Radon transform
  • the theoretical framework of the generalized Radon transform in the angle domain is:
  • the wave field of the acoustic wave medium can be expressed as a generalized Radon transform (GRT) form (frequency domain):
  • p′(r, s, ⁇ ) is the acoustic pressure scattering wave field
  • a s and ⁇ s respectively represent the amplitude and travel time of the source end
  • a r and ⁇ r are the amplitude and travel time of the receiving end respectively
  • F( ⁇ ) is the source wavelet.
  • ⁇ s ⁇ s (x) is the direction angle from s to the x-ray at x
  • f(x, ⁇ ) is an angle domain model
  • ⁇ 0 (x) and ⁇ 0 (x) are the background models of density reciprocal and compressibility, respectively, and ⁇ ' and ⁇ ' are the corresponding disturbance parameters.
  • the GRT integral transformation is equivalent to the integral transformation from the angle domain model f(x, ⁇ ) to the acoustic pressure scattering wavefield p′(r,s,t), where s and r are only distributed on the boundary line of the 2D disturbance region.
  • both f(x, ⁇ ) and p'(r, s, t) have a three-dimensional distribution, therefore, the corresponding inverse transformation from p'(r, s, t) to f(x, ⁇ ) Also buildable.
  • AD-GRT Angle-Domain Generalized Radon Inverse Transform
  • J s (y) and J r (y) are the Jacobians, and can be expressed as and ⁇ s (y) and ⁇ r (y) are the direction angles of the mappable source s and receiver point r.
  • J s (y) and J r (y) can be represented by geometric diffusion from y to s and r:
  • ⁇ is the sampling interval of scattering angle ⁇ 0.
  • the integral term of formula (4) is accumulated to the ⁇ 0 sampling point of f(y, ⁇ 0 ) , and finally multiplied by 1/ ⁇ to average.
  • this conventional algorithm will lead to very unbalanced calculation of f(y, ⁇ 0 ) with different ⁇ 0 , and it is prone to saw-tooth fluctuations distributed with ⁇ 0 .
  • the present application provides a discrete method, system, terminal and medium of angle-domain generalized Radon transform, aiming to solve one of the above-mentioned technical problems in the prior art at least to a certain extent.
  • a discretization method of generalized Radon transform in angle domain comprising:
  • Calculate the nearest sampling point ⁇ 0 of the scattering angle and its partition boundary ⁇ is the scattering angle interval, and calculate the distance l from the scattering angle to the partition boundary ⁇ bound , judge whether the partition boundary ⁇ bound splits the scattering angle discrete surface according to the distance l, if so, use the discrete element splitting algorithm to divide the partition
  • the discrete unit near the boundary ⁇ bound is split into two units, and the area of the split two units is calculated according to the discrete interval;
  • the two split units are used to calculate the imaging accumulation at the sampling point ⁇ 0 and ⁇ 0 ⁇ respectively.
  • the technical solution adopted in the embodiment of the present application also includes: before reading the seismic trace data, it also includes:
  • the shot point coordinates and the corresponding ray travel time, amplitude and direction table are set, and the receiving point coordinates and the corresponding ray travel time, amplitude and direction table are set.
  • the technical solution adopted in the embodiment of the present application also includes: the travel time, amplitude and direction table of reading the coordinates of each imaging point are specifically:
  • the technical solution adopted in the embodiment of the present application also includes: the calculation of the scattering angle and the discrete interval according to the travel time, amplitude and direction table corresponding to the coordinates of each imaging point is specifically:
  • the scattering angle ⁇ mid ⁇ s - ⁇ r ;
  • the technical solution adopted in the embodiment of the present application further includes: the discrete unit splitting algorithm includes a rectangular splitting algorithm and a parallel splitting algorithm.
  • the rectangle splitting algorithm is specifically:
  • the parallel splitting algorithm is specifically:
  • a discrete system of angle-domain generalized Radon transform including:
  • Data reading module used to read the seismic trace data, set the imaging point coordinates of the seismic trace data, and read the travel time, amplitude and direction table of each imaging point coordinate;
  • Scattering angle calculation module used to calculate the scattering angle and discrete interval according to the travel time, amplitude and direction table corresponding to the coordinates of each imaging point;
  • Discrete unit splitting module used to calculate the nearest sampling point ⁇ 0 of the scattering angle and its partition boundary ⁇ is the scattering angle interval, and calculate the distance l from the scattering angle to the partition boundary ⁇ bound , judge whether the partition boundary ⁇ bound splits the scattering angle discrete surface according to the distance l, if so, use the discrete element splitting algorithm to divide the partition
  • the discrete unit near the boundary ⁇ bound is split into two units, and the area of the split two units is calculated according to the discrete interval;
  • Imaging module used to calculate the imaging accumulation at the sampling point ⁇ 0 and ⁇ 0 ⁇ using the two split units respectively.
  • a terminal includes a processor and a memory coupled to the processor, wherein,
  • the memory is stored with program instructions for realizing the discrete method of the generalized Radon transform in the angle domain;
  • the processor is configured to execute the program instructions stored in the memory to control the discretization of the generalized Radon transform in the angle domain.
  • a storage medium storing program instructions executable by a processor, and the program instructions are used to execute the discrete method of the generalized Radon transform in the angle domain.
  • the beneficial effect produced by the embodiment of the present application lies in that the discrete method, system, terminal and storage medium of the angle-domain generalized Radon transform in the embodiment of the present application are designed by combining the theoretical framework and physical meaning of AD-GRT transform
  • Two angle-domain discrete difference algorithms conforming to the AD-GRT transformation theory by quantitatively splitting the discrete unit area, the discrete discrete unit is reasonably split into two units, which realizes the accurate calculation of each discrete unit and ensures that the algorithm is efficient and convenient
  • the smooth continuity of the angle domain amplitude is ensured, and the amplitude oscillation problem of the conventional angle domain discrete algorithm is solved.
  • the embodiment of the present application can obtain continuous and smooth angle domain inversion results, and realize efficient and stable angle domain information extraction.
  • Fig. 1 is the flowchart of the discrete method of the angle domain generalized Radon transform of the embodiment of the present application
  • FIG. 2 is a schematic diagram of the distribution of rectangular discrete units and scattering angle contours in an embodiment of the present application
  • FIG. 3 is a schematic diagram of the distribution of parallelogram discrete units and scattering angle contours in an embodiment of the present application
  • Figure 4a is a schematic diagram of a horizontal single-interface model
  • Figure 4(b) is a schematic diagram of a synthetic single-shot seismic record
  • FIG. 7 is a schematic structural diagram of a discrete system of an angle-domain generalized Radon transform according to an embodiment of the present application.
  • FIG. 8 is a schematic structural diagram of a terminal according to an embodiment of the present application.
  • FIG. 9 is a schematic structural diagram of a storage medium according to an embodiment of the present application.
  • FIG. 1 is a flow chart of the discrete method of the angle-domain generalized Radon transform according to the embodiment of the present application.
  • the discrete method of the angle domain generalized Radon transform of the embodiment of the present application comprises the following steps:
  • S1 Set the operating environment parameters such as shot point discrete interval ⁇ s, receiving point discrete interval ⁇ r, and scattering angle interval ⁇ ;
  • S2 Set the coordinates of the shot point and the corresponding ray travel time, amplitude and direction table, and set the coordinates of the receiving point and the corresponding ray travel time, amplitude and direction table;
  • the embodiment of the present application further analyzes the integral summation related to the scattering angle ⁇ , and designs a reasonable numerical solution scheme.
  • the discretized micro-element satisfies the relational expression:
  • ⁇ s and ⁇ r are the discrete intervals of the shot point and the receiving point, respectively, and ⁇ s and ⁇ r are the ray direction intervals of rays s and r propagating to y, respectively, which can be obtained by calculating the Jacobian coefficients J s and J r . Therefore, a discrete sampling point of (s, r) corresponds to a discrete unit in y, and the area of the discrete unit is ⁇ s ⁇ r . When the sampling point is very close to the partition boundary ⁇ 0 ⁇ /2, some of the discrete units should belong to the partition near the sampling point ⁇ 0 . Based on the above, in the embodiment of the present application, the discrete unit near the partition boundary ⁇ bound is split into two units, and the area sizes S 1 and S 2 of the two split units are calculated.
  • FIG. 2 it is a schematic diagram of the distribution of rectangular discrete units and scattering angle contours in the embodiment of the present application.
  • the size of the discrete unit is split into two parts, S 1 and S 2 , and instead of the J s J r ⁇ s ⁇ r part in formula (4), they are accumulated to ⁇ 0 and its In the inversion result of f(y, ⁇ 0 ) nearby.
  • S 2 is:
  • Figure 3 it is a schematic diagram of the distribution of parallelogram discrete units and scattering angle contours .
  • This method can also control the discrete unit independently through ⁇ s , which can be realized by exchanging all d r and d s in the formula.
  • the parallel splitting algorithm is simple and convenient to calculate, and is suitable for future high-dimensional discrete calculations.
  • S8 Use two units S 1 and S 2 to calculate the imaging accumulation at the sampling point ⁇ 0 ⁇ and the scattering angle sampling point ⁇ 0 respectively;
  • S10 Determine whether there are remaining receiving points, if there are remaining receiving points, set the next receiving point, and re-execute S2 to S9; otherwise, execute S11;
  • S11 Determine whether there are remaining shot points, if there are remaining shot points, set the next shot point, and re-execute S2 to S10; otherwise, end.
  • there is obvious oblique noise at the edge of the inversion results (see area A), which is related to the data truncation at the boundary of the maximum offset of the seismic record. Record boundaries are characterized by data discontinuities that do not meet the local approximation assumptions of the provided AD-GRT.
  • the experimental results show that the traditional discretization method tends to form obvious saw-tooth oscillations in the angle domain information extraction, especially in the small angle range.
  • Both the rectangular splitting algorithm and the parallel splitting algorithm proposed in the embodiment of the present application can provide a very continuous and smooth amplitude distribution effect, and the inversion amplitude distribution and the true value distribution of the model can be highly fitted within the effective angle range.
  • the discretization method of the generalized Radon transform in the angle domain of the embodiment of the present application combines the theoretical framework and physical meaning of the AD-GRT transform, and designs two discretization difference algorithms in the angle domain that conform to the AD-GRT transform theory.
  • Unit area the discretized discrete unit is reasonably split into two units, which realizes the accurate calculation of each discrete unit, ensures the efficiency and convenience of the algorithm, and ensures the smooth continuity of the angle domain amplitude, which solves the problem of conventional angle domain discretization
  • the amplitude oscillation problem of the algorithm The embodiment of the present application can obtain continuous and smooth angle domain inversion results, and realize efficient and stable angle domain information extraction.
  • FIG. 7 is a schematic structural diagram of a discrete system of the angle-domain generalized Radon transform according to an embodiment of the present application.
  • the discrete system 40 of the angle-domain generalized Radon transform of the embodiment of the present application includes:
  • Data reading module 41 used to read the seismic trace data, set the imaging point coordinates of the seismic trace data, and read the travel time, amplitude and direction table of each imaging point coordinate;
  • Scattering angle calculation module 42 used to calculate the scattering angle and discrete interval according to the travel time, amplitude and direction table corresponding to the coordinates of each imaging point;
  • Discrete unit splitting module 43 used to calculate the sampling point ⁇ 0 with the closest scattering angle and its partition boundary ⁇ is the scattering angle interval, and the distance l from the scattering angle to the partition boundary ⁇ bound is calculated. According to the distance l , it is judged whether the partition boundary ⁇ bound splits the scattering angle discrete surface. The unit is split into two units, and the area of the split two units is calculated according to the discrete interval;
  • the imaging module 44 is used to calculate the imaging accumulation at the sampling point ⁇ 0 and ⁇ 0 ⁇ using the two split units respectively.
  • FIG. 8 is a schematic diagram of a terminal structure according to an embodiment of the present application.
  • the terminal 50 includes a processor 51 and a memory 52 coupled to the processor 51 .
  • the memory 52 stores program instructions for realizing the discrete method of the above-mentioned angle-domain generalized Radon transform.
  • the processor 51 is used to execute the program instructions stored in the memory 52 to control the discretization of the generalized Radon transform in the angle domain.
  • the processor 51 may also be referred to as a CPU (Central Processing Unit, central processing unit).
  • the processor 51 may be an integrated circuit chip with signal processing capabilities.
  • the processor 51 can also be a general-purpose processor, a digital signal processor (DSP), an application-specific integrated circuit (ASIC), an off-the-shelf programmable gate array (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components .
  • a general-purpose processor may be a microprocessor, or the processor may be any conventional processor, or the like.
  • FIG. 9 is a schematic structural diagram of a storage medium according to an embodiment of the present application.
  • the storage medium of the embodiment of the present application stores a program file 61 capable of realizing all the above-mentioned methods, wherein the program file 61 can be stored in the above-mentioned storage medium in the form of a software product, and includes several instructions to make a computer device (which can It is a personal computer, a server, or a network device, etc.) or a processor (processor) that executes all or part of the steps of the methods in various embodiments of the present invention.
  • a computer device which can It is a personal computer, a server, or a network device, etc.
  • processor processor
  • the aforementioned storage media include: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), magnetic disk or optical disc, etc., which can store program codes. , or terminal devices such as computers, servers, mobile phones, and tablets.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Image Processing (AREA)

Abstract

本申请涉及一种角度域广义Radon变换的离散方法、系统、终端及介质。所述方法包括：读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；计算所述散射角最近的采样点θ 0及其分区边界 (I) Δθ为散射角间隔，并计算所述散射角到分区边界θ bound的距离l，根据所述距离l判断分区边界θ bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；分别使用分裂后的两个单元计算采样点θ 0以及θ 0±Δθ处的成像累加。本申请可以获得连续平滑的角度域反演结果。

Description

角度域广义Radon变换的离散方法、系统、终端及介质 技术领域
本申请属于角度域信息提取技术领域，特别涉及一种角度域广义Radon变换的离散方法、系统、终端及介质。
背景技术
发明人提出且建立了角度域广义Radon变换(AD-GRT)反演方法(栗学磊等,2020)，在非迭代情况下即可实现声波和弹性介质等多参数反演。其中，角度域广义Radon变换理论框架为：
基于Born近似，声波介质波场可表示为一种广义Radon变换(GRT)形式(频率域)：
p′(r,s,ω)＝ω 2∫dx 1dx 3A sA rexp[iω(φ s+φ r)]κ 0F(ω)f(x,θ)       (1)
其中，p′(r,s,ω)为声波压力散射波场，A s和φ s分别代表震源端振幅和走时，A r和φ r分别为接收端振幅和走时，F(ω)为震源子波。θ＝θ(s,x,r)为散射角，是s到x射线和r到x射线在x处的夹角，且满足：
θ＝α s-α r         (2)
其中，α s＝α s(x)是从s到x射线在x处的方向角，α r＝α r(x)是从r到x射线在x处的方向角。
f(x,θ)为一种角度域模型：
Figure PCTCN2022138195-appb-000001
σ 0(x)和κ 0(x)分别为密度倒数和压缩率的背景模型，σ′和κ′为对应的扰动参数。
GRT积分变换相当于角度域模型f(x,θ)到声波压力散射波场p′(r,s,t)的积分变换，此处的s和r仅分布在2D扰动区域的边界线上。理论上，f(x,θ)和p′(r,s,t)均具有三个维度分布，因此，对应的从p′(r,s,t)到f(x,θ)的逆变换也可搭建。通过建立对应的角度域广 义Radon逆变换(AD-GRT)，以支持从p′(r,s,t)到f(x,θ)的反向积分变换，反向积分变换表达式如下：
Figure PCTCN2022138195-appb-000002
其中，J s(y)和J r(y)为雅可比行列式，且可表示为
Figure PCTCN2022138195-appb-000003
和
Figure PCTCN2022138195-appb-000004
α s(y)和α r(y)为可映射震源s和接收点r的方向角。J s(y)和J r(y)可由从y到s和r的几何扩散表示：
Figure PCTCN2022138195-appb-000005
但是式(4)中的δ(θ-θ 0)并不能直接参与计算，需要利用其他的近似插值或求和算法代替。并且，θ 0离散采样点与震源s和接收点r的离散分布并不能直接对应，给AD-GRT的精确数值计算带来一定麻烦。为此，需要在常规角度域离散算法的基础上，设计AD-GRT的准确离散算法，反向积分变换表达式(4)可以近似表示为：
Figure PCTCN2022138195-appb-000006
其中，Δθ为散射角θ 0的采样间隔，在传统角度域离散算法中，如果
Figure PCTCN2022138195-appb-000007
在成像点y处对应的散射角θ属于[θ 0-Δθ/2,θ 0+Δθ/2)区间，则式(4)的积分项累加到f(y,θ 0)的θ 0采样点上，最后乘以1/Δθ来平均。然而，由于离散的整数化特征，这种常规算法会导致不同θ 0的f(y,θ 0)计算很不均衡，易出现随θ 0分布的锯齿状震荡起伏。尤其是震源点和接收点分布比较稀疏，或者θ 0取值采样比较密集时，这种震荡起伏非常严重，导致难以提取准确的f(y,θ 0)分布关系，致使角度域反演结果无效。即使一些平滑因子可以减缓的震荡，但依然存在明显起伏，并且平滑因子会带来额外的信息损失。
发明内容
本申请提供了一种角度域广义Radon变换的离散方法、系统、终端及介质，旨在至少在一定程度上解决现有技术中的上述技术问题之一。
为了解决上述问题，本申请提供了如下技术方案：
一种角度域广义Radon变换的离散方法，包括：
读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；
根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；
计算所述散射角最近的采样点θ 0及其分区边界
Figure PCTCN2022138195-appb-000008
Δθ为散射角间隔，并计算所述散射角到分区边界θ bound的距离l，根据所述距离l判断分区边界θ bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；
分别使用分裂后的两个单元计算采样点θ 0以及θ 0±Δθ处的成像累加。
本申请实施例采取的技术方案还包括：所述读取地震道数据之前还包括：
设定炮点离散间隔Δs、接收点离散间隔Δr以及散射角间隔Δθ；
设定所述炮点坐标和对应的射线走时、幅值以及方向表，并设定所述接收点坐标和对应的射线走时、幅值以及方向表。
本申请实施例采取的技术方案还包括：所述读取每一个成像点坐标的走时、幅值以及方向表具体为：
所述方向表中包括方向角α s和α r，其中，α s＝α s(x)为从s到x射线在x处的方向角，α r＝α r(x)为从r到x射线在x处的方向角。
本申请实施例采取的技术方案还包括：所述根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔具体为：
所述散射角θ mid＝α s-α r；
所述离散间隔d s＝J sΔs和d r＝J rΔr，其中J s和J r为雅可比行列式，Δs为炮点离散间隔，Δr为接收点离散间隔。
本申请实施例采取的技术方案还包括：所述离散单元分裂算法包括矩形分裂算法和平行分裂算法。
本申请实施例采取的技术方案还包括：所述矩形分裂算法具体为：
设d s＝Δα s，d r＝Δα r，l＝|θ bound-θ mid|，Δα s和Δα r分别为射线s和r传播到y处的射线方向间隔，Δα s≈J sΔs,Δα r≈J rΔr；
如果
Figure PCTCN2022138195-appb-000009
则表示分区边界θ bound穿过离散单元，则求解分裂后的两个单元S 1和S 2：
如果d s≥d r：
Figure PCTCN2022138195-appb-000010
如果d s<d r：
Figure PCTCN2022138195-appb-000011
S 2＝d sd r-S 1.
本申请实施例采取的技术方案还包括：所述平行分裂算法具体为：
如果
Figure PCTCN2022138195-appb-000012
表示分区边界θ bound穿过离散单元，则求解S 1和S 2：
Figure PCTCN2022138195-appb-000013
Figure PCTCN2022138195-appb-000014
S 1+S 2＝d sd r。
本申请实施例采取的另一技术方案为：一种角度域广义Radon变换的离散系统，包括：
数据读取模块：用于读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；
散射角计算模块：用于根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；
离散单元分裂模块：用于计算所述散射角最近的采样点θ 0及其分区边界
Figure PCTCN2022138195-appb-000015
Figure PCTCN2022138195-appb-000016
Δθ为散射角间隔，并计算所述散射角到分区边界θ bound的距离l，根据所述距离l判断分区边界θ bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；
成像模块：用于分别使用分裂后的两个单元计算采样点θ 0以及θ 0±Δθ处的成像累加。
本申请实施例采取的又一技术方案为：一种终端，所述终端包括处理器、与所述处理器耦接的存储器，其中，
所述存储器存储有用于实现所述角度域广义Radon变换的离散方法的程序指令；
所述处理器用于执行所述存储器存储的所述程序指令以控制角度域广义Radon变换的离散。
本申请实施例采取的又一技术方案为：一种存储介质，存储有处理器可运行的程序指令，所述程序指令用于执行所述角度域广义Radon变换的离散方法。
相对于现有技术，本申请实施例产生的有益效果在于：本申请实施例的角度域广义Radon变换的离散方法、系统、终端以及存储介质通过结合AD-GRT变换理论框架和物理意义，设计了两种符合AD-GRT变换理论的角度域离散差值算法，通过定量分裂离散单元面积，将离散后的离散单元合理分裂为两个单元，实现了每一个离散单元的精确计算，确保算法高效便捷的同时，确保了角度域幅值的平滑连续性，解决了常规角度域离散算法的幅值震荡问题。本申请实施例可以获得连续平滑的角度域反演结果，实现高效稳定的角度域信息提取。
附图说明
图1是本申请实施例的角度域广义Radon变换的离散方法的流程图；
图2为本申请实施例的矩形离散单元与散射角等值线分布示意图；
图3为本申请实施例的平行四边形离散单元与散射角等值线分布示意图；
图4a为水平单界面模型示意图，图4(b)为合成单炮地震记录示意图；
图5为利用AD-GRT反演方法求解的角度域模型f(x,θ)反演结果(x＝2000m)示意图；
图6为在水平线z＝1000m上的角度域反演赋值分布示意图，其中(a)、(b)、(c)分别为使用传统离散方法、矩形分裂算法、平行分裂算法的真实值和反演值；
图7为本申请实施例的角度域广义Radon变换的离散系统结构示意图；
图8为本申请实施例的终端结构示意图；
图9为本申请实施例的存储介质的结构示意图。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本申请，并不用于限定本申请。
请参阅图1，是本申请实施例的角度域广义Radon变换的离散方法的流程图。本申请实施例的角度域广义Radon变换的离散方法包括以下步骤：
S1：设定炮点离散间隔Δs、接收点离散间隔Δr以及散射角间隔Δθ等运行环境参数；
S2：设定炮点坐标和对应的射线走时、幅值以及方向表，并设定接收点坐标和对应的射线走时、幅值以及方向表；
其中，方向表中包括方向角α s和α r，其中，α s＝α s(x)是从s到x射线在x处的方向角，α r＝α r(x)是从r到x射线在x处的方向角。
S3：读取地震道数据；
S4：设定地震道数据的成像点坐标，并读取每一个成像点坐标对应的走时、幅值以及方向表；
S5：根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角θ mid＝α s-α r及其离散间隔d s＝J sΔs和d r＝J rΔr，其中J s和J r为雅可比行列式；
S6：计算散射角θ mid最近的散射角采样点θ 0及其分区边界
Figure PCTCN2022138195-appb-000017
并计算散射角θ mid到分区边界θ bound的距离l＝|θ bound-θ mid|；
S7：根据距离l＝|θ bound-θ nid|判断分区边界θ bound是否分裂散射角θ mid离散面，如果是，利用离散单元分裂算法将分区边界θ bound附近的离散单元分裂为两个单元S 1和S 2，并计算分裂后的两个单元的面积大小；
本步骤中，为了解决传统角度域算法的离散问题，本申请实施例对散射角θ相关的积分求和做进一步分析，设计合理的数值求解方案。离散后的微元满足关系式：
Δα s≈J sΔs,Δα r≈J rΔr          (7)
其中Δs和Δr分别为炮点和接收点的离散间隔，Δα s和Δα r分别为射线s和r传播到y的射线方向间隔，可通过雅克比系数J s和J r计算求得。因此，一个(s,r)的离散采样点在y对应的是一个离散单元，离散单元的面积大小为Δα s×Δα r。当采样点非常靠近分区边界θ 0±Δθ/2时，该离散单元中有一部分应属于采样点θ 0附近分区。基于上述，本申请实施例将分区边界θ bound附近的离散单元分裂为两个单元，并求解分裂后的两个单元的面积大小S 1和S 2。
具体的，如图2所示，为本申请实施例的矩形离散单元与散射角等值线分布示意图。矩形分裂算法具体为：每一个(α s,α r)采样点的离散单元由一个矩形表示，矩形面积大小为Δα s×Δα r，采样点位于长和宽的中间位置，θ bound＝θ 0±Δθ/2为θ 0分区边界，θ mid为(α s,α r)采样点对应的θ值。当分区边界θ bound穿过一个离散单元时，将该离散单元大小拆分为S 1和S 2两部分，并代替式(4)中的J sJ rΔsΔr部分，分别累加到θ 0及其附近的f(y,θ 0)反演结果中。
设d s＝Δα s，d r＝Δα r，l＝|θ bound-θ mid|；如果
Figure PCTCN2022138195-appb-000018
则分区边界θ bound穿过该离散单元，需要求解分裂后的两个单元S 1和S 2，求解方法如下：
如果d s≥d r：
Figure PCTCN2022138195-appb-000019
如果d s<d r：
Figure PCTCN2022138195-appb-000020
基于S 1，S 2为：
S 2＝d sd r-S 1.          (10)
平行分裂算法具体为：散射角θ满足
Figure PCTCN2022138195-appb-000021
如果设置Δα s＝0(或α s为常数)，则式(2)中的θ积分范围可通过α r单独控制.如图3所示，为平行四边形离散单元与散射角等值线分布示意图。通过将(α s,α r)采样点的离散单元由一个平行四边形表示(即，其中两条边在α s常数线上，另外两条边平行于θ等值线)，此时的离散单元分裂只与α r离散情况有关。
如果
Figure PCTCN2022138195-appb-000022
则θ bound穿过该离散单元，需要求解S 1和S 2，并且S 1和S 2满足：
Figure PCTCN2022138195-appb-000023
Figure PCTCN2022138195-appb-000024
此处同样满足S 1+S 2＝d sd r.该方法也可以通过α s单独控制离散单元，调换公式中的所有d r和d s即可实现。平行分裂算法计算简洁方便，且适用于未来高维离散计算。
S8：分别使用两个单元S 1和S 2计算采样点θ 0±Δθ以及散射角采样点θ 0处的成像累加；
S9：判断地震道数据的成像点是否循环结束，如果没有循环结束，则获取下一个成像点，并重新执行S4至S8；否则，执行S10；
S10：判断是否存在剩余接收点，如果存在剩余接收点，则设置下一个接收点，并重新执行S2至S9；否则执行S11；
S11：判断是否存在剩余炮点，如果存在剩余炮点，则设置下一个炮点，并重新执行S2至S10；否则，结束。
数据实验:
为了验证本申请实施例的可行性和有效性，以下利用简单模型数据来验证本申请实施例的计算效果。具体为：水平单界面模型如图4a所示，水平界面在深度z＝1000m，并且界面上下声波波速分别为2000m/s和2100m/s，依照Gardner关系式ρ＝0.31×c0.25(波速单位:m/s,密度单位:g/cm3)设定模型密度。模型在水平和垂直方向的网格间隔均为5m，炮点和检波点均分布于边界线z＝0m上，且炮点间隔和检波点间隔均为10m，从炮点到检波点的最大偏移距为2000m。该模型地震记录使用声波波动方程有限差分方法合成，以验证准确波场数据反演结果的有效性。用于AD-GRT的背景波速和密度模型通过平滑原始模型来获得。图4(b)为合成单炮地震记录示意图，显示了震源x＝2000m的单炮记录，这里只显示了主要反射信息，排除了直达波信息。
图5给出了利用AD-GRT反演方法求解的角度域模型f(x,θ)反演结果(x＝2000m)。为排除滤波导致的幅值震荡，此处设置震源子波为F(ω)＝1，忽略子波反滤波操作。反演的主要信息分布在水平线z＝1000m附近，且同相轴具有明显的水平分布特征。同时，在反演结果边缘位置存在明显的倾斜噪音(见区域A)，这与地震记录最大偏移距边界处的数据截断有关。记录边界具有数据不连续特征，而不连续性不符合提供的AD-GRT的局部近似假设。
为了验证本申请实施例的离散单元分裂法在数值计算方面的有效性，以下分别利用传统离散方法、矩形分裂算法和平行分裂算法的数值计算反演结果，具体如图6所示，为在水平线z＝1000m上的角度域反演赋值分布示意图，其中(a)、(b)、(c)分别为使用传 统离散方法、矩形分裂算法、平行分裂算法的真实值和反演值。图6的横轴方向使用角度的余弦值表示，以显示角度域模型的线性倾斜分布特征。其中的真实值是界面z＝1000m上下扰动值的差值(跳跃值)，并且为方便对比，反演值乘了一个常数。实验结果表明，传统离散方法在角度域信息提取中易形成明显的锯齿状震荡，尤其在小角度范围内。而本申请实施例提出的矩形分裂算法和平行分裂算法均能提供非常连续平滑的幅值分布效果，并且反演幅值分布与模型真实值分布在有效角度范围内均能高度拟合。
基于上述，本申请实施例的角度域广义Radon变换的离散方法通过结合AD-GRT变换理论框架和物理意义，设计了两种符合AD-GRT变换理论的角度域离散差值算法，通过定量分裂离散单元面积，将离散后的离散单元合理分裂为两个单元，实现了每一个离散单元的精确计算，确保算法高效便捷的同时，确保了角度域幅值的平滑连续性，解决了常规角度域离散算法的幅值震荡问题。本申请实施例可以获得连续平滑的角度域反演结果，实现高效稳定的角度域信息提取。
请参阅图7，为本申请实施例的角度域广义Radon变换的离散系统结构示意图。本申请实施例的角度域广义Radon变换的离散系统40包括：
数据读取模块41：用于读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；
散射角计算模块42：用于根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；
离散单元分裂模块43：用于计算散射角最近的采样点θ 0及其分区边界
Figure PCTCN2022138195-appb-000025
Figure PCTCN2022138195-appb-000026
Δθ为散射角间隔，并计算散射角到分区边界θ bound的距离l，根据距离l判断分区边界θ bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将分区边界θ bound附近的离散单元分裂为两个单元，并根据离散间隔计算分裂后的两个单元的面积；
成像模块44：用于分别使用分裂后的两个单元计算采样点θ 0以及θ 0±Δθ处的成像累加。
请参阅图8，为本申请实施例的终端结构示意图。该终端50包括处理器51、与处理器51耦接的存储器52。
存储器52存储有用于实现上述角度域广义Radon变换的离散方法的程序指令。
处理器51用于执行存储器52存储的程序指令以控制角度域广义Radon变换的离散。
其中，处理器51还可以称为CPU(Central Processing Unit，中央处理单元)。处理器51可能是一种集成电路芯片，具有信号的处理能力。处理器51还可以是通用处理器、数字信号处理器(DSP)、专用集成电路(ASIC)、现成可编程门阵列(FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
请参阅图9，为本申请实施例的存储介质的结构示意图。本申请实施例的存储介质存储有能够实现上述所有方法的程序文件61，其中，该程序文件61可以以软件产品的形式存储在上述存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)或处理器(processor)执行本发明各个实施方式方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质，或者是计算机、服务器、手机、平板等终端设备。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本申请。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本申请中所定义的一般原理可以在不脱离本申请的精神或范围的情况下，在其它实施例中实现。因此，本申请将不会被限制于本申请所示的这些实施例，而是要符合与本申请所公开的原理和新颖特点相一致的最宽的范围。

Claims (10)

  1. 一种角度域广义Radon变换的离散方法，其特征在于，包括：
    读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；
    根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；
    计算所述散射角最近的采样点θ 0及其分区边界
    Figure PCTCN2022138195-appb-100001
    Δθ为散射角间隔，并计算所述散射角到分区边界θ bound的距离l，根据所述距离l判断分区边界θ bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；
    分别使用分裂后的两个单元计算采样点θ 0以及θ 0±Δθ处的成像累加。
  2. 根据权利要求1所述的角度域广义Radon变换的离散方法，其特征在于，所述读取地震道数据之前还包括：
    设定炮点离散间隔Δs、接收点离散间隔Δr以及散射角间隔Δθ；
    设定所述炮点坐标和对应的射线走时、幅值以及方向表，并设定所述接收点坐标和对应的射线走时、幅值以及方向表。
  3. 根据权利要求2所述的角度域广义Radon变换的离散方法，其特征在于，所述读取每一个成像点坐标的走时、幅值以及方向表具体为：
    所述方向表中包括方向角α s和α r，其中，α s＝α s(x)为从s到x射线在x处的方向角，α r＝α r(x)为从r到x射线在x处的方向角。
  4. 根据权利要求3所述的角度域广义Radon变换的离散方法，其特征在于，所述根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔具体为：
    所述散射角θ mid＝α s-α r；
    所述离散间隔d s＝J sΔs和d r＝J rΔr，其中J s和J r为雅可比行列式，Δs为炮点离散间隔，Δr为接收点离散间隔。
  5. 根据权利要求1至4任一项所述的角度域广义Radon变换的离散方法，其特征在于，所述离散单元分裂算法包括矩形分裂算法和平行分裂算法。
  6. 根据权利要求5所述的角度域广义Radon变换的离散方法，其特征在于，所述矩形分裂算法具体为：
    设d s＝Δα s，d r＝Δα r，l＝|θ bound-θ mid|，Δα s和Δα r分别为射线s和r传播到y处的射线方向间隔，Δα s≈J sΔs,Δα r≈J rΔr；
    如果
    Figure PCTCN2022138195-appb-100002
    则表示分区边界θ bound穿过离散单元，则求解分裂后的两个单元S 1和S 2：
    如果d s≥d r：
    Figure PCTCN2022138195-appb-100003
    如果d s<d r：
    Figure PCTCN2022138195-appb-100004
    S 2＝d sd r-S 1.
  7. 根据权利要求5所述的角度域广义Radon变换的离散方法，其特征在于，所述平行分裂算法具体为：
    如果
    Figure PCTCN2022138195-appb-100005
    表示分区边界θ bound穿过离散单元，则求解S 1和S 2：
    Figure PCTCN2022138195-appb-100006
    Figure PCTCN2022138195-appb-100007
    S 1+S 2＝d sd r。
  8. 一种角度域广义Radon变换的离散系统，其特征在于，包括：
    数据读取模块：用于读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；
    散射角计算模块：用于根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；
    离散单元分裂模块：用于计算所述散射角最近的采样点θ 0及其分区边界
    Figure PCTCN2022138195-appb-100008
    Figure PCTCN2022138195-appb-100009
    Δθ为散射角间隔，并计算所述散射角到分区边界θ bound的距离l，根据所述距离l判断分区边界θ bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；
    成像模块：用于分别使用分裂后的两个单元计算采样点θ 0以及θ 0±Δθ处的成像累加。
  9. 一种终端，其特征在于，所述终端包括处理器、与所述处理器耦接的存储器，其中，
    所述存储器存储有用于实现权利要求1-7任一项所述的角度域广义Radon变换的离散方法的程序指令；
    所述处理器用于执行所述存储器存储的所述程序指令以控制角度域广义Radon变换的离散。
  10. 一种存储介质，其特征在于，存储有处理器可运行的程序指令，所述程序指令用于执行权利要求1至7任一项所述角度域广义Radon变换的离散方法。
PCT/CN2022/138195 2021-12-30 2022-12-09 角度域广义Radon变换的离散方法、系统、终端及介质 WO2023124922A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202111660036.9A CN114417660A (zh) 2021-12-30 2021-12-30 角度域广义Radon变换的离散方法、系统、终端及介质
CN202111660036.9 2021-12-30

Publications (1)

Publication Number Publication Date
WO2023124922A1 true WO2023124922A1 (zh) 2023-07-06

Family

ID=81272075

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2022/138195 WO2023124922A1 (zh) 2021-12-30 2022-12-09 角度域广义Radon变换的离散方法、系统、终端及介质

Country Status (2)

Country Link
CN (1) CN114417660A (zh)
WO (1) WO2023124922A1 (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114417660A (zh) * 2021-12-30 2022-04-29 深圳先进技术研究院 角度域广义Radon变换的离散方法、系统、终端及介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070064530A1 (en) * 2005-08-26 2007-03-22 Ian Moore Method for processing a record of seismic traces
US20120004849A1 (en) * 2010-03-22 2012-01-05 Schlumberger Technology Corporation Efficient windowed radon transform
CN106772583A (zh) * 2017-01-10 2017-05-31 中国科学院地质与地球物理研究所 一种地震绕射波分离方法和装置
CN108415073A (zh) * 2018-03-06 2018-08-17 中国科学院测量与地球物理研究所 角度域逆散射偏移成像方法及装置
CN112698400A (zh) * 2020-12-04 2021-04-23 中国科学院深圳先进技术研究院 反演方法、反演装置、计算机设备和计算机可读存储介质
CN114417660A (zh) * 2021-12-30 2022-04-29 深圳先进技术研究院 角度域广义Radon变换的离散方法、系统、终端及介质

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070064530A1 (en) * 2005-08-26 2007-03-22 Ian Moore Method for processing a record of seismic traces
US20120004849A1 (en) * 2010-03-22 2012-01-05 Schlumberger Technology Corporation Efficient windowed radon transform
CN106772583A (zh) * 2017-01-10 2017-05-31 中国科学院地质与地球物理研究所 一种地震绕射波分离方法和装置
CN108415073A (zh) * 2018-03-06 2018-08-17 中国科学院测量与地球物理研究所 角度域逆散射偏移成像方法及装置
CN112698400A (zh) * 2020-12-04 2021-04-23 中国科学院深圳先进技术研究院 反演方法、反演装置、计算机设备和计算机可读存储介质
CN114417660A (zh) * 2021-12-30 2022-04-29 深圳先进技术研究院 角度域广义Radon变换的离散方法、系统、终端及介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
TEYFOURI NILOUFAR, RABBANI HOSSEIN, KAFIEH RAHELEH, JABBARI IRAJ: "An Exact and Fast CBCT Reconstruction via Pseudo-Polar Fourier Transform-Based Discrete Grangeat’s Formula", IEEE TRANSACTIONS ON IMAGE PROCESSING, IEEE, USA, vol. 29, 1 January 2020 (2020-01-01), USA, pages 5832 - 5847, XP093075620, ISSN: 1057-7149, DOI: 10.1109/TIP.2020.2985874 *
ZHANG, CHUAN-QIANG; LI, LING-YUN: "Multiple Suppression Based on the Improved Radon Transform in the Frequency Domain", PETROLEUM GEOPHYSICS, vol. 13, no. 1, 26 January 2015 (2015-01-26), pages 18 - 22, XP009553127 *

Also Published As

Publication number Publication date
CN114417660A (zh) 2022-04-29

Similar Documents

Publication Publication Date Title
WO2023124922A1 (zh) 角度域广义Radon变换的离散方法、系统、终端及介质
Tagliasacchi et al. Mean curvature skeletons
Wu et al. A highly accurate finite-difference method with minimum dispersion error for solving the Helmholtz equation
Trefethen et al. Exponential node clustering at singularities for rational approximation, quadrature, and PDEs
Salman et al. Feature preserving mesh generation from 3d point clouds
WO2019242045A1 (zh) 一种虚拟震源二维波前构建地震波走时计算方法
Javidrad et al. Contour curve reconstruction from cloud data for rapid prototyping
TW201610730A (zh) 點雲網格簡化系統及方法
Boehm et al. Time-domain spectral-element ultrasound waveform tomography using a stochastic quasi-Newton method
Hendin et al. Tsunami and acoustic-gravity waves in water of constant depth
Yang et al. Development of 3D PUFEM with linear tetrahedral elements for the simulation of acoustic waves in enclosed cavities
Marburg et al. Cat's eye radiation with boundary elements: Comparative study on treatment of irregular frequencies
Kim Extraction of ridge and valley lines from unorganized points
Fan et al. Modeling, analysis, and validation of an active T-shaped noise barrier
Liu et al. Variational progressive-iterative approximation for RBF-based surface reconstruction
Chernokozhin et al. Scattering by thin shells in fluids: Fast solver and experimental validation
Gorthi et al. Windowed high-order ambiguity function method for fringe analysis
Bruno et al. Efficient evaluation of doubly periodic Green functions in 3D scattering, including Wood anomaly frequencies
Salomons et al. Sound propagation in a turbulent atmosphere near the ground: An approach based on the spectral representation of refractive-index fluctuations
Aristégui et al. New results for isotropic point scatterers: Foldy revisited
Lee et al. Computation of scattering of a plane wave from multiple prolate spheroids using the collocation multipole method
Marcon et al. Fast PDE approach to surface reconstruction from large cloud of points
Blackmore et al. Spatial filters suppress ripple artifacts in the computation of acoustic fields with the angular spectrum method
Antuono et al. Validation of a three-dimensional depth-semi-averaged model
Mehrabi et al. On estimating the curvature attributes and strain invariants of deformed surface through radial basis functions

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 22914177

Country of ref document: EP

Kind code of ref document: A1