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

Abstract

本申请涉及一种角度域广义Radon变换的离散方法、系统、终端以及存储介质。所述方法包括：读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；计算散射角及离散间隔；计算所述散射角最近的采样点θ₀及其分区边界

并计算所述散射角到分区边界θ_bound的距离l，根据所述距离l判断分区边界θ_bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ_bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；分别使用分裂后的两个单元计算采样点θ₀以及θ₀±Δθ处的成像累加。本申请可以获得连续平滑的角度域反演结果。

Description

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

技术领域

本申请属于角度域信息提取技术领域，特别涉及一种角度域广义Radon变换的离散方法、系统、终端以及存储介质。

背景技术

发明人提出且建立了角度域广义Radon变换(AD-GRT)反演方法(栗学磊等,2020)，在非迭代情况下即可实现声波和弹性介质等多参数反演。其中，角度域广义Radon变换理论框架为：

基于Born近似，声波介质波场可表示为一种广义Radon变换(GRT)形式(频率域)：

p′(r,s,ω)＝ω²∫dx₁dx₃A_sA_rexp[iω(φ_s+φ_r)]κ⁰F(ω)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,θ)为一种角度域模型：

σ⁰(x)和κ⁰(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,θ)的反向积分变换，反向积分变换表达式如下：

其中，J_s(y)和J_r(y)为雅可比行列式，且可表示为

和

α_s(y)和α_r(y)为可映射震源s和接收点r的方向角。J_s(y)和J_r(y)可由从y到s和r的几何扩散表示：

但是式(4)中的δ(θ-θ₀)并不能直接参与计算，需要利用其他的近似插值或求和算法代替。并且，θ₀离散采样点与震源s和接收点r的离散分布并不能直接对应，给AD-GRT的精确数值计算带来一定麻烦。为此，需要在常规角度域离散算法的基础上，设计AD-GRT的准确离散算法，反向积分变换表达式(4)可以近似表示为：

其中，Δθ为散射角θ₀的采样间隔，在传统角度域离散算法中，如果

在成像点y处对应的散射角θ属于[θ₀-Δθ/2,θ₀+Δθ/2)区间，则式(4)的积分项累加到f(y,θ₀)的θ₀采样点上，最后乘以1/Δθ来平均。然而，由于离散的整数化特征，这种常规算法会导致不同θ₀的f(y,θ₀)计算很不均衡，易出现随θ₀分布的锯齿状震荡起伏。尤其是震源点和接收点分布比较稀疏，或者θ₀取值采样比较密集时，这种震荡起伏非常严重，导致难以提取准确的f(y,θ₀)分布关系，致使角度域反演结果无效。即使一些平滑因子可以减缓的震荡，但依然存在明显起伏，并且平滑因子会带来额外的信息损失。

发明内容

本申请提供了一种角度域广义Radon变换的离散方法、系统、终端以及存储介质，旨在至少在一定程度上解决现有技术中的上述技术问题之一。

为了解决上述问题，本申请提供了如下技术方案：

一种角度域广义Radon变换的离散方法，包括：

读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；

根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；

计算所述散射角最近的采样点θ₀及其分区边界

Δθ为散射角间隔，并计算所述散射角到分区边界θ_bound的距离l，根据所述距离l判断分区边界θ_bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ_bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；

分别使用分裂后的两个单元计算采样点θ₀以及θ₀±Δθ处的成像累加。

本申请实施例采取的技术方案还包括：所述读取地震道数据之前还包括：

设定炮点离散间隔Δ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；

如果

则表示分区边界θ_bound穿过离散单元，则求解分裂后的两个单元S₁和S₂：

如果d_s≥d_r：

如果d_s＜d_r：

S₂＝d_sd_r-S₁

本申请实施例采取的技术方案还包括：所述平行分裂算法具体为：

如果

表示分区边界θ_bound穿过离散单元，则求解S₁和S₂：

S₁+S₂＝d_sd_r。

本申请实施例采取的另一技术方案为：一种角度域广义Radon变换的离散系统，包括：

数据读取模块：用于读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；

散射角计算模块：用于根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；

离散单元分裂模块：用于计算所述散射角最近的采样点θ₀及其分区边界

Δθ为散射角间隔，并计算所述散射角到分区边界θ_bound的距离l，根据所述距离l判断分区边界θ_bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ_bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；

成像模块：用于分别使用分裂后的两个单元计算采样点θ₀以及θ₀±Δθ处的成像累加。

本申请实施例采取的又一技术方案为：一种终端，所述终端包括处理器、与所述处理器耦接的存储器，其中，

所述存储器存储有用于实现所述角度域广义Radon变换的离散方法的程序指令；

所述处理器用于执行所述存储器存储的所述程序指令以控制角度域广义Radon变换的离散。

本申请实施例采取的又一技术方案为：一种存储介质，存储有处理器可运行的程序指令，所述程序指令用于执行所述角度域广义Radon变换的离散方法。

相对于现有技术，本申请实施例产生的有益效果在于：本申请实施例的角度域广义Radon变换的离散方法、系统、终端以及存储介质通过结合AD-GRT变换理论框架和物理意义，设计了两种符合AD-GRT变换理论的角度域离散差值算法，通过定量分裂离散单元面积，将离散后的离散单元合理分裂为两个单元，实现了每一个离散单元的精确计算，确保算法高效便捷的同时，确保了角度域幅值的平滑连续性，解决了常规角度域离散算法的幅值震荡问题。本申请实施例可以获得连续平滑的角度域反演结果，实现高效稳定的角度域信息提取。

附图说明

图1是本申请实施例的角度域广义Radon变换的离散方法的流程图；

图2为本申请实施例的矩形离散单元与散射角等值线分布示意图；

图3为本申请实施例的平行四边形离散单元与散射角等值线分布示意图；

图4a为水平单界面模型示意图，图4(b)为合成单炮地震记录示意图；

图5为利用AD-GRT反演方法求解的角度域模型f(x,θ)反演结果示意图；

图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最近的散射角采样点θ₀及其分区边界

并计算散射角θ_mid到分区边界θ_bound的距离l＝|θ_bound-θ_mid|；

S7：根据距离l＝|θ_bound-θ_mid|判断分区边界θ_bound是否分裂散射角θ_mid离散面，如果是，利用离散单元分裂算法将分区边界θ_bound附近的离散单元分裂为两个单元S₁和S₂，并计算分裂后的两个单元的面积大小；

本步骤中，为了解决传统角度域算法的离散问题，本申请实施例对散射角θ相关的积分求和做进一步分析，设计合理的数值求解方案。离散后的微元满足关系式：

Δα_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。当采样点非常靠近分区边界θ₀±Δθ/2时，该离散单元中有一部分应属于采样点θ₀附近分区。基于上述，本申请实施例将分区边界θ_bound附近的离散单元分裂为两个单元，并求解分裂后的两个单元的面积大小S₁和S₂。

具体的，如图2所示，为本申请实施例的矩形离散单元与散射角等值线分布示意图。矩形分裂算法具体为：每一个(α_s,α_r)采样点的离散单元由一个矩形表示，矩形面积大小为Δα_s×Δα_r，采样点位于长和宽的中间位置，θ_bo_und＝θ₀±Δθ/2为θ₀分区边界，θ_mid为(α_s,α_r)采样点对应的θ值。当分区边界θ_bound穿过一个离散单元时，将该离散单元大小拆分为S₁和S₂两部分，并代替式(4)中的J_sJ_rΔsΔr部分，分别累加到θ₀及其附近的f(y,θ₀)反演结果中。

设d_s＝Δα_s，d_r＝Δα_r，l＝|θ_bound-θ_mid|；如果

则分区边界θ_bound穿过该离散单元，需要求解分裂后的两个单元S₁和S₂，求解方法如下：

如果d_s≥d_r：

如果d_s＜d_r：

基于S₁，S₂为：

S₂＝d_sd_r-S₁. (10)

平行分裂算法具体为：散射角θ满足

如果设置Δα_s＝0(或α_s为常数)，则式(2)中的θ积分范围可通过α_r单独控制.如图3所示，为平行四边形离散单元与散射角等值线分布示意图。通过将(α_s,α_r)采样点的离散单元由一个平行四边形表示(即，其中两条边在α_s常数线上，另外两条边平行于θ等值线)，此时的离散单元分裂只与α_r离散情况有关。

如果

则θ_bound穿过该离散单元，需要求解S₁和S₂，并且S₁和S₂满足：

此处同样满足S₁+S₂＝d_sd_r.该方法也可以通过α_s单独控制离散单元，调换公式中的所有d_r和d_s即可实现。平行分裂算法计算简洁方便，且适用于未来高维离散计算。

S8：分别使用两个单元S₁和S₂计算采样点θ₀±Δθ以及散射角采样点θ₀处的成像累加；

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：用于计算散射角最近的采样点θ₀及其分区边界

Δθ为散射角间隔，并计算散射角到分区边界θ_bound的距离l，根据距离l判断分区边界θ_bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将分区边界θ_bound附近的离散单元分裂为两个单元，并根据离散间隔计算分裂后的两个单元的面积；

成像模块44：用于分别使用分裂后的两个单元计算采样点θ₀以及θ₀±Δθ处的成像累加。

请参阅图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变换的离散方法，其特征在于，包括：

读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；

根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；

计算所述散射角最近的采样点θ₀及其分区边界

Δθ为散射角间隔，并计算所述散射角到分区边界θ_bound的距离l，根据所述距离l判断分区边界θ_bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ_bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；

分别使用分裂后的两个单元计算采样点θ₀以及θ₀±Δθ处的成像累加。

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；

如果

则表示分区边界θ_bound穿过离散单元，则求解分裂后的两个单元S₁和S₂：

如果d_s≥d_r：

如果d_s＜d_r：

S₂＝d_sd_r-S₁。

7.根据权利要求5所述的角度域广义Radon变换的离散方法，其特征在于，所述平行分裂算法具体为：

如果

表示分区边界θ_bound穿过离散单元，则求解S₁和S₂：

S₁+S₂＝d_sd_r。

8.一种角度域广义Radon变换的离散系统，其特征在于，包括：

数据读取模块：用于读取地震道数据，设定地震道数据的成像点坐标，并读取每一个成像点坐标的走时、幅值以及方向表；

散射角计算模块：用于根据每一个成像点坐标对应的走时、幅值以及方向表计算散射角及离散间隔；

离散单元分裂模块：用于计算所述散射角最近的采样点θ₀及其分区边界

Δθ为散射角间隔，并计算所述散射角到分区边界θ_bound的距离l，根据所述距离l判断分区边界θ_bound是否分裂散射角离散面，如果是，利用离散单元分裂算法将所述分区边界θ_bound附近的离散单元分裂为两个单元，并根据所述离散间隔计算分裂后的两个单元的面积；

成像模块：用于分别使用分裂后的两个单元计算采样点θ₀以及θ₀±Δθ处的成像累加。

9.一种终端，其特征在于，所述终端包括处理器、与所述处理器耦接的存储器，其中，

所述存储器存储有用于实现权利要求1-7任一项所述的角度域广义Radon变换的离散方法的程序指令；

所述处理器用于执行所述存储器存储的所述程序指令以控制角度域广义Radon变换的离散。

10.一种存储介质，其特征在于，存储有处理器可运行的程序指令，所述程序指令用于执行权利要求1至7任一项所述角度域广义Radon变换的离散方法。

Images