CN113805233B - 一种点扩散函数的计算方法 - Google Patents

一种点扩散函数的计算方法 Download PDF

Info

Publication number
CN113805233B
CN113805233B CN202010536572.7A CN202010536572A CN113805233B CN 113805233 B CN113805233 B CN 113805233B CN 202010536572 A CN202010536572 A CN 202010536572A CN 113805233 B CN113805233 B CN 113805233B
Authority
CN
China
Prior art keywords
point
coordinates
calculating
spread function
space
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.)
Active
Application number
CN202010536572.7A
Other languages
English (en)
Other versions
CN113805233A (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.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
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 China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN202010536572.7A priority Critical patent/CN113805233B/zh
Publication of CN113805233A publication Critical patent/CN113805233A/zh
Application granted granted Critical
Publication of CN113805233B publication Critical patent/CN113805233B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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 (9)

1.一种点扩散函数的计算方法，其特征在于，包括以下步骤：
(1)读取速度模型和观测系统，选取任意炮点和检波点并分别确定其空间坐标；
(2)计算所选取的炮点和检波点对地下空间的射线追踪；
(3)通过给定震源函数，对地下目标空间位置进行循环，根据下述公式在每个目标空间位置的局部邻域内计算：
其中，H代表Hessian矩阵，S代表震源函数，表示炮点空间坐标，/>表示检波点空间坐标，/>表示离散激励的存在坐标，/>是表示炮点空间坐标/>出发的射线在离散激励的存在坐标/>位置的方向矢量，/>是表示检波点空间坐标/>出发的射线在离散激励的存在坐标/>位置的方向矢量，/>是空间位移矢量，/> 是炮点空间位移矢量，是检波点空间位移矢量，t表示时间，dt表示时间的微分；
其中，所述Hessian矩阵的每个单元即为点扩散函数；
在步骤(1)中，建立Hessian矩阵的表达式：
其中，H代表Hessian矩阵，G是Green函数，uI表示入射波场，表示炮点空间坐标，/>表示检波点空间坐标，/>表示空间坐标，/>表示离散激励的存在坐标，t0和t0′表示时间，/>是炮点空间位移矢量，/>是检波点空间位移矢量。
2.根据权利要求1所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，对所建立的Hessian矩阵进行高频近似处理：
其中，S代表震源函数，A代表振幅场，τ代表走时场，δ表示狄拉克抽样函数，带入Hessian表达式，得到高频近似条件下的Hessian矩阵的表达式：
3.根据权利要求2所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，使进一步近似的Hessian矩阵为：
4.根据权利要求3所述的点扩散函数的计算方法，其特征在于，在步骤(2)中，将进一步近似的Hessian矩阵带入射线追踪系统：
其中，是表示炮点空间坐标/>出发的射线在空间坐标/>位置的方向矢量，/>是空间位移矢量，得出Hessian矩阵为：
5.根据权利要求4所述的点扩散函数的计算方法，其特征在于，在步骤(3)中，忽略振幅项，得出Hessian矩阵为：
6.根据权利要求1至5中任一项所述的点扩散函数的计算方法，其特征在于，在步骤(3)中，所述地下目标空间位置包括全部成像空间。
7.根据权利要求1至5中任一项所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，采用3000m/s进行点扩散函数计算。
8.根据权利要求1至5中任一项所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，空间分布为起始点在(200m，200m)的位置。
9.根据权利要求1至5中任一项所述的点扩散函数的计算方法，其特征在于，在步骤(1)中，x方向的增量为400m。
CN202010536572.7A 2020-06-12 2020-06-12 一种点扩散函数的计算方法 Active CN113805233B (zh)

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 CN113805233A (zh) 2021-12-17
CN113805233B true 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)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112630824B (zh) * 2019-10-09 2024-03-22 中国石油化工股份有限公司 一种地震成像中的离散点扩散函数生成方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013176579A1 (ru) * 2012-05-23 2013-11-28 Закрытое акционерное общество "Научно-инженерный центр "СИНАПС" Измерение координат и параметров очагов при микросейсмическом мониторинге
KR101459388B1 (ko) * 2014-04-18 2014-11-07 한국해양대학교 산학협력단 지하 속도정보 도출 방법
CN111158049A (zh) * 2019-12-27 2020-05-15 同济大学 一种基于散射积分法的地震逆时偏移成像方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9176930B2 (en) * 2011-11-29 2015-11-03 Exxonmobil Upstream Research Company Methods for approximating hessian times vector operation in full wavefield inversion
US10036818B2 (en) * 2013-09-06 2018-07-31 Exxonmobil Upstream Research Company Accelerating full wavefield inversion with nonstationary point-spread functions
US9754360B2 (en) * 2014-12-22 2017-09-05 The Research Foundation For The State University Of New York Determination of spatial distribution of charged particle beams
CN108369289A (zh) * 2015-12-18 2018-08-03 埃克森美孚上游研究公司 使用全波场反演点扩展函数分析设计地球物理勘测的方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013176579A1 (ru) * 2012-05-23 2013-11-28 Закрытое акционерное общество "Научно-инженерный центр "СИНАПС" Измерение координат и параметров очагов при микросейсмическом мониторинге
KR101459388B1 (ko) * 2014-04-18 2014-11-07 한국해양대학교 산학협력단 지하 속도정보 도출 방법
CN111158049A (zh) * 2019-12-27 2020-05-15 同济大学 一种基于散射积分法的地震逆时偏移成像方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于波动方程正演的地震观测系统设计;李万万;;石油地球物理勘探(02);全文 *
目标函数叠前保幅偏移方法与应用;刘伊克;常旭;卢孟夏;;地球物理学报(04);全文 *

Also Published As

Publication number Publication date
CN113805233A (zh) 2021-12-17

Similar Documents

Publication Publication Date Title
Zhou et al. Tracking the direction-of-arrival of multiple moving targets by passive arrays: Algorithm
US8436912B2 (en) Range measurement using multiple coded apertures
CN103503022B (zh) 提供构建目标物体的区域的图像的图像数据的方法和装置
US8330852B2 (en) Range measurement using symmetric coded apertures
US20110267485A1 (en) Range measurement using a coded aperture
CN112949134B (zh) 基于非结构有限元方法的地-井瞬变电磁反演方法
CA2961572C (en) Velocity tomography using property scans
WO2008069364A1 (en) 4-d inversion of geophysical data and 4-d imaging method of geologic structure using it
CN109633749B (zh) 基于散射积分法的非线性菲涅尔体地震走时层析成像方法
JP6396037B2 (ja) データ解析装置及び方法
Lagovsky et al. Image Restoration of Two-dimensional Signal Sources with Superresolution.
CN113805233B (zh) 一种点扩散函数的计算方法
CN110765631B (zh) 基于有效成像像素的红外辐射特性测量小目标判断方法
Ravetta et al. Noise source localization and optimization of phased-array results
KR101814023B1 (ko) 유한차분격자 자료 자동 보정 장치 및 방법
Gurbuz et al. Sparse ground-penetrating radar imaging method for off-the-grid target problem
CN116449305A (zh) 基于可控变分自编码器的稠密时变阵列构建方法及系统
Rajani et al. Direction of arrival estimation by using artificial neural networks
US20100092106A1 (en) Quasi-quadratic interpolation of image data
CN110118968A (zh) 倾斜四反射板镜像综合孔径辐射计及成像方法
CN109061592B (zh) 基于压缩感知的多点发射毫米波雷达测向方法
CN112630840A (zh) 基于统计特征参数的随机反演方法及处理器
Ozdemir et al. ACSAR-antenna coupling synthetic aperture radar imaging algorithm
US20240125935A1 (en) Arithmetic operation system, training method, and non-transitory computer readable medium storing training program
US20240126953A1 (en) Arithmetic operation system, training method, and non-transitory computer readable medium storing training program

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