点扩散函数快速空间褶积方法与流程

1.本发明涉及石油地球物理勘探地震数据处理技术领域，特别是涉及到一种点扩散函数快速空间褶积方法。


背景技术：

2.作为油气勘探的关键技术环节之一，偏移成像技术在地震勘探理论研究和实际应用中具有重要的地位，在日趋复杂和精细化的地下介质成像中发挥着越来越重要的作用。然而，传统的偏移成像方法只是born正演算子的共轭转置，对于有限观测系统所采集的带限地震数据，只能产生模糊的构造成像结果。此外，复杂的地下介质构造和非规则的地震数据空间采样，还会导致偏移假象和非均匀的成像照明，严重影响成像振幅的可靠性。
3.基于线性化地震反演理论的最小二乘偏移技术是解决上述问题的有效途径。该技术将地震偏移成像视为线性反问题，并利用最优化方法进行迭代求解，理论上能够消除非规则采集、带限子波等因素对成像结果的不利影响，提高成像分辨率和照明度。然而，该技术的理论优势并没有转化成为实际的应用效果，其主要原因在于计算成本过于高昂。经典的数据域最小二乘偏移每次迭代均需一次常规正演和偏移运算，完整的反演过程往往需要几十倍于常规偏移的计算量，因此庞大的计算成本严重制约了该方法的实用性。成像域最小二乘偏移则需要计算和存储庞大的hessian矩阵，即便可以通过数据编码等技术优化相应的计算过程，其计算和存储的开销依然昂贵，难以应用于大规模三维实际地震数据处理。
4.在申请号：cn202010476294.0的中国专利申请中，涉及到一种深度域成像模拟方法及系统。该深度域成像模拟方法包括：将初始密度模型划分为窗口，在各个窗口中设置密度扰动点；根据密度扰动点影响下的初始密度模型确定扰动点密度模型；根据扰动点密度模型和初始密度模型得到点散射波场数据；对点散射波场数据进行偏移处理，得到点扩散函数；根据点扩散函数和反射率模型得到深度域模拟成像。该发明可以得到与实际偏移成像匹配且稳定连续的模拟成像。
5.在申请号：cn201210290287.7的中国专利申请中，涉及到一种用于图像复原的基于kalman模型的点扩散函数估计方法。该发明首先建立点扩散函数估计的对称全平面kalman状态方程和观测方程；其次根据选择图像平稳区域，确定点扩散函数模型的观测矩阵；最后根据平稳区域进行kalman迭代，估计出图像的点扩散函数。该发明有效改善了点扩散函数估计不准确和非自适应的缺点，可以根据图像自适应的估计点扩散函数，估计的点扩散函数数值更准确。为后续图像复原提供了更准确的点扩散函数。
6.在申请号：cn201910427576.9的中国专利申请中，涉及到一种点扩散函数重建方法，包括：(1)输入目标点扩散函数；(2)针对输入的目标点扩散函数进行采样；(3)将采样得到的每个点看作一个偶极子进行初始相位和初始振幅估计；(4)计算得到偶极子在瞳面上的辐射场；(5)通过辐射场得到瞳面上的振幅和相位分布；(6)通过得到的振幅和相位生成点扩散函数；(7)将点扩散函数与目标点扩散函数比较，若满足要求，则输出得到的振幅和相位分布；若不满足要求，则更新初始相位和初始振幅，返回步骤(4)。该发明直接由目标点
扩散函数得到瞳面上的相位及振幅分布，大大提高了计算效率。并且，由于偶极子辐射的光是向随机方向传播的，该方法可以扩展到多个具有多个光瞳面的显微系统。
7.以上现有技术均与本发明有较大区别，未能解决我们想要解决的技术问题，为此我们发明了一种新的点扩散函数快速空间褶积方法。


技术实现要素：

8.本发明的目的是提供一种提高成像域最小二乘偏移的计算效率的点扩散函数快速空间褶积方法。
9.本发明的目的可通过如下技术措施来实现：点扩散函数快速空间褶积方法，该点扩散函数快速空间褶积方法包括：
10.步骤1，计算成像剖面和点扩散函数场；
11.步骤2，对点扩散函数进行循环，确定其中心位置的空间位置坐标；
12.步骤3，确定插值方法，并计算针对点扩散函数有贡献的局部成像剖面的插值系数；
13.步骤4，进行波数域快速空间褶积；
14.步骤5，滤波后将剖面乘以相应的插值系数；
15.步骤6，累加局部成像剖面的滤波结果，得到最终的滤波剖面。
16.本发明的目的还可通过如下技术措施来实现：
17.在步骤1中，利用地震偏移和点扩散函数计算求取偏移成像剖面和位于稀疏空间位置的点扩散函数场。
18.在步骤2中，开始计算循环，读取每个点扩散函数，并确定其中心点的空间位置。
19.在步骤3中，根据选择的插值方法，确定该点扩散函数有贡献的局部成像范围，以及对应每个成像样点的插值系数。
20.在步骤4中，利用快速傅里叶变换，计算点扩散函数和局部成像剖面波数谱的乘积。
21.在步骤4中，计算点扩散函数和局部成像剖面波数谱的乘积的公式为：
22.m(k)＝h(k)r(k)
23.其中，h(k)，r(k)分别为点扩散函数和局部成像剖面的波数谱，m(k)为点扩散函数滤波后局部成像剖面波数谱。
24.在步骤5中，将波数谱反变换到空间域得到空间褶积滤波后的成像剖面，将其乘以相应的插值系数。
25.在步骤5中，滤波后将剖面乘以相应的插值系数的公式为：
26.m(x)＝cof(x)*ifft(m(k))
27.其中，m(x)是应用插值系数后的空间域滤波成像剖面，cof(x)为空间位置x处的插值系数，ifft代表反傅里叶变换。
28.在步骤6中，重复上述步骤(2)-步骤(5)，并将褶积滤波后的成像剖面进行累加，即可得到最终的点扩散函数褶积滤波结果。
29.在步骤6中，计算点扩散函数空间滤波成像剖面的公式为：
30.m
′
(x)＝m
′
(x)+m(x)
31.其中，m
′
(x)为最终的点扩散函数空间滤波成像剖面；m(x)是应用插值系数后的空间域滤波成像剖面。
32.本发明中的点扩散函数快速空间褶积方法，首先利用fft进行快速的空间褶积滤波，然后应用相应的插值系数后将滤波结果进行累加，得到同空间褶积滤波相同的计算结果。本发明提供了一种快速的点扩散函数空间褶积滤波方法，基于插值方法线性累加的特点，采用fft进行点扩散函数同局部成像剖面的褶积，再应用插值系数进行计算求和得到褶积滤波结果。相对于传统方法，该方法可以在保证计算精度的同时，大幅提升点扩散函数空间褶积滤波过程的计算效率，降低计算成本。
附图说明
33.图1为本发明的点扩散函数快速空间褶积方法的一具体实施例的流程图；
34.图2为本发明的一具体实施例中marmousi模型对应的反射率剖面的示意图；
35.图3为本发明的一具体实施例中位于marmousi模型稀疏空间位置的点扩散函数场的示意图；
36.图4为本发明的一具体实施例中利用传统的点扩散函数空间褶积滤波的结果的示意图；
37.图5为本发明的一具体实施例中利用本发明所给出的快速点扩散函数空间褶积滤波的结果的示意图。
具体实施方式
38.应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
39.需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作和/或它们的组合。
40.本发明的点扩散函数快速空间褶积方法包括以下步骤：
41.(1)利用地震偏移和点扩散函数计算程序求取偏移成像剖面和位于稀疏空间位置的点扩散函数场；
42.(2)对点扩散函数进行循环，确定其中心位置的空间位置坐标；
43.(3)根据选择的插值方法，确定该点扩散函数有贡献的局部成像范围，以及对应每个成像样点的插值系数；
44.(4)利用快速傅里叶变换(fft)，计算点扩散函数和局部成像剖面波数谱的乘积；
45.(5)将波数谱反变换到空间域得到空间褶积滤波后的成像剖面，将其乘以相应的插值系数；
46.(6)重复上述步骤(2)-步骤(5)，并将褶积滤波后的成像剖面进行累加，即可得到最终的点扩散函数褶积滤波结果。
47.以下为应用本发明的几个具体实施例
48.实施例1
49.在应用本发明的一具体实施例1中，图1是本发明所描述的点扩散函数空间褶积滤波方法的一个实施例的流程图，具体实施步骤包括：
50.(1)计算成像剖面和点扩散函数场
51.利用地震偏移和点扩散函数计算程序求取偏移成像剖面和位于稀疏空间位置的点扩散函数场；
52.(2)开始计算循环，读取每个点扩散函数，并确定其中心点的空间位置
53.对点扩散函数进行循环，确定其中心位置的空间位置坐标；
54.(3)确定插值方法，并计算针对点扩散函数有贡献的局部成像剖面的插值系数
55.根据选择的插值方法，确定该点扩散函数有贡献的局部成像范围，以及对应每个成像样点的插值系数；
56.(4)进行波数域快速空间褶积
57.利用快速傅里叶变换(fft)，计算点扩散函数和局部成像剖面波数谱的乘积；
58.(5)滤波后剖面乘以相应的插值系数
59.将波数谱反变换到空间域得到空间褶积滤波后的成像剖面，将其乘以相应的插值系数；
60.(6)累加局部成像剖面的滤波结果，得到最终的滤波剖面
61.重复上述步骤(2)-步骤(5)，并将褶积滤波后的成像剖面进行累加，即可得到最终的点扩散函数褶积滤波结果。
62.实施例2：
63.在应用本发明的具体实施例2中，图2-图5为本发明所描述的点扩散函数空间褶积滤波方法的应用效果：
64.1)利用marmousi模型速度场计算对应的反射率剖面，结果如图2所示；
65.2)利用点扩散函数计算方法，计算marmousi模型位于稀疏空间位置的点扩散函数，结果如图3所示；
66.3)利用常规的点扩散函数空间褶积方法，求取点扩散函数同marmousi模型反射率剖面的褶积滤波剖面，结果如图4所示，耗时约3.6s；
67.4)利用本发明求取点扩散函数同marmousi模型反射率的褶积滤波剖面，结果如图5所示，耗费时间约0.5s；
68.通过对比发现，可以大幅度提升点扩散函数空间褶积滤波过程的计算效率，本发明计算成本相对于传统方法减小7倍以上。
69.实施例3：
70.在应用本发明的具体实施例3中，本发明所描述的点扩散函数空间褶积滤波方法具体包括了以下步骤：
71.(1)计算成像剖面和点扩散函数场
72.利用地震偏移和点扩散函数计算程序求取偏移成像剖面和位于稀疏空间位置的点扩散函数场；成像剖面和点扩散函数的计算利用常规算法实现即可。
73.(2)开始计算循环，读取每个点扩散函数，并确定其中心点的空间位置
74.对点扩散函数进行循环，确定其中心位置的空间位置坐标；
75.(3)确定插值方法，并计算针对点扩散函数有贡献的局部成像剖面的插值系数
76.根据选择的插值方法，确定该点扩散函数有贡献的局部成像范围，以及对应每个成像样点的插值系数；插值方法的选取取决于所需要的精度，但是一般来说常规的线性插值即可满足精度要求。
77.(4)进行波数域快速空间褶积
78.利用快速傅里叶变换(fft)，计算点扩散函数和局部成像剖面波数谱的乘积；
79.m(k)＝h(k)r(k)
80.其中，h(k)，r(k)分别为点扩散函数和局部成像剖面的波数谱，m(k)为点扩散函数滤波后局部成像剖面波数谱。
81.(5)滤波后剖面乘以相应的插值系数
82.将波数谱反变换到空间域得到空间褶积滤波后的成像剖面，将其乘以相应的插值系数；
83.m(x)＝cof(x)*ifft(m(k))
84.其中，m(x)是应用插值系数后的空间域滤波成像剖面，cof(x)为空间位置x处的插值系数，ifft代表反傅里叶变换。
85.(6)累加局部成像剖面的滤波结果，得到最终的滤波剖面
86.重复上述步骤(2)-步骤(5)，并将褶积滤波后的成像剖面进行累加，即可得到最终的点扩散函数褶积滤波结果。
87.m
′
(x)＝m
′
(x)+m(x)
88.其中，m
′
(x)为最终的点扩散函数空间滤波成像剖面。
89.最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域技术人员来说，其依然可以对前述实施例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
90.除说明书所述的技术特征外，均为本专业技术人员的已知技术。

技术特征：
1.点扩散函数快速空间褶积方法，其特征在于，该点扩散函数快速空间褶积方法包括：步骤1，计算成像剖面和点扩散函数场；步骤2，对点扩散函数进行循环，确定其中心位置的空间位置坐标；步骤3，确定插值方法，并计算针对点扩散函数有贡献的局部成像剖面的插值系数；步骤4，进行波数域快速空间褶积；步骤5，滤波后将剖面乘以相应的插值系数；步骤6，累加局部成像剖面的滤波结果，得到最终的滤波剖面。2.根据权利要求1所述的点扩散函数快速空间褶积方法，其特征在于，在步骤1中，利用地震偏移和点扩散函数计算求取偏移成像剖面和位于稀疏空间位置的点扩散函数场。3.根据权利要求1所述的点扩散函数快速空间褶积方法，其特征在于，在步骤2中，开始计算循环，读取每个点扩散函数，并确定其中心点的空间位置。4.根据权利要求1所述的点扩散函数快速空间褶积方法，其特征在于，在步骤3中，根据选择的插值方法，确定该点扩散函数有贡献的局部成像范围，以及对应每个成像样点的插值系数。5.根据权利要求1所述的点扩散函数快速空间褶积方法，其特征在于，在步骤4中，利用快速傅里叶变换，计算点扩散函数和局部成像剖面波数谱的乘积。6.根据权利要求5所述的点扩散函数快速空间褶积方法，其特征在于，在步骤4中，计算点扩散函数和局部成像剖面波数谱的乘积的公式为：m(k)＝h(k)r(k)其中，h(k)，r(k)分别为点扩散函数和局部成像剖面的波数谱，m(k)为点扩散函数滤波后局部成像剖面波数谱。7.根据权利要求1所述的点扩散函数快速空间褶积方法，其特征在于，在步骤5中，将波数谱反变换到空间域得到空间褶积滤波后的成像剖面，将其乘以相应的插值系数。8.根据权利要求7所述的点扩散函数快速空间褶积方法，其特征在于，在步骤5中，滤波后将剖面乘以相应的插值系数的公式为：m(x)＝cof(x)*ifft(m(k))其中，m(x)是应用插值系数后的空间域滤波成像剖面，cof(x)为空间位置x处的插值系数，ifft代表反傅里叶变换。9.根据权利要求1所述的点扩散函数快速空间褶积方法，其特征在于，在步骤6中，重复上述步骤(2)-步骤(5)，并将褶积滤波后的成像剖面进行累加，即可得到最终的点扩散函数褶积滤波结果。10.根据权利要求9所述的点扩散函数快速空间褶积方法，其特征在于，在步骤6中，计算点扩散函数空间滤波成像剖面的公式为：m
′
(x)＝m
′
(x)+m(x)其中，m
′
(x)为最终的点扩散函数空间滤波成像剖面；m(x)是应用插值系数后的空间域滤波成像剖面。

技术总结
本发明提供一种点扩散函数快速空间褶积方法，该点扩散函数快速空间褶积方法包括：步骤1，计算成像剖面和点扩散函数场；步骤2，对点扩散函数进行循环，确定其中心位置的空间位置坐标；步骤3，确定插值方法，并计算针对点扩散函数有贡献的局部成像剖面的插值系数；步骤4，进行波数域快速空间褶积；步骤5，滤波后将剖面乘以相应的插值系数；步骤6，累加局部成像剖面的滤波结果，得到最终的滤波剖面。该点扩散函数快速空间褶积方法基于插值方法线性累加的特点，可以大幅度提升点扩散函数空间褶积滤波过程的计算效率，降低计算成本。降低计算成本。降低计算成本。


技术研发人员：秦宁 毕丽飞 梁鸿贤 杨晓东 曹晓莉 谢伟 孙兴刚
受保护的技术使用者：中国石油化工股份有限公司胜利油田分公司物探研究院
技术研发日：2022.03.23
技术公布日：2023/10/7