CN113805233A - 一种点扩散函数的计算方法 - Google Patents
一种点扩散函数的计算方法 Download PDFInfo
- Publication number
- CN113805233A CN113805233A CN202010536572.7A CN202010536572A CN113805233A CN 113805233 A CN113805233 A CN 113805233A CN 202010536572 A CN202010536572 A CN 202010536572A CN 113805233 A CN113805233 A CN 113805233A
- Authority
- CN
- China
- Prior art keywords
- point
- calculating
- spatial
- spread function
- coordinates
- 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
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 27
- 239000011159 matrix material Substances 0.000 claims abstract description 37
- 238000000034 method Methods 0.000 claims abstract description 34
- 239000000523 sample Substances 0.000 claims abstract description 11
- 238000006073 displacement reaction Methods 0.000 claims description 15
- 238000003384 imaging method Methods 0.000 claims description 14
- 238000001514 detection method Methods 0.000 claims description 9
- 230000005284 excitation Effects 0.000 claims description 6
- 230000005855 radiation Effects 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 3
- 238000004613 tight binding model Methods 0.000 claims description 3
- 229920013655 poly(bisphenol-A sulfone) Polymers 0.000 description 8
- 238000007796 conventional method Methods 0.000 description 5
- 238000013508 migration Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 230000005012 migration Effects 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000013506 data mapping Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
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
本发明提供一种点扩散函数的计算方法，包括以下步骤：(1)读取速度模型和观测系统，选取任意炮点和检波点并分别确定其空间坐标；(2)计算所选取的炮点和检波点对地下空间的射线追踪；(3)通过给定震源函数，对地下目标空间位置进行循环，根据下述公式在每个目标空间位置的局部邻域内计算；其中，所述Hessian矩阵的每个单元即为点扩散函数。本发明通过高频近似，获取了PSF的解析表达式，提出一种不需要偏移，也不需要插值，可以计算任意空间点位PSF的方法。本方法在保证计算精度的同时，大幅提高了PSF的计算效率。
Description
技术领域
本发明涉及地震资料处理领域，具体涉及地震资料偏移成像领域。
背景技术
点扩散函数(Point Scatter Function，PSF)的物理意义是地下单一绕射点的地震成像响应，是构成地震成像数据的基础，求取点扩散函数可以进行通过地下模型直接获取成像结果地震正演，也可以插值求取Hessian矩阵(黑塞矩阵)进而进行反演，从而求取反射系数。
现有的地震数据和地下模型之间的关系可以表达为如下方程：
d＝Lm，
其中，d是观测数据，m是地下模型，通常指反射系数，L是模型到数据的映射矩阵。在地震勘探中，这个矩阵代表波场传播效应。通过这个公式可以看出，地震数据和地下模型之间的关系是一个大型的线性系统。
地震偏移成像采用与这个过程类似的过程，实际上也就是对数据作用上L算子的转置：
I＝LTd＝LTLm＝Hm，
其中，I是偏移获得的像，H成为Hessian矩阵。通过上述公式可以得出，成像结果就是地下介质作用上Hessian矩阵。Hessian矩阵的意义是地下每一个离散点存在单位激励时产生的偏移成像响应，因此它的每一行成为点扩散函数。
Hessian矩阵是沟通偏移成像结果和地下介质的桥梁，而PSF是Hessian矩阵的单元，计算出所有PSF也就计算出了Hessian矩阵。目前计算PSF的方法是通过反偏移和偏移算法实现的，具体做法如下：
(1)读取速度模型和观测系统，设计分布在成像空间的稀疏的离散点，离散点的最小间距不能使PSF在计算的过程中发生重叠；
(2)计算正演过程(反偏移)
y＝Lx，
其中，x是上一步构建的离散点阵，y是离散点阵的正演数据；
(3)计算偏移过程
P＝LTy，
其中，P是离散点阵的PSF；
(4)如果需要离散点以外的PSF，使用空间插值算法利用临近的已知点进行空间插值。
这个流程中采用反偏移-偏移算法计算离散点的PSF，偏移本身是一个计算量极大的处理环节，反偏移与偏移计算量相当，计算一次离散点的PSF需要两次深度偏移的计算量。得到的离散PSF考虑到不能重叠，所以稀疏性往往是成像空间的几千甚至几万分之一，导致后期的插值误差较大，无法保证计算精度。
发明内容
本发明的目的在于提供一种点扩散函数的计算方法，在保证计算精度的同时，大幅提高PSF的计算效率。
为达上述目的，本发明提供一种点扩散函数的计算方法，其包括以下步骤：
(1)读取速度模型和观测系统，选取任意炮点和检波点并分别确定其空间坐标；
(2)计算所选取的炮点和检波点对地下空间的射线追踪；
(3)通过给定震源函数，对地下目标空间位置进行循环，根据下述公式在每个目标空间位置的局部邻域内计算：
其中，H代表Hessian矩阵，S代表震源函数，
表示炮点空间坐标，
表示检波点空间坐标，
表示离散激励的存在坐标，
是表示炮点空间坐标
出发的射线在离散激励的存在坐标
位置的方向矢量，
是表示检波点空间坐标
出发的射线在离散激励的存在坐标
位置的方向矢，
是空间位移矢量，
是炮点空间位移矢量，
是检波点空间位移矢量；
其中，所述Hessian矩阵的每个单元即为点扩散函数。
在优选实施例中，在步骤(1)中，建立Hessian矩阵的表达式：
其中，H代表Hessian矩阵，G是Green函数，uI表示入射波场，
表示炮点空间坐标，
表示检波点空间坐标，
表示空间坐标，
表示离散激励的存在坐标，t0和t0′表示时间，
是炮点空间位移矢量，
是检波点空间位移矢量。
在优选实施例中，在步骤(1)中，对所建立的Hessian矩阵进行高频近似处理：
其中，S代表震源函数，A代表振幅场，τ代表走时场，δ表示狄拉克抽样函数，带入Hessian表达式，得到高频近似条件下的Hessian矩阵的表达式：
在优选实施例中，在步骤(1)中，使
进一步近似的Hessian矩阵为：
在优选实施例中，在步骤(2)中，将进一步近似的Hessian矩阵带入射线追踪系统：
其中，
是表示炮点空间坐标
出发的射线在空间坐标
位置的方向矢量，
是空间位移矢量，得出Hessian矩阵为：
在优选实施例中，在步骤(3)中，忽略振幅项，得出Hessian矩阵为：
在优选实施例中，在步骤(3)中，所述地下目标空间位置包括全部成像空间。
在优选实施例中，在步骤(1)中，采用3000m/s进行点扩散函数计算。
在优选实施例中，在步骤(1)中，空间分布为起始点在(200m，200m)的位置。
在优选实施例中，在步骤(1)中，x方向的增量为400m。
本发明的有益效果在于：本发明通过高频近似，获取了PSF的解析表达式，提出一种不需要偏移，也不需要插值，可以计算任意空间点位PSF的方法。本方法在保证计算精度的同时，大幅提高了PSF的计算效率。
附图说明
图1为现有技术中反偏移-偏移方法计算的稀疏空间位置PSF；
图2为根据本发明的计算方法计算的稀疏空间位置PSF；
图3为空间坐标(8000m，1200m)处常规方法计算的单个PSF(a)与本发明方法计算的单个PSF(b)的对比图；
图4为空间坐标(8000m，2000m)处常规方法计算的单个PSF(a)与本发明方法计算的单个PSF(b)的对比图。
图5为空间坐标(8000m，4000m)处常规方法计算的单个PSF(a)与本发明方法计算的单个PSF(b)的对比图。
图6为现有方法和本发明方法计算效率的对比图。
具体实施方式
下面将结合附图对本发明作进一步说明。
首先，读取速度模型和观测系统，建立Hessian矩阵的表达式：
其中，H代表Hessian矩阵，G是Green函数，uI表示入射波场，
表示炮点空间坐标，
表示检波点空间坐标，
表示空间坐标，
表示离散激励的存在坐标，t0和t0′表示时间，
是炮点空间位移矢量，
是检波点空间位移矢量。
应用高频近似：
其中，S代表震源函数，A代表振幅场，τ代表走时场，δ表示狄拉克抽样函数，带入Hessian表达式，可以得到高频近似条件下的Hessian矩阵的表达式：
对于计算PSF，考虑到其分布半径较小，因此
可以进一步近似：
将上述方程带入射线追踪系统：
其中，
是表示炮点空间坐标
出发的射线在空间坐标
位置的方向矢量，
是表示检波点空间坐标
出发的射线在空间坐标
位置的方向矢量，
是空间位移矢量，可得出：
其中，
对于能量分布范围局限在一个小邻域内的PSF，振幅可以看成一个系统差，代表球面扩散的效应，由于相对于PSF而言相位更加重要，因此忽略振幅项，得出：
从而获取了一个PSF的解析表达式，通过这个表达式计算PSF，也就是求取一个特定空间位置
的Hessian矩阵，仅需要已知震源函数，并求取空间波传播的方向矢量。
根据上述公式，快速计算任意位置的PSF基本步骤为：
(1)读取速度模型和观测系统，计算任意炮点和检波点对地下空间的射线追踪，求取每个空间位置的射线方向，也就是波的传播方向P；
(2)给定震源函数，对地下目标空间位置(也可以是全部成像空间)进行循环，在每个空间位置的局部邻域内计算：
所述Hessian矩阵的每个单元即为点扩散函数PSF；
(3)输出PSF。
实施例
(1)使用3000m/s的常速进行PSF计算测试，成像空间为二维空间，对比测试使用的稀疏点采用等间距分布，空间分布为起始点在(200m，200m)的位置，x和z方向的增量分别为400m。
(2)分别使用反偏移-偏移方法和本发明方法计算离散点位的PSF，计算结果如图1、图2所示。
(3)对(2)所示结果选取几个典型位置进行对比，如图3、图4、图5所示，通过对比可以看出，本发明方法计算的PSF与传统方法的形态完全一致。
(4)对比本发明方法和传统反偏移-偏移计算时间，结果如图6所示，通过对比可以看出，本发明的方法在计算与传统方法相同数量的PSF时计算效率提高了11倍，而在单个PSF计算的实验中，由于传统方法同样要计算相同的反偏移-偏移次数，计算时间比计算多个PSF时并不减少，而本发明方法优于使用解析方法逐个计算PSF，当计算少量PSF时，计算时间程线性减少。
综上所述，本发明的有益效果在于：本发明通过高频近似，获取了PSF的解析表达式，提出一种不需要偏移，也不需要插值，可以计算任意空间点位PSF的方法。本方法在保证计算精度的同时，大幅提高了PSF的计算效率。
虽然已经参考优选实施例对本发明进行了描述，但在不脱离本发明的范围的情况下，可以对其进行各种改进并且可以用等效物替换其中的部件。尤其是，只要不存在结构冲突，各个实施例中所提到的各项技术特征均可以任意方式组合起来。本发明并不局限于文中公开的特定实施例，而是包括落入权利要求的范围内的所有技术方案。
Claims (10)
1.一种点扩散函数的计算方法，其特征在于，包括以下步骤：
(1)读取速度模型和观测系统，选取任意炮点和检波点并分别确定其空间坐标；
(2)计算所选取的炮点和检波点对地下空间的射线追踪；
(3)通过给定震源函数，对地下目标空间位置进行循环，根据下述公式在每个目标空间位置的局部邻域内计算：
其中，H代表Hessian矩阵，S代表震源函数，
表示炮点空间坐标，
表示检波点空间坐标，
表示离散激励的存在坐标，
是表示炮点空间坐标
出发的射线在离散激励的存在坐标
位置的方向矢量，
是表示检波点空间坐标
出发的射线在离散激励的存在坐标
位置的方向矢量，
是空间位移矢量，
是炮点空间位移矢量，
是检波点空间位移矢量；
其中，所述Hessian矩阵的每个单元即为点扩散函数。
2.根据权利要求1所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，建立Hessian矩阵的表达式：
其中，H代表Hessian矩阵，G是Green函数，uI表示入射波场，
表示炮点空间坐标，
表示检波点空间坐标，
表示空间坐标，
表示离散激励的存在坐标，t0和t′0表示时间，
是炮点空间位移矢量，
是检波点空间位移矢量。
3.根据权利要求2所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，对所建立的Hessian矩阵进行高频近似处理：
其中，S代表震源函数，A代表振幅场，τ代表走时场，δ表示狄拉克抽样函数，带入Hessian表达式，得到高频近似条件下的Hessian矩阵的表达式：
4.根据权利要求3所述的点扩散函数的计算方法，其特征在于，在步骤(1) 中，使
进一步近似的Hessian矩阵为：
5.根据权利要求4所述的点扩散函数的计算方法，其特征在于，在步骤(2)中，将进一步近似的Hessian矩阵带入射线追踪系统：
其中，
是表示炮点空间坐标
出发的射线在空间坐标
位置的方向矢量，
是空间位移矢量，得出Hessian矩阵为：
6.根据权利要求5所述的点扩散函数的计算方法，其特征在于，在步骤(3)中，忽略振幅项，得出Hessian矩阵为：
7.根据权利要求1至6中任一项所述的点扩散函数的计算方法，其特征在于，在步骤(3)中，所述地下目标空间位置包括全部成像空间。
8.根据权利要求1至6中任一项所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，采用3000m/s进行点扩散函数计算。
9.根据权利要求1至6中任一项所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，空间分布为起始点在(200m，200m)的位置。
10.根据权利要求1至6中任一项所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，x方向的增量为400m。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010536572.7A CN113805233B (zh) | 2020-06-12 | 2020-06-12 | 一种点扩散函数的计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010536572.7A CN113805233B (zh) | 2020-06-12 | 2020-06-12 | 一种点扩散函数的计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113805233A true CN113805233A (zh) | 2021-12-17 |
| CN113805233B CN113805233B (zh) | 2024-04-09 |
Family
ID=78892188
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010536572.7A Active CN113805233B (zh) | 2020-06-12 | 2020-06-12 | 一种点扩散函数的计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113805233B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112630824A (zh) * | 2019-10-09 | 2021-04-09 | 中国石油化工股份有限公司 | 一种地震成像中的离散点扩散函数生成方法及系统 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130138408A1 (en) * | 2011-11-29 | 2013-05-30 | Sunwoong Lee | Methods for Approximating Hessian Times Vector Operation in Full Wavefield Inversion |
| WO2013176579A1 (ru) * | 2012-05-23 | 2013-11-28 | Закрытое акционерное общество "Научно-инженерный центр "СИНАПС" | Измерение координат и параметров очагов при микросейсмическом мониторинге |
| KR101459388B1 (ko) * | 2014-04-18 | 2014-11-07 | 한국해양대학교 산학협력단 | 지하 속도정보 도출 방법 |
| US20150073755A1 (en) * | 2013-09-06 | 2015-03-12 | Yaxun Tang | Accelerating Full Wavefield Inversion with Nonstationary Point-Spread Functions |
| US20160180190A1 (en) * | 2014-12-22 | 2016-06-23 | The Research Foundation For The State University Of New York | Determination of spatial distribution of charged particle beams |
| US20170176613A1 (en) * | 2015-12-18 | 2017-06-22 | William A. Burnett | Method To Design Geophysical Surveys Using Full Wavefield Inversion Point- Spread Function Analysis |
| CN111158049A (zh) * | 2019-12-27 | 2020-05-15 | 同济大学 | 一种基于散射积分法的地震逆时偏移成像方法 |
-
2020
- 2020-06-12 CN CN202010536572.7A patent/CN113805233B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130138408A1 (en) * | 2011-11-29 | 2013-05-30 | Sunwoong Lee | Methods for Approximating Hessian Times Vector Operation in Full Wavefield Inversion |
| WO2013176579A1 (ru) * | 2012-05-23 | 2013-11-28 | Закрытое акционерное общество "Научно-инженерный центр "СИНАПС" | Измерение координат и параметров очагов при микросейсмическом мониторинге |
| US20150073755A1 (en) * | 2013-09-06 | 2015-03-12 | Yaxun Tang | Accelerating Full Wavefield Inversion with Nonstationary Point-Spread Functions |
| KR101459388B1 (ko) * | 2014-04-18 | 2014-11-07 | 한국해양대학교 산학협력단 | 지하 속도정보 도출 방법 |
| US20160180190A1 (en) * | 2014-12-22 | 2016-06-23 | The Research Foundation For The State University Of New York | Determination of spatial distribution of charged particle beams |
| US20170176613A1 (en) * | 2015-12-18 | 2017-06-22 | William A. Burnett | Method To Design Geophysical Surveys Using Full Wavefield Inversion Point- Spread Function Analysis |
| CN111158049A (zh) * | 2019-12-27 | 2020-05-15 | 同济大学 | 一种基于散射积分法的地震逆时偏移成像方法 |
Non-Patent Citations (2)
| Title |
|---|
| 刘伊克;常旭;卢孟夏;: "目标函数叠前保幅偏移方法与应用", 地球物理学报, no. 04 * |
| 李万万;: "基于波动方程正演的地震观测系统设计", 石油地球物理勘探, no. 02 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112630824A (zh) * | 2019-10-09 | 2021-04-09 | 中国石油化工股份有限公司 | 一种地震成像中的离散点扩散函数生成方法及系统 |
| CN112630824B (zh) * | 2019-10-09 | 2024-03-22 | 中国石油化工股份有限公司 | 一种地震成像中的离散点扩散函数生成方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113805233B (zh) | 2024-04-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107193003B (zh) | 一种稀疏奇异值分解扫描雷达前视成像方法 | |
| Zhou et al. | Tracking the direction-of-arrival of multiple moving targets by passive arrays: Algorithm | |
| CN105677942B (zh) | 一种重复轨道星载自然场景sar复图像数据快速仿真方法 | |
| Rawlinson et al. | Seismic ray tracing and wavefront tracking in laterally heterogeneous media | |
| SA519410263B1 (ar) | توليد تجميع صورة مشتركة باستخدام فصل مجال الموجة | |
| CN108072892B (zh) | 一种自动化的地质构造约束层析反演方法 | |
| CN105319589B (zh) | 一种利用局部同相轴斜率的全自动立体层析反演方法 | |
| CN112949134B (zh) | 基于非结构有限元方法的地-井瞬变电磁反演方法 | |
| CN106483559B (zh) | 一种地下速度模型的构建方法 | |
| CN105353405B (zh) | 一种全波形反演方法和系统 | |
| WO2020230214A1 (ja) | 深度推定装置、深度推定モデル学習装置、深度推定方法、深度推定モデル学習方法、及び深度推定プログラム | |
| CN109633749B (zh) | 基于散射积分法的非线性菲涅尔体地震走时层析成像方法 | |
| JP6396037B2 (ja) | データ解析装置及び方法 | |
| Han et al. | DiLO: Direct light detection and ranging odometry based on spherical range images for autonomous driving | |
| CN111665556B (zh) | 地层声波传播速度模型构建方法 | |
| CN113805233A (zh) | 一种点扩散函数的计算方法 | |
| CN109085652B (zh) | 基于改进迭代法的地空时间域电磁系统高精度下延拓方法 | |
| CN107479091B (zh) | 一种提取逆时偏移角道集的方法 | |
| CN109975869B (zh) | 一种沿地层走向光滑约束的反射波波形反演方法 | |
| AU739128B2 (en) | A method of seismic processing, and in particular a 3D seismic prospection method implementing seismic data migration | |
| Chen et al. | Potential field data interpolation by Taylor series expansion | |
| CN115267673B (zh) | 考虑重建网格偏移的稀疏声源成像方法、系统 | |
| CN109655888B (zh) | 地震数据处理中光滑浮动基准面的定量选择方法及系统 | |
| Elsherbini et al. | Image distortion effects in SAR subsurface imaging and a new iterative approach for refocusing and coregistration | |
| Lagovsky et al. | Increasing the angular resolution of control and measurement systems in signal processing |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |