CN115061197A - 二维海面鬼波水体成像测量方法、系统、终端及测流设备 - Google Patents

二维海面鬼波水体成像测量方法、系统、终端及测流设备 Download PDF

Info

Publication number
CN115061197A
CN115061197A CN202210592015.6A CN202210592015A CN115061197A CN 115061197 A CN115061197 A CN 115061197A CN 202210592015 A CN202210592015 A CN 202210592015A CN 115061197 A CN115061197 A CN 115061197A
Authority
CN
China
Prior art keywords
ghost
reflection
frequency domain
wave field
going
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.)
Pending
Application number
CN202210592015.6A
Other languages
English (en)
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.)
Ocean University of China
Original Assignee
Ocean University of China
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 Ocean University of China filed Critical Ocean University of China
Priority to CN202210592015.6A priority Critical patent/CN115061197A/zh
Publication of CN115061197A publication Critical patent/CN115061197A/zh
Pending legal-status Critical Current

Links

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/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • 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/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/38—Seismology; Seismic or acoustic prospecting or detecting specially adapted for water-covered areas
    • G01V1/3808—Seismic data acquisition, e.g. survey design
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/10—Aspects of acoustic signal generation or detection
    • G01V2210/14—Signal detection
    • G01V2210/142—Receiver location
    • G01V2210/1423—Sea
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/63—Seismic attributes, e.g. amplitude, polarity, instant phase

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Oceanography (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明属于图像处理技术领域，公开了二维海面鬼波水体成像测量方法、系统、终端及测流设备。获取地震数据的频率域正变换数据；对于任意频率分量及射线参数，计算得到虚反射与一次反射的延迟时间，求解得到第一次的下行波场频率域Radon变换结果；获取射线参数所对应的虚反射与一次反射的平均延迟时间和平均振幅差异系数并再次计算，不断交替迭代直至满足收敛条件，得到最终的虚反射对应的频率域Radon变换结果，再进行频率域Radon反变换、时间域傅里叶反变换，得到最终的虚反射预测数据。本发明能够对ADCP接收到的虚反射信号等有用信息进行全波形反演来获得表层海水的性质信息，对水体实现精确成像。

Description

二维海面鬼波水体成像测量方法、系统、终端及测流设备
技术领域
本发明属于图像处理技术领域，尤其涉及二维海面鬼波水体成像测量方法、系统、终端及测流设备。
背景技术
ADCP经过30年的发展已形成多种形式和不同用途的系列化产品，从换能器阵结构上区分有采用四个圆形换能器组成的常规阵和由近千个小换能器组成的相控阵。从用途上区分有自容式ADCP、直读/走航式ADCP(走航式ADCP全部采用相控阵技术)和下放式ADCP(简称LowADCP)。
以信号检测技术为特征，ADCP的研究和发展可分为三个阶段。第一阶段是以锁相式频率跟踪器为特征，其测流的空间分辨率和精度均较低。第二阶段是以复相关技术作为计算信号多普勒频移的手段，如RD公司的窄带声学多普勒海流剖面仪(NBADCP)。它的优点是可同步测量多层海流、体积小、设备简单，稳定性好，可长期工作。缺点是受频率和脉宽的限制，剖面厚度不能太小，单次测量的方差大。第三阶段是采用宽带信号测流技术的宽带声学多普勒海流剖面仪(BBADCP)，它以宽带脉冲对相干技术为特征，发射一个伪随机编码脉冲，该脉冲是由N个码元对组成。剖面测量深度取决于整个脉冲的长度和频率，测量层厚度取决于码元的长度。由于码元的宽度远小于脉冲的宽度，每一个码元对应给出一个海流测量值，因而BBADCP的测层厚度大大减小，单次测量的方差减小到NBADCP单次测量方差的1N。因此，BBADCP受到海洋界的广泛欢迎。
ADCP就是利用声波的多普勒效应发展起来的一种新型测流设备，它既能测量相对水底速度，又可兼顾测得相对水流速度。然而，现有技术中并没有基于ADCP技术进行震源与海面之间水体的探测成像方法。
通过上述分析，现有技术存在的问题及缺陷为：对于虚反射信号的处理反演分辨率以及准确度存在进一步优化的空间。现有技术对ADCP接收到的虚反射信号等有用信息不能进行反演，使得获得表层海水的水体不能精确成像。
解决以上问题及缺陷的难度为：将ADCP信息系统与常规虚反射建立联系，并应用虚反射信息准确反演成像。
解决以上问题及缺陷的意义为：将常规处理是屏蔽滤波的虚反射信息加以提取利用，提高了信号的信噪比，丰富的信号信息有助于高精度分析。
发明内容
为克服相关技术中存在的问题，本发明公开实施例提供了一种二维海面鬼波水体成像测量方法、系统、终端及测流设备。具体涉及一种海洋勘探中的水流剖面成像方法。
所述技术方案如下：一种二维海面鬼波水体成像测量方法，包括：
对测流设备ADCP接收到的富含水体信息的虚反射信号进行全波形反演，获得表层海水的性质信息，对水体进行精确成像。
在一实施例中，所述二维海面鬼波水体成像测量方法具体包括：
(1)对地震记录S(xi，t)做时间方向的傅里叶变换，得到向量S(x，ω)，利用高精度Radon变换，得到地震数据的频率域Radon正变换数据Us(pk，ω)；
(2)对于任意频率分量ω及射线参数pk，给定虚反射与一次反射的振幅差异系数Rk＝-1，并计算得到虚反射与一次反射的延迟时间Δτk，构建出矩阵G(xi，pk，ω)，求解得到第一次的下行波场频率域Radon变换结果
Figure BDA0003665716580000031
(3)通过交替迭代反演目标函数，得到射线参数pk所对应的虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000032
和平均振幅差异系数
Figure BDA0003665716580000033
(4)再利用步骤(3)得到
Figure BDA0003665716580000034
和
Figure BDA0003665716580000035
利用步骤(2)得到下行波场
Figure BDA0003665716580000036
(5)再次计算虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000037
和平均振幅差异系数
Figure BDA0003665716580000038
不断交替迭代直至满足收敛条件，得到最终的虚反射对应的频率域Radon变换结果UDown(p，ω)；
(6)对UDown(p，ω)进行频率域Radon反变换，并进行时间域傅里叶反变换，则得到最终的虚反射预测数据uDown(x，t)。
在一实施例中，所述步骤(3)得到射线参数pk所对应的虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000039
和平均振幅差异系数
Figure BDA00036657165800000310
包括：
在采集中，第i个检波器的深度为zi，偏移距为xi；p为与检波器相匹配的射线参数；
其中，第k个射线参数
Figure BDA0003665716580000041
θk为地震波传播至海平面的出射角，虚反射的入射角；Vw为海水速度；该检波器处接收到的地震波场为上行波场和下行波场之和，则
uS(pk，τ)＝uUp(pk，τ)+uDown(pk，τ)
Figure BDA0003665716580000042
其中，uS为总波场，uUp为上行波场，uDown为下行波场；τ为地震波在偏移距为0处的时间，Δτ为上行波场与下行波场的时差，Rk为鬼波与一次反射波振幅的比值系数；
根据地震波传播的几何关系，下行波场与上行波场之间的时间差Δτ为：
Figure BDA0003665716580000043
又因为接收的总波场由Nm个水平慢度p表示：
Figure BDA0003665716580000044
则总波场与虚反射对应的下行波场在频率域的表达式为：
Figure BDA0003665716580000045
S(xi，ω)为地震数据的时间方向的傅里叶变换结果，若：
Figure BDA0003665716580000046
则虚反射波场的最小二乘解为：
UDown(p，ω)＝(GH(x，p，ω)G(x，p，ω))-1G(x，p，ω)HUs(x，ω)
忽略炮检距误差，未知参数包括下行波场与一次波场的时差Δτ、虚反射与一次反射能量比值系数R，以及待求的下行波场UDown(p，ω)；为求得上述三个参数，利用交替迭代的反演方法，固定下行波场与一次波场的时差Δτ、虚反射与一次反射能量比值系数R，以及待求的下行波场UDown(p，ω)中的两个参数，反演一个参数，然后再将反演得到参数作为约束值，反演另外一个参数，由此不断反复得到最优解。
在一实施例中，在开始进行交替迭代反演过程前，需确定一个反演的目标函数；上行波场与下行波场之间的关系用两种表达形式：
uS(pk，τ)＝uUp(pk，τ)+uDown(pk，τ)
Figure BDA0003665716580000051
进行傅里叶变换得：
UUp(pk，ω)＝Us(pk，ω)-UDown(pk，ω)
Figure BDA0003665716580000052
其中，Us(pk，ω)为总波场的频率域正变换；则交替反演的目标利用反演得到下行波场，通过上述
Figure BDA0003665716580000053
公式，所得到的一次波场相等，即
Figure BDA0003665716580000054
在一实施例中，所述二维海面鬼波水体成像测量方法进一步包括：
虚反射包含地下构造的地震响应，将检波器假象移动至真实检波器相对海平面的海面以上的镜像位置，将海面上的速度等价看作海水速度；将虚反射作为新的一次反射，利用常规偏移成像方法进行偏移成像；同时还将该虚反射成像结果与一次反射成像结果进行匹配相关，得到虚反射与一次反射的联合成像结果。
本发明的另一目的在于提供一种实施所述二维海面鬼波水体成像测量方法的二维海面鬼波水体成像测量系统，所述二维海面鬼波水体成像测量系统包括：
地震数据频率域正变换数据获取模块，用于对地震记录S(xi，t)做时间方向的傅里叶变换，得到向量S(x，ω)，利用高精度Radon变换，得到地震数据的频率域Radon正变换数据Us(pk，ω)；
第一次下行波场频率域变换结果获取模块，用于对于任意频率分量ω及射线参数pk，给定虚反射与一次反射的振幅差异系数Rk＝-1，并计算得到虚反射与一次反射的延迟时间Δτk，构建出矩阵G(xi，pk，ω)，求解得到第一次的下行波场频率域Radon变换结果
Figure BDA0003665716580000061
预测参数获取模块，用于通过交替迭代反演目标函数，得到射线参数pk所对应的虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000062
和平均振幅差异系数
Figure BDA0003665716580000063
下行波场获取模块，用于再利用步骤(3)得到
Figure BDA0003665716580000064
和
Figure BDA0003665716580000065
利用步骤(2)得到下行波场
Figure BDA0003665716580000066
虚反射对应的频率域变换结果获取模块，用于再次计算虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000067
和平均振幅差异系数
Figure BDA0003665716580000068
不断交替迭代直至满足收敛条件，得到最终的虚反射对应的频率域Radon变换结果UDown(p，ω)；
最终虚反射预测数据获取模块，用于对UDown(p，ω)进行频率域Radon反变换，并进行时间域傅里叶反变换，则得到最终的虚反射预测数据uDown(x，t)。
本发明的另一目的在于提供一种信息数据处理终端，所述信息数据处理终端包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行所述二维海面鬼波水体成像测量方法。
本发明的另一目的在于提供一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行所述二维海面鬼波水体成像测量方法。
本发明的另一目的在于提供一种测流设备ADCP，所述测流设备ADCP执行所述二维海面鬼波水体成像测量方法。
本发明的另一目的在于提供一种存储在计算机可读介质上的计算机程序产品，包括计算机可读程序，供于电子装置上执行时，提供用户输入接口以实施所述二维海面鬼波水体成像测量方法。
结合上述的所有技术方案，本发明所具备的优点及积极效果为：创新性的利用了虚反射信息，丰富勘探有效信号成分，既可以佐证前期信号处理结果，又提高信噪比，提供了丰富的信号信息。
虚反射信号中富含水体信息，本发明能够对ADCP接收到的虚反射信号等有用信息进行全波形反演来获得表层海水的性质信息，对水体实现精确成像。
当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本发明的公开。
附图说明
此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本公开的实施例，并与说明书一起用于解释本公开的原理。
图1是本发明实施例提供的二维海面鬼波水体成像测量方法流程图。
图2是本发明实施例提供的ADCP发射声波的JANUS结构图。
图3是本发明实施例提供的具有JANUS结构的ADCP工作示意图.
具体实施方式
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似改进，因此本发明不受下面公开的具体实施的限制。
如图1所示，本发明提供一种二维海面鬼波水体成像测量方法包括：
S101,对地震记录S(xi，t)做时间方向的傅里叶变换，得到向量S(x，ω)，利用高精度Radon变换，得到地震数据的频率域Radon正变换数据Us(pk，ω)；
S102,对于任意频率分量ω及射线参数pk，给定虚反射与一次反射的振幅差异系数Rk＝-1，并计算得到虚反射与一次反射的延迟时间Δτk，构建出矩阵G(xi，pk，ω)，求解得到第一次的下行波场频率域Radon变换结果
Figure BDA0003665716580000081
S103,通过交替迭代反演目标函数，得到射线参数pk所对应的虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000091
和平均振幅差异系数
Figure BDA0003665716580000092
S104,再次利用步骤S103得到新的
Figure BDA0003665716580000093
和
Figure BDA0003665716580000094
然后利用步骤S102得到第二次下行波场频率域Radon变换结果
Figure BDA0003665716580000095
S105,再次计算虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000096
和平均振幅差异系数
Figure BDA0003665716580000097
不断交替迭代直至满足收敛条件，得到最终的虚反射对应的频率域Radon变换结果UDown(p，ω)；
S106,对虚反射对应的频率域Radon变换结果UDown(p，ω)进行频率域Radon反变换，并进行时间域傅里叶反变换，则得到最终的虚反射预测数据uDown(x，t)。
本发明的另一目的在于提供一种实施所述二维海面鬼波水体成像测量方法的二维海面鬼波水体成像测量系统，所述二维海面鬼波水体成像测量系统包括：
地震数据频率域正变换数据获取模块，用于对地震记录S(xi，t)做时间方向的傅里叶变换，得到向量S(x，ω)，利用高精度Radon变换，得到地震数据的频率域Radon正变换数据Us(pk，ω)；
第一次下行波场频率域变换结果获取模块，用于对于任意频率分量ω及射线参数pk，给定虚反射与一次反射的振幅差异系数Rk＝-1，并计算得到虚反射与一次反射的延迟时间Δτk，构建出矩阵G(xi，pk，ω)，求解得到第一次的下行波场频率域Radon变换结果
Figure BDA0003665716580000098
预测参数获取模块，用于通过交替迭代反演目标函数，得到射线参数pk所对应的虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000099
和平均振幅差异系数
Figure BDA00036657165800000910
下行波场获取模块，用于再利用步骤(3)得到
Figure BDA0003665716580000101
和
Figure BDA0003665716580000102
利用步骤(2)得到下行波场
Figure BDA0003665716580000103
虚反射对应的频率域变换结果获取模块，用于再次计算虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000104
和平均振幅差异系数
Figure BDA0003665716580000105
不断交替迭代直至满足收敛条件，得到最终的虚反射对应的频率域Radon变换结果UDown(p，ω)；
最终虚反射预测数据获取模块，用于对UDown(p，ω)进行频率域Radon反变换，并进行时间域傅里叶反变换，则得到最终的虚反射预测数据uDown(x，t)。
下面结合具体实施例对本发明的技术方案作进一步描述。
由于本发明提出的基于ADCP装置接收到的有用信息全波形反演来进行水体精确成像等相关内容，下面对相关内容进行介绍，以便本领域技术人员对本发明的内容更加清楚。
一、鬼波的产生原理
为了增加激发能量、降低海水涌浪噪声，海上地震勘探往往将震源与检波器布置在海水中。但海水阻抗近似为1，空气阻抗为0，前者阻抗大于后者，两者之间会出现一个反射系数近似为-1的反射界面。当震源在海平面下激发后，会沿上、下两个方向传播，当向上传播至海平面时，会出现极性翻转的地震反射，即虚反射，也称为鬼波。鬼波的传播路径大致有3类，按照产生机制通常可定义为：震源虚反射、检波器虚反射以及震源和检波器虚反射。其中，震源虚反射是指震源激发后向上传播的能量经海水-空气阻抗界面反射，再经海底反射至检波器；检波器虚反射是指震源部分能量向下传播，经地下界面反射后，再经海水-空气波阻抗界面反射至检波器；震源和检波器虚反射是指包含震源能量经海水-空气阻抗界面反射，向下传播至海底，又反射至海水-空气波阻抗界面，最后由检波器接收。
二、ADCP测流原理
ADCP—声学多普勒海流剖面仪是目前观测上层海流剖面的最有效方法。其基本原理是测定声波入射到海水中微颗粒后向散射在频率上的多普勒频移，从而得到不同水层水体的运动速度。其特点是精度高、分辨率高，操作方便。
ADCP根据测定水体中微颗粒声后向散射的多普勒频移来测量水体速度。在仪器坐标下，每一声束方向上每一深度段上的水体的流速分量可以根据在该段水体上测得的多普勒频率(D)，由下面方程计算出:
Figure BDA0003665716580000111
式中Ft为发射频率，单位为kHz。C为水中的声速，可按下列公式计算:
C＝1449.2+4.6T-5.5*10-3T3+2.9*10-4T3+(1.34-10-2T)
式中:
D＝深度，单位为m；
S＝盐度，单位为10-3；
T＝温度，单位为℃。
利用四声束(至少三束)测得的水体回反散射的多普勒频移，便可以求得三维流速，Vx，Vy和Vz。并且可以转换为地球坐标下的U(东分量)，V(北分量)和W(垂直分量)。由于声速在一定水域中，在一定深度范围内(表层至几百米深度)的水体中的传播速度基本是不变的，根据由声波发射到接收的时间差，便可以确定深度。利用不断发射的声脉冲，确定一定的发射时间间隔及滞后，通过对多普勒频移及谱宽度的估计运算，便可以得到整个水体剖面逐层段上水体的流速。
一般情况下，海流速度在一定范围内可以认为是不变化的，所以，可将海流速度定义为东西方向、南北方向和垂线方向的三维速度矢量。为测量这三个方向上的速度矢量，ADCP通常同时向海水中的东、西、南、北四个方向上同时发射与垂线成夹角为φ的四个脉冲信号，如图2所示，该波束结构称为JANUS结构。其工作示意图如图3所示。
对于海上地震处理来说，通常采用消除虚反射，之后利用常规偏移成像方法对一次反射波进行偏移成像。而虚反射同样包含水层的地震响应，对其单独进行成像可采用镜像电缆方法，即将检波器假象移动至真实检波器相对海平面的海面以上的镜像位置，此时将海面上的速度等价看作海水速度，此时，便可将虚反射看作一类新的一次反射，利用常规偏移成像方法对其进行偏移成像。
三、ADCP波形求取
1、算法原理
传统的FWI算法是基于求取使用模拟数据与实测数据之间的残差构建的目标函数公式(1)的极小值实现的，模拟合成记录通常是使用当前速度模型由波动方程正演得到。
Figure BDA0003665716580000131
该极值求解过程是从目标函数构建一个梯度函数公式(2)来完成，目前常用的方法有共轭梯度法、最速下降法、高斯牛顿法等。与全局算法相反，各种局部极值算法的最大区别在于计算机时成本。
Figure BDA0003665716580000132
2、实现过程
2.1子波求取
从ADCP波形求取的基本原理可以看出，模拟合成记录与实测记录之间的残差驱动了反演算法，因此求取一个用于正演的准确子波是ADCP波形求取技术的基础也是非常关键的一步。常见的子波求取算法有：用气枪阵列参数模拟子波，从ADCP资料直达波中提取子波，使用近场记录(NFH)合成的远场子波等方法。本方法此处通过ADCP装置接收到的信息中提取子波。
2.2ADCP数据准备
ADCP数据准备是指对用于ADCP波形求取的数据做的一些基础性处理，也可称为ADCP数据预处理。ADCP波形求取技术一个潜在的问题是周期跳跃，也就是寻优过程陷入局部极值的问题。
有效的低频信号对于反演的成败非常关键。因此获得低频成分丰富的ADCP资料是基础。对于所获资料来说，首先应用一个低截滤波，消除低频直流分量以及强能量涌浪噪声等。多次波也是一般海洋资料中主要的干扰信号，是否对数据多次波压制也是讨论比较多的一个方面，需要配合其他参数综合考虑，最主要的是选用何种震源子波。因此在数据的准备过程中，只需要做一些必须的基础工作即可，相对于常规生产处理来说，工作量明显减少，这也是为什么业界有人提出反演可以缩短生产周期的说法。
2.3初始模型建立
ADCP波形求取严格要求一个平滑但是相对准确的速度模型，好的初始模型对于ADCP波形求取更新至关重要，如果使用初始模型正演得到的模拟数据与实测数据误差过大，反演过程很容易陷入局部极值，因此需要常规的速度建模手段对速度模型做适当优化后作为反演的初始模型，并且在输入ADCP波形求取使用之前做仔细的检查分析。初始模型可通过井资料，时间域处理以及旅行时反演等方法获取。
本发明方法使用的是时间域ADCP波形求取算法，采用分频段逐步提高反演频率的策略。利用尽可能低的有效信号能量开始反演更新，确保反演的结果稳定收敛，然后逐步提高反演频率，因此反演的起始频率是一个重要参数。
其他关键参数包括最大反演频率等，ADCP波形求取是一个迭代寻优的过程，每个迭代都包含有两个波动方程算法。
2.4虚反射提取方法
虚反射与一次反射之间的时差与波形特征在时间域剖面中仅存在着微弱的差异，直接计算比较难以区分这种差异性，从而预测得到虚反射地震记录。为此，滤波法多次波消除多将时间域剖面转换至FK域或Radon域进行处理，从而得到性质更优良的特征，利用一次波与虚反射之间的时差关系将其分离。同时，一次波虚反射间时差和虚反射反射系数的估计也会影响虚反射预测的精度，以适应海底构造起伏、速度变化给虚反射预测、消除带来的复杂性。
2.4.1基于Majorization-Minimization框架的高精度Radon变换
Radon变换是一种地震处理中常用的多次波消除方法，任何方程、图像都可以用其自身的投影来重建。该理论最先被应用于CT扫描，通过设置X射线的路径来计算某一个角度下物体的积分投影值，对这些值进行反变换从而恢复原图像。假设射线的位置为r＝xsin θ+ycos θ，则在该角度下射线投影的总衰减为：
Figure BDA0003665716580000151
其中，f(x，y)待重建的二维数据。
2.4.2Radon变换及其存在的关键问题分析
时间域的Radon变换是按照特定的射线轨迹对地震剖面进行线性积分，其表达式为：
Figure BDA0003665716580000152
其中，m(τ，p)为时间域数据d(t，x)在Radon正变换后的转换数据，t＝φ(τ，p，x)为积分路径，t，τ，p，x表示时间、偏移距、双程时间截距以及曲率参数。由Radon域数据m(τ，p)恢复d(t，x)的变换公式为：
Figure BDA0003665716580000153
对上述公式进行离散化，沿着偏移距xi(i＝1，…，Nx)进行地震振幅求和，则Radon变换的时间域具体公式为：
Figure BDA0003665716580000154
设射线参数p的数目为NM个，对应的反变换形式为：
Figure BDA0003665716580000161
将Radon正变换进行矩阵表示，可得
m＝LTd (3-6)
通常为了求取m，首先定义反变换算子L，
d＝Lm (3-7)
构建反演目标函数
Figure BDA0003665716580000162
为保证m可以恢复原始数据，利用最优化原理可得m的最小二乘解
m＝(LLT)-1LTd (3-9)
可看出，时间域变换中二维地震数据d的维数为Nx×Nt，m的维数与离散的射线参数有关，则m的大小为NM×Nt，所以LT的大小应该(NM×Nt)×(Nx×Nt)，根据公式矩阵表达形式，满足矩阵相乘定律的LLT的大小是(NM×Nt)×(Nx×Nt)。
为了恢复原始数据，要保证射线参数足够丰富，导致M需要足够大，使得LLT规模过大，这将导致LLT需要极大的存储空间，而LTd也需要(NM×Nt)次(Nx×Nt)维向量点积和向量求和运算，计算开销大。为此，Radon域变换往往需要在频率域中计算，但前提是所用的积分路径形式应满足时不变特征。
将t-x数据转换为频率域，可按每道进行时间域的傅里叶正变换。
则，Radon正变换的频率域表达公式为：
Figure BDA0003665716580000171
对应的反变换形式为：
Figure BDA0003665716580000172
其中，ω＝2πf为角频率，f为频率。
其矩阵表达为：
Figure BDA0003665716580000173
对于任意一个频率来说，LH(ω)大小为Nx×Nt。
此时，目标函数变为
Figure BDA0003665716580000174
则根据最小二乘求解方法：
M＝(LLH)-1LHd (3-14)
则Radon正变换可在频率域里分频点进行，由于该求解过程属于多解问题，为了保证求解的稳定性，通常施加阻尼项u。
M＝(LLH+μI)-1LHd (3-15)
根据上述分析可得Radon频率域正变换的计算流程：
Step1：对t-x地震数据进行时间域的傅里叶变换；
Step2：对每个频率分别利用公式求解M；
Step3：对求解得到M施加时间域傅里叶逆变换获取Radon正变换τ-p数据。
从上述计算过程中不难看出，射线路径形式、射线参数p的取值范围、阻尼项大小以及反演求解算法的选择都会影响Radon变换应用的精度和分辨率，而这对后续处理来说是至关重要的。
Radon变换在地震勘探中有着广泛的应用，出现了大量的改进算法。在射线路径选择方面，主要采用了线性、抛物线、双曲线等，分别针对平面波分解、速度叠加及正常时差校正之后的多次波消除等实际应用。
2.5基于交替反演的Radon域虚反射预测关键参数计算方法
2.5.1基本原理
虽然Radon域处理算法可以较好地预测虚反射，但虚反射的预测过程需要将原始地震数据转换为Radon域数据，而实际数据形态特征复杂Radon域数据的计算过程中射线参数、反演约束参数(阻尼因子、稀疏性约束权重)选择不当都会导致变换结果出现偏差，从而使得虚反射预测结果偏离预期效果。
为不失一般性，假设为变深度缆采集，其中第i个检波器的埋深为zi，偏移距为xi。p为与之相匹配的射线参数。
其中，
Figure BDA0003665716580000181
θk为地震波传播至海平面的出射角，虚反射的入射角；Vw为海水速度。同时可知该检波器处接收到的地震波场为上行波场和下行波场之和，则
Figure BDA0003665716580000182
其中，uS为总波场，uUp为上行波场，uDown为下行波场；τ为地震波在偏移距为0处的时间，Δτ为上行波场与下行波场的时差，Rk为鬼波与一次反射波振幅的比值系数；
根据地震波传播的几何关系，可知下行波场与上行波场之间的时间差Δτ为：
Figure BDA0003665716580000191
又因为接收的总波场可由Nm个水平慢度p表示：
Figure BDA0003665716580000192
则总波场与虚反射对应的下行波场在频率域的表达式为：
Figure BDA0003665716580000193
S(xi，ω)为地震数据的时间方向的傅里叶变换结果，假设：
Figure BDA0003665716580000194
则虚反射波场的最小二乘解为：
UDown(p，ω)＝(GH(x，p，ω)G(x，p，ω))-1G(x，p，ω)HUs(x，ω) (3-21)
忽略炮检距误差，未知参数包括下行波场与一次波场的时差Δτ、虚反射与一次反射能量比值系数R，以及待求的下行波场UDown(p，ω)。为求得上述三个参数，需要利用交替迭代的反演方法，即固定其中两个参数，反演一个参数，然后再将反演得到参数作为约束值，反演另外一个参数，由此不断反复得到最优解。
在开始进行交替迭代反演过程之前，首先需要确定一个反演的目标函数。由上述论述可知，上行波场与下行波场之间的关系可用两种表达形式：
Figure BDA0003665716580000195
对其进行傅里叶变换可得：
Figure BDA0003665716580000201
其中，Us(pk，ω)为总波场的频率域正变换。则交替反演的目标应保证利用反演得到下行波场，通过上述公式，所得到一次波场相等，即
Figure BDA0003665716580000202
根据之前分析的基于Majorization-Minimization优化框架的Radon域虚反射预测方法，当Δτ固定时，可采用L1范数获得高分辨的Radon域变换结果。因此，基于交替反演的Radon域虚反射预测流程为：
(1)对地震记录S(xi，t)做时间方向的傅里叶变换，得到向量S(x，ω)，利用高精度Radon变换，得到地震数据的频率域Radon正变换数据Us(pk，ω)；
(2)对于任意频率分量ω及射线参数pk，给定虚反射与一次反射的振幅差异系数Rk＝-1，并计算得到虚反射与一次反射的延迟时间Δτk，从而构建矩阵G(xi，pk，ω)，求解得到第一次的下行波场频率域Radon变换结果
Figure BDA0003665716580000203
(3)通过交替迭代反演目标函数，得到射线参数pk所对应的虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000204
和平均振幅差异系数
Figure BDA0003665716580000205
(4)再利用步骤(3)得到
Figure BDA0003665716580000206
和
Figure BDA0003665716580000207
利用步骤(2)得到下行波场
Figure BDA0003665716580000208
(5)再次计算虚反射与一次反射的平均延迟时间
Figure BDA0003665716580000209
和平均振幅差异系数
Figure BDA0003665716580000211
不断交替迭代直至满足收敛条件，得到最终的虚反射对应的频率域Radon变换结果UDown(p，ω)；
(6)对UDown(p，ω)进行频率域Radon反变换，并进行时间域傅里叶反变换，则得到最终的虚反射预测数据uDown(x，t)。
由上述过程反演得到虚反射数据能量与地震数据的能量是相匹配的，即地震数据减去预测得到虚反射数据即为一次反射数据。
2.6虚反射成像基本思路
对于海上地震处理来说，通常采用消除虚反射，之后利用常规偏移成像方法对一次反射波进行偏移成像。而虚反射同样包含地下构造的地震响应，将检波器假象移动至真实检波器相对海平面的海面以上的镜像位置，此时将海面上的速度等价看作海水速度，此时，便可将虚反射看作一类新的一次反射，利用常规偏移成像方法对其进行偏移成像。同时还将该虚反射成像结果与一次反射成像结果进行匹配相关，得到虚反射与一次反射的联合成像结果，进一步提高地下构造成像、提高精度。
本领域技术人员在考虑说明书及实践这里公开的公开后，将容易想到本公开的其它实施方案。本申请旨在涵盖本公开的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本公开的一般性原理并包括本公开未公开的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本公开的真正范围和精神由所附的权利要求指出。
应当理解的是，本公开并不局限于上面已经描述并在附图中示出的精确结构，并且可以在不脱离其范围进行各种修改和改变。本公开的范围应由所附的权利要求来限制。

Claims (10)

1.一种二维海面鬼波水体成像测量方法，其特征在于，该二维海面鬼波水体成像测量方法通过测流设备ADCP接收富含水体信息，提取所过滤的虚反射的有效信息，并进行全波形反演，获得表层海水的性质信息，对水体信号处理成像，具体包括以下步骤：
(1)对地震记录S(xi，t)做时间方向的傅里叶变换，得到向量S(x，ω)，利用Radon变换，得到地震数据的频率域Radon正变换数据Us(pk，ω)；
(2)对于任意频率分量ω及射线参数pk，给定虚反射与一次反射的振幅差异系数Rk＝-1，并计算得到虚反射与一次反射的延迟时间Δτk，构建出矩阵G(xi，pk，ω)，求解得到第一次的下行波场频率域Radon变换结果
Figure FDA0003665716570000011
(3)通过交替迭代反演目标函数，得到射线参数pk所对应的虚反射与一次反射的平均延迟时间
Figure FDA0003665716570000012
和平均振幅差异系数
Figure FDA0003665716570000013
(4)通过步骤(3)得到
Figure FDA0003665716570000014
和
Figure FDA0003665716570000015
利用步骤(2)的计算得到下行波场
Figure FDA0003665716570000016
(5)再次计算虚反射与一次反射的平均延迟时间
Figure FDA0003665716570000017
和平均振幅差异系数
Figure FDA0003665716570000018
不断交替迭代直至满足收敛条件，得到最终的虚反射对应的频率域Radon变换结果
Figure FDA0003665716570000019
(6)对UDown(p，ω)进行频率域Radon反变换，并进行时间域傅里叶反变换，则得到最终的虚反射预测数据uDown(x，t)。
2.根据权利要求1所述的二维海面鬼波水体成像测量方法，其特征在于，在步骤(3)中，得到射线参数pk所对应的虚反射与一次反射的平均延迟时间
Figure FDA0003665716570000021
和平均振幅差异系数
Figure FDA0003665716570000022
包括：
通过具有独立信号接收的检波器接收信号并输入到测流设备ADCP进行处理，在采集中，第i个检波器的深度为zi，偏移距为xi；p为与检波器相匹配的射线参数；
其中，第k个射线参数
Figure FDA0003665716570000023
θk为地震波传播至海平面的出射角，虚反射的入射角；Vw为海水速度；该检波器处接收到的地震波场为上行波场和下行波场之和，则
Figure FDA0003665716570000024
Figure FDA0003665716570000025
其中，uS为总波场，uUp为上行波场，uDown为下行波场；τ为地震波在偏移距为0处的时间，Δτ为上行波场与下行波场的时差，Rk为鬼波与一次反射波振幅的比值系数；
根据地震波传播的几何关系，下行波场与上行波场之间的时间差Δτ为：
Figure FDA0003665716570000026
接收的总波场由Nm个水平慢度p表示：
Figure FDA0003665716570000027
则总波场与虚反射对应的下行波场在频率域的表达式为：
Figure FDA0003665716570000028
S(xi，ω)为地震数据的时间方向的傅里叶变换结果，若：
Figure FDA0003665716570000029
则虚反射波场的最小二乘解为：
UDown(p，ω)＝(GH(x，p，ω)G(x，p，ω))-1G(x，p，ω)HUS(x，ω)。
3.根据权利要求2所述的二维海面鬼波水体成像测量方法，其特征在于，求得下行波场与一次波场的时差Δτ、虚反射与一次反射能量比值系数R，待求的下行波场UDown(p，ω)中，利用交替迭代的反演方法，固定下行波场与一次波场的时差Δτ、虚反射与一次反射能量比值系数R，以及待求的下行波场UDown(p，ω)中的两个参数，反演一个参数，然后再将反演得到参数作为约束值，反演另外一个参数，由此不断反复得到最优解。
4.根据权利要求2所述二维海面鬼波水体成像测量方法，其特征在于，在开始进行交替迭代反演过程前，需确定一个反演的目标函数；上行波场与下行波场之间的关系表达形式包括：
uS(pk，τ)＝uUp(pk，τ)+uDowm(pk，τ)
Figure FDA0003665716570000031
进行傅里叶变换得：
Figure FDA0003665716570000032
Figure FDA0003665716570000033
其中，Us(pk，ω)为总波场的频率域正变换；则交替反演的目标利用反演得到下行波场，通过上述
Figure FDA0003665716570000034
公式，得到的一次波场相等，为：
Figure FDA0003665716570000035
5.根据权利要求1所述二维海面鬼波水体成像测量方法，其特征在于，所述二维海面鬼波水体成像测量方法还包括：
虚反射包含地下构造的地震响应，将检波器假象移动至真实检波器相对海平面的海面以上的镜像位置，将海面上的速度等价看作海水速度；将虚反射作为新的一次反射，利用常规偏移成像方法进行偏移成像；同时还将该虚反射成像结果与一次反射成像结果进行匹配相关，得到虚反射与一次反射的联合成像结果。
6.一种实施权利要求1-5任意一项所述二维海面鬼波水体成像测量方法的二维海面鬼波水体成像测量系统，其特征在于，所述二维海面鬼波水体成像测量系统包括：
地震数据频率域正变换数据获取模块，用于对地震记录S(xi，t)做时间方向的傅里叶变换，得到向量S(x，ω)，利用高精度Radon变换，得到地震数据的频率域Radon正变换数据Us(pk，ω)；
第一次下行波场频率域变换结果获取模块，用于对于任意频率分量ω及射线参数pk，给定虚反射与一次反射的振幅差异系数Rk＝-1，并计算得到虚反射与一次反射的延迟时间Δτk，构建出矩阵G(xi，pk，ω)，求解得到第一次的下行波场频率域Radon变换结果
Figure FDA0003665716570000041
预测参数获取模块，用于通过交替迭代反演目标函数，得到射线参数pk所对应的虚反射与一次反射的平均延迟时间
Figure FDA0003665716570000042
和平均振幅差异系数
Figure FDA0003665716570000043
下行波场获取模块，通过得到
Figure FDA0003665716570000044
和
Figure FDA0003665716570000045
计算下行波场
Figure FDA0003665716570000046
虚反射对应的频率域变换结果获取模块，用于再次计算虚反射与一次反射的平均延迟时间
Figure FDA0003665716570000047
和平均振幅差异系数
Figure FDA0003665716570000048
不断交替迭代直至满足收敛条件，得到最终的虚反射对应的频率域Radon变换结果UDown(p，ω)；
最终虚反射预测数据获取模块，用于对UDown(p，ω)进行频率域Radon反变换，并进行时间域傅里叶反变换，则得到最终的虚反射预测数据uDown(x，t)。
7.一种信息数据处理终端，其特征在于，所述信息数据处理终端包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行权利要求1-5任意一项所述二维海面鬼波水体成像测量方法。
8.一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行权利要求1-5任意一项所述二维海面鬼波水体成像测量方法。
9.一种测流设备ADCP，其特征在于，所述测流设备ADCP执行权利要求1-5任意一项所述二维海面鬼波水体成像测量方法。
10.一种存储在计算机可读介质上的计算机程序产品，包括计算机可读程序，供于电子装置上执行时，提供用户输入接口以实施权利要求1-5任意一项所述二维海面鬼波水体成像测量方法。
CN202210592015.6A 2022-05-27 2022-05-27 二维海面鬼波水体成像测量方法、系统、终端及测流设备 Pending CN115061197A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210592015.6A CN115061197A (zh) 2022-05-27 2022-05-27 二维海面鬼波水体成像测量方法、系统、终端及测流设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210592015.6A CN115061197A (zh) 2022-05-27 2022-05-27 二维海面鬼波水体成像测量方法、系统、终端及测流设备

Publications (1)

Publication Number Publication Date
CN115061197A true CN115061197A (zh) 2022-09-16

Family

ID=83198090

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210592015.6A Pending CN115061197A (zh) 2022-05-27 2022-05-27 二维海面鬼波水体成像测量方法、系统、终端及测流设备

Country Status (1)

Country Link
CN (1) CN115061197A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117148443A (zh) * 2023-10-27 2023-12-01 胜利信科(山东)勘察测绘有限公司 基于鬼波提取与转换的浅剖数据信噪比增强方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117148443A (zh) * 2023-10-27 2023-12-01 胜利信科(山东)勘察测绘有限公司 基于鬼波提取与转换的浅剖数据信噪比增强方法
CN117148443B (zh) * 2023-10-27 2024-03-19 胜利信科(山东)勘察测绘有限公司 基于鬼波提取与转换的浅剖数据信噪比增强方法

Similar Documents

Publication Publication Date Title
US11428834B2 (en) Processes and systems for generating a high-resolution velocity model of a subterranean formation using iterative full-waveform inversion
CN102298156B (zh) 用于反虚反射地震数据的方法和装置
Shin et al. Waveform inversion using a logarithmic wavefield
US7961551B2 (en) Determining directional propagation attributes of a seismic event
US9632194B2 (en) Estimating and using slowness vector attributes in connection with a multi-component seismic gather
MX2011006036A (es) Uso de inversion de forma de onda para determinar las propiedades de un medio en el subsuelo.
US8760967B2 (en) Generating an angle domain common image gather
US20170031045A1 (en) Method and apparatus for modeling and separation of primaries and multiples using multi-order green's function
MX2014003060A (es) Sistemas y metodos de filtracion de dominio de frecuencia y discriminacion de dominio de espacio-tiempo de datos sismicos.
EA032186B1 (ru) Сейсмическая адаптивная фокусировка
US10795039B2 (en) Generating pseudo pressure wavefields utilizing a warping attribute
CN111239827A (zh) 基于局部相似系数的三维地震数据多次波压制方法
US9310504B2 (en) Systems and methods for detecting swell noise in a seismic gather
CN115061197A (zh) 二维海面鬼波水体成像测量方法、系统、终端及测流设备
CN112462427B (zh) 多分量地震资料保幅角度域共成像点道集提取方法及系统
Guo et al. Becoming effective velocity-model builders and depth imagers, Part 2—The basics of velocity-model building, examples and discussions
GB2503640A (en) Quality Assurance in a Full Waveform Inversion Process
CN113866821B (zh) 一种基于照明方向约束的被动源干涉偏移成像方法和系统
Rabinovich et al. Location technology for construction of seismic images
Pinson et al. The image‐source method: a tool for geoacoustic inversion
CN112462428B (zh) 一种多分量地震资料偏移成像方法及系统
US20220390632A1 (en) Method and system for reflection-based travel time inversion using segment dynamic image warping
Wang et al. Converted-wave reflection traveltime inversion with free-surface multiples for ocean-bottom-node data
Thiel et al. 2d Acoustic Full Waveform Inversion of Submarine Salt Layer Using Dual Sensor Streamer Data
CN117518251A (zh) 自由界面多次波预测方法、装置、计算设备及存储介质

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