CN109541692A - 一种基于地震资料共振成像中的时频分析方法 - Google Patents
一种基于地震资料共振成像中的时频分析方法 Download PDFInfo
- Publication number
- CN109541692A CN109541692A CN201811635998.7A CN201811635998A CN109541692A CN 109541692 A CN109541692 A CN 109541692A CN 201811635998 A CN201811635998 A CN 201811635998A CN 109541692 A CN109541692 A CN 109541692A
- Authority
- CN
- China
- Prior art keywords
- data
- parameter
- length
- frequency
- resonance image
- 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
- 238000004458 analytical method Methods 0.000 title claims abstract description 23
- 238000000034 method Methods 0.000 claims abstract description 22
- 230000004888 barrier function Effects 0.000 claims description 9
- 102000008297 Nuclear Matrix-Associated Proteins Human genes 0.000 claims description 3
- 108010035916 Nuclear Matrix-Associated Proteins Proteins 0.000 claims description 3
- 210000000299 nuclear matrix Anatomy 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 abstract description 11
- 238000004422 calculation algorithm Methods 0.000 abstract description 8
- 238000004364 calculation method Methods 0.000 abstract description 4
- 238000005457 optimization Methods 0.000 abstract description 3
- 238000010183 spectrum analysis Methods 0.000 description 6
- 238000011160 research Methods 0.000 description 4
- 238000001228 spectrum Methods 0.000 description 4
- 230000009977 dual effect Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000002123 temporal effect Effects 0.000 description 3
- 230000009466 transformation Effects 0.000 description 3
- 238000005553 drilling Methods 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 238000002203 pretreatment Methods 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000002547 anomalous effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 238000007596 consolidation process Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000009659 non-destructive testing Methods 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
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
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
本发明涉及一种基于地震资料共振成像中的时频分析方法。该方法基于稀疏优化算法和Tikhonov正则化理论,通过引入拉格朗日乘子和控制步长,采用Newton迭代计算的方式获取频率系数的最优解。该方法提高了成像分析的速度、以及成像的分辨率和准确性,该方法能有效的应用于地下空间精细成像中。
Description
技术领域
本发明涉及勘探技术领域,具体地说,涉及一种基于地震资料共振成像中的时频分析方法。
背景技术
自然界中,任何物质都有其自身的固有频率,包括地下各地质体。当有一个宽频的震动传播到该地质体,特征固有频率能量将被放大,通过观测被放大的特征频率信号,对特征频率信号成像获得地下精细成像效果。目前,共振成像法应用广泛,例如大型土木工程方面的查明地下介质结构并进行分层,评价地基的加固效果,建筑工程的抗震设计,公路、机场跑道的质量无损检测;环境地质方面的勘查岩溶洞穴的分布范围,滑坡活动面等地质灾害和环境地质的调查;地下空间精细探测;以及考古发现等。
基于地震频谱的动力学特征,1989年,Nakamura首次提出了基于地震资料的单点谱比法。该方法其中的关键技术之一就是地震数据的频谱分析。常用的地震数据频谱分析方法主要包括傅里叶变换(FFT)、短时傅里叶变换(STFT)、连续小波分析(CWT)、匹配追踪算法(MP)、以及魏格纳分布(WVD)等。其中,FFT是最早提出并用于信号频谱分析的技术,但由于其是对数据进行整体分析,缺乏时间分辨率,不适用于非平稳的地震信号分析;Potter(1947)对FFT加上窗函数,首次提出STFT的概念,进而能获得时间-频域信息,但由于每次计算的窗函数固定,不能兼顾时间分辨率和频率分辨率;Morlet(1982)提出的CWT是一种多分辨率的时频分析方法,但是计算量大,耗时长;Mallet等(1999)提出的MP方法是借助子波字典获得由匹配参数构成的视频谱,该方法十分依赖于所选择的子波字典,并且计算量大;WVD方法弱化窗函数的影响,提高了时间分辨率,但由于其是非线性变换会造成交叉项干扰是的频谱失真。
发明内容
在地震数据时频分析过程中,噪音干扰,资料品质差、先验信息缺乏、数据量不足等条件都会导致资料时频分析反演时出现解不存在、解不唯一、或者解不稳定等问题。为了克服地震面波数据时频分析反问题中可能出现的不适当性、计算过程的病态性以及大规模计算等问题。提出了一种基于稀疏优化算法和正则化理论的时频分析中地震频率反问题求解算法。
基于地震资料共振成像中的时频分析方法的具体步骤如下:
步骤1:设置初始频率系数f1,f1>0;迭代容忍度εf>0,ελ>0,εμ>0,初始控制步长τ1=1,迭代步长τ的修正因子θ∈(0,1),衰减参数β∈[0,1],迭代次数阈值imax=30;
步骤2:初始化(f1,λ1,s1)T,其中,拉格朗日乘子λ1>0,中间参数s1>0,迭代次数i=1;
步骤3:计算余量:其中M为实数或者复数正弦基函数的核矩阵,d为带噪音的地震面波数据,N是数据d的长度,e为一个所有分量都为1的向量,μ为障碍函数参数,若或i=imax,结束迭代,否则,进入步骤4;
步骤4:取衰减参数β∈(0,1],计算参数矫正量(Δfi,Δλi,Δsi);
步骤5:计算取θ∈(0,1),令τi+1=min{θτmax,1},更新迭代,i=i+1,返回步骤3,继续判断。
其中,步骤2中所述的拉格朗日乘子满足以下条件,其中,μ为障碍函数参数。
其中,障碍函数参数μ满足以下条件,
其中,步骤2中所述的中间参数s=μD―1e,其中μ为障碍函数参数,d为带噪音的地震面波数据,N是数据d的长度,e为一个所有分量都为1的向量。
其中,步骤计4中所述的参数矫正量(Δfi,Δλi,Δsi)由以下公式获得其中
该方法基于稀疏优化算法和Tikhonov正则化理论,通过引入拉格朗日乘子和控制步长,采用Newton迭代计算的方式获取频率系数的最优解。该方法提高了成像分析的速度、以及成像的分辨率和准确性,该方法能有效的应用于地下空间精细成像中。
附图说明
图1:模拟数据分析图
图2:预处理后数据整体分布图
图3:桩号34点预处理后数据分布图
图4:桩号34点资料的频谱分析结果
图5:共振成像剖面与钻井资料的对比分析图
图6:共振成像剖面
图7:共振成像剖面的局部分析图(管道)
图8:共振成像剖面的局部分析图(塌陷)
具体实施方式
我们可以把地震面波时频分析问题归结为如下方程形式:
Mf=d(1)
其中,f表示频率系数,d表示带噪音的地震面波数据,M表示实数或者复数正弦基函数的核矩阵,
M(t,x)=cos(2π·k·Δx·n·Δt)+isin(2π·k·Δx·n·Δt) (2)
其中,k表示频率域采样指数,Δx表示频率增量,n表示时间采样指数,Δt表示时间增量。
对于上式的求解一般可通过如下公式(3)、(4)求得:
M*Mf=M*f (3)
f=(M*M)―1M*d (4)
如果M为一正交矩阵,则M*M=I,这时公式(1)的最小二乘解可由f=M*d获得,即离散傅里叶变换(DFT)。但是,在实际地震资料中,M的元素通常是相关的,M具有非正交性。这将导致公式(1)解的不唯一,以及直接求解公式(4)的病态性。
考虑地震观测系统的限制,地球物理反问题一般是离散的、稀疏的,以及频率系数f的非负性等特征。为此,本文利用Tikhonov正则化思想,对地震时频分析中求解频率f的反问题建立如下基于l1范数的稀疏优化模型:
公式(5)的原始对偶形式为
min eTf,s.t.Mf=d,f≥0 (6)
其中,f有N个分量,f=[f1,f2,…fN];N是数据d的长度;e是一个所有分量都为1的向量。
利用障碍函数(μ>0)来求解公式(6)的原始对偶问题:
公式(7)中,当f<0时,表达式无意义;也就是借助障碍函数将不等式约束(公式6)转化为无约束的最小化问题(公式7)。当对偶间隙μ→0,公式(7)的解会不断逼近原问题的全局最优解。
引入拉格朗日乘子λ,利用拉格朗日乘子法对公式(7)的最小化问题进行求解,
由于
gradfL(f,λ)=e―μD―1e-MTλ (9)
gradλL(f,λ)=d―Mf (10)
令s=μD―1e,可得到公式(7)的最优性条件:
s+MTλ=e,Mf=d,Ds=μe,f>0 (11)
其中,
该系统的一个最优解(f,λ,s)T即是原始的可行解又是对偶可行的。在线性规划领域,向量(f,λ,s)T称为中心路径。
考虑到对偶间隙
eTf―dTλ=fTs=μN (13)
因此,数值上,求解非线性方程组
其中,β是衰减参数,β∈[0,1],控制下降速度使得迭代点的轨迹在中心路径附近。
借助Newton迭代法对上式进行求解,假设第i次迭代的(fi,λi,si)T给定,对公式(14)求导,其梯度如下:
其中,
Newton迭代方法中的参数校正量(Δfi,Δλi,Δsi)T可由以下公式获得:
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
利用图1中的模拟信号进行测试分析,令εf=0.0001,ελ=0.0001,εμ=0.00001,频率系数均为0.1,μ=4,θ=0.85,β=0.9,s1可由s=μD―1e获得,λ1可由公式(9)获得。对比常规STFT,本算法在不降低计算效率的情况下,能清晰的区分10-15Hz段中的11.108Hz、12.573Hz、13.062Hz三个频率信息,改善了资料的频域分辨率,证实算法的有效性。进一步,根据浅层勘探(地下200米以上)中共振频率-深度的经验关系,,A为一常数,取值范围是50-300。以A=180为例,三个频率分别对应地下8.499米、7.264米、6.921米,也因此表明该算法能清晰区分、并刻画此地下三个不同深度的异常体,证实了该算法在浅层勘探共振成像中的有效性。
研究区位于成都市空港区牧华路沿线,共振激发方式为天然源,采用的是高精度2.5HZ三分量低频检波器,采集时间长度为18分钟,采样率为4毫秒,图2预处理后数据整体分布图,图3是桩号34点预处理后数据分布图。
同样,在本研究区的数据测试中,令εf=0.0001,ελ=0.0001,εμ=0.00001,频率系数均为0.1,μ=4,θ=0.85,β=0.9,s1可由s=μD―1e获得,λ1可由公式(9)获得。对预处理之后的数据,进行逐点单分量的频谱分析。图4是桩号34点处地震资料的频谱分析结果。由图5分析可见,本次研究的共振成像剖面与相应的钻井资料对比,成像深度误差小于0.5米,证实了本方法结果的可靠性与准确性。
从成都空港区的共振成像剖面(图6)可见清晰的鹅卵石层、粘土层、强风化基岩和弱风化基岩等地质现象,为该研究区地下地质空间建模提供高精度成像。同时,从共振成像局部剖面(图7、8)的精细对比分析可见沿测线的下水管道、雨水管道、塌陷等,成像刻画清晰,且成像分辨率较高,这可为城市地下管线探测,特别是深埋管线探测,以及道路灾害探测,提供高精度成像。此外,从分析结果可推断地面下陷原因与下覆雨水管道掏空造成地层疏松有关,这可以为城市建设规划提供重要的参考。
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解,本发明不受上述实施例的限制,上述实施例和说明书中描述的只是说明本发明的原理,在不脱离本发明精神和范围的前提下,本发明还会有各种变化和改进,这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
Claims (5)
1.一种基于地震资料共振成像中的时频分析方法,其特征在于,具体包括如下步骤:
步骤1:设置初始频率系数f1,f1>0;迭代容忍度εf>0,ελ>0,εμ>0,初始控制步长τ1=1,迭代步长τ的修正因子θ∈(0,1),衰减参数β∈[0,1],迭代次数阈值imax=30;
步骤2:初始化(f1,λ1,s1)T,其中,拉格朗日乘子λ1>0,中间参数s1>0,迭代次数i=1;
步骤3:计算余量:其中M为实数或者复数正弦基函数的核矩阵,d为带噪音的地震面波数据,N是数据d的长度,e为一个所有分量都为1的向量,μ为障碍函数参数,若或i=imax,结束迭代,否则,进入步骤4;
步骤4:取衰减参数β∈(0,1],计算参数矫正量(Δfi,Δλi,Δsi);
步骤5:计算取θ∈(0,1),令τi+1=min{θτmax,1},更新迭代,i=i+1,返回步骤3,继续判断。
2.根据权利要求1所述的方法,其特征在于,所述步骤2中所述的拉格朗日乘子满足以下条件,其中,μ为障碍函数参数。
3.根据权利要求2所述的方法,其特征在于,所述障碍函数参数μ满足以下条件,s.t.Mf=d。
4.根据权利要求1所述的方法,其特征在于,所述步骤2中所述的中间参数s=μD-1e,其中μ为障碍函数参数,d为带噪音的地震面波数据,N是数据d的长度,e为一个所有分量都为1的向量。
5.根据权利要求1所述的方法,其特征在于,所述步骤计4中所述的参数矫正量(Δfi,Δλi,Δsi)由以下公式获得其中
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811635998.7A CN109541692B (zh) | 2018-12-29 | 2018-12-29 | 一种基于地震资料共振成像中的时频分析方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811635998.7A CN109541692B (zh) | 2018-12-29 | 2018-12-29 | 一种基于地震资料共振成像中的时频分析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109541692A true CN109541692A (zh) | 2019-03-29 |
CN109541692B CN109541692B (zh) | 2020-06-05 |
Family
ID=65831223
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811635998.7A Active CN109541692B (zh) | 2018-12-29 | 2018-12-29 | 一种基于地震资料共振成像中的时频分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109541692B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113759422A (zh) * | 2021-09-10 | 2021-12-07 | 刘俊伟 | 一种地下异常体探测方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20040186667A1 (en) * | 2003-03-18 | 2004-09-23 | Lee Ki Ha | Source-independent full waveform inversion of seismic data |
CN102077108A (zh) * | 2008-04-28 | 2011-05-25 | 康奈尔大学 | 分子mri中的磁敏度精确量化 |
CN106291700A (zh) * | 2016-09-28 | 2017-01-04 | 西安交通大学 | 基于同步挤压变换的地震加权平均瞬时频率提取方法 |
CN106918838A (zh) * | 2017-03-06 | 2017-07-04 | 中国科学院地质与地球物理研究所 | 起伏地表条件下高斯束偏移成像方法及装置 |
CN107728211A (zh) * | 2017-08-31 | 2018-02-23 | 电子科技大学 | 基于张量核范数正则化的地震信号恢复算法 |
-
2018
- 2018-12-29 CN CN201811635998.7A patent/CN109541692B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20040186667A1 (en) * | 2003-03-18 | 2004-09-23 | Lee Ki Ha | Source-independent full waveform inversion of seismic data |
CN102077108A (zh) * | 2008-04-28 | 2011-05-25 | 康奈尔大学 | 分子mri中的磁敏度精确量化 |
CN106291700A (zh) * | 2016-09-28 | 2017-01-04 | 西安交通大学 | 基于同步挤压变换的地震加权平均瞬时频率提取方法 |
CN106918838A (zh) * | 2017-03-06 | 2017-07-04 | 中国科学院地质与地球物理研究所 | 起伏地表条件下高斯束偏移成像方法及装置 |
CN107728211A (zh) * | 2017-08-31 | 2018-02-23 | 电子科技大学 | 基于张量核范数正则化的地震信号恢复算法 |
Non-Patent Citations (1)
Title |
---|
王彦飞 等: ""地震时频分析的加权l1范数稀疏正则化及交替方向乘子算"", 《中国科学:数学》 * |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113759422A (zh) * | 2021-09-10 | 2021-12-07 | 刘俊伟 | 一种地下异常体探测方法 |
CN113759422B (zh) * | 2021-09-10 | 2023-11-21 | 刘俊伟 | 一种地下异常体探测方法 |
Also Published As
Publication number | Publication date |
---|---|
CN109541692B (zh) | 2020-06-05 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Mi et al. | Horizontal resolution of multichannel analysis of surface waves | |
CN103454685B (zh) | 利用测井约束波阻抗反演预测砂体厚度的方法和装置 | |
CN104977618A (zh) | 一种评价页岩气储层及寻找甜点区的方法 | |
CN108152854B (zh) | 一种基于微动功率谱密度的无损探测方法及其应用 | |
CN104853822A (zh) | 一种评价页岩气储层及寻找甜点区的方法 | |
CN109490963B (zh) | 裂缝储层岩石物理建模方法及系统 | |
Molnar et al. | Overview of local site effects and seismic microzonation mapping in Metropolitan Vancouver, British Columbia, Canada | |
CN106154323A (zh) | 基于地震拓频处理的相控随机反演薄储层预测方法 | |
CN105277982A (zh) | 一种泥页岩总有机碳含量地震预测方法 | |
Liang et al. | The application of HVSR method in detecting sediment thickness in karst collapse area of Pearl River Delta, China | |
Xue et al. | Extracting the virtual reflected wavelet from TEM data based on regularizing method | |
Montone et al. | Constraints on the structure of the shallow crust in central Italy from geophysical log data | |
Passeri et al. | The Polito Surface Wave flat-file Database (PSWD): statistical properties of test results and some inter-method comparisons | |
EP4089446A1 (en) | A method for seismic frequency resonance exploration technology | |
Xiaoxia et al. | Multi-scale fracture prediction and characterization method of a fractured carbonate reservoir | |
CN109541692A (zh) | 一种基于地震资料共振成像中的时频分析方法 | |
Sisman et al. | Evaluation of site response with alternative methods: A case study for engineering implications | |
Peng et al. | Joint tomography of multi-cross-hole and borehole-to-surface seismic data for karst detection | |
Mohamed et al. | Near-surface site characterization at Quriyat City, Sultanate of Oman using HVSR and MASW techniques | |
Cheng et al. | Experimental study of small fixed-loop transient electromagnetic method for characterizing water-bearing structures in tunnels | |
CN104345347B (zh) | 一种用于致密含气砂岩储层预测的测井曲线恢复方法 | |
Dong et al. | Diagnosis of concentrated leakage channel embedded in dam base by means of hydraulic tomography | |
Danilov et al. | Study of deep structure of the kimberlite pipe named after M. Lomonosov of the Arkhangelsk diamondiferous province obtained by joint using of passive seismic and radiometric methods | |
Lou et al. | A new method and application of groundwater prediction along the direction of tunnel excavation in karst strata | |
Lu et al. | A rapid imaging method of the seismic back‐scattered wavefield for urban road near‐surface anomalous structures |
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 |