CN104730572A

CN104730572A - A Diffraction Wave Imaging Method and Device Based on L0 Half-Norm

Info

Publication number: CN104730572A
Application number: CN201510107121.0A
Authority: CN
Inventors: 于彩霞; 王彦飞
Original assignee: Institute of Geology and Geophysics of CAS
Current assignee: Institute of Geology and Geophysics of CAS
Priority date: 2015-03-11
Filing date: 2015-03-11
Publication date: 2015-06-24
Anticipated expiration: 2035-03-11
Also published as: CN104730572B

Abstract

The invention discloses a diffraction wave imaging method based on the L0 semi-norm, which acquires seismic data from which reflected waves have been removed as input data; discretizes the underground imaging space, and selects any imaging point as the current diffraction imaging point, and according to the relationship between the shot point and the receiver point in the input seismic data and the given velocity model, calculate the Green's function of any diffraction imaging point above; execute other imaging points that are not selected in a loop, and calculate their corresponding Green's function, Until the Green's function of all imaging points in the underground imaging space is obtained; a diffraction wave imaging model based on the L0 semi-norm is constructed; the model is solved by the homotopy analysis iterative algorithm, and the diffraction wave imaging result is obtained. The invention also discloses a diffraction wave imaging device based on the L0 half norm. The invention can improve the imaging resolution of seismic data and enhance the signal-to-noise ratio of diffracted wave imaging, so that it is easier to identify small-scale geological units related to the spatial connectivity of reservoirs.

Description

A Diffraction Wave Imaging Method and Device Based on L0 Half-Norm

技术领域technical field

本发明属于勘探地震技术领域，涉及一种基于L0半范数的绕射波成像方法，本发明还涉及一种基于L0半范数的绕射波成像装置。The invention belongs to the field of exploration seismic technology, and relates to a diffraction wave imaging method based on the L0 half norm, and also relates to a diffraction wave imaging device based on the L0 half norm.

背景技术Background technique

在油气资源勘探开发过程中，有效识别构造和岩性异常体，如溶洞、裂缝、地层尖灭、风化壳等对储层认识至关重要。上述小尺度地质体在空间展布上通常小于或远小于地震子波波长，因此在地震数据中以绕射波的形式存在。在野外作业过程中，检波器接收到的地震信号包含反射波与绕射波，而反射波成像条件假定无限大光滑界面，因此一般只能反映大尺度地质背景。相对而言，绕射波则是地质细节的反映，是提高地震分辨率的重要信息载体。但地质体产生的绕射能量约为反射能量0.1-0.01倍，为弱信号，常淹没在反射波成像结果中，因此需要去除反射波，常用的方法包括信号分解的有Harlan变换(Harlan et at.,1984)、Radon变换(Zhang,2005)、平面波破坏滤波(Taner et al.,2006；Fomel et al.,2006,2007)等。In the process of exploration and development of oil and gas resources, effective identification of structural and lithological anomalies, such as karst caves, fractures, formation pinch-outs, and weathering crusts, is crucial to reservoir understanding. The above-mentioned small-scale geological bodies are usually smaller or much smaller than the wavelength of the seismic wavelet in terms of spatial distribution, so they exist in the form of diffracted waves in seismic data. During field operations, the seismic signals received by the geophone include reflected waves and diffracted waves, and the reflected wave imaging condition assumes an infinite smooth interface, so generally it can only reflect large-scale geological backgrounds. Relatively speaking, diffracted waves reflect geological details and are an important information carrier for improving seismic resolution. However, the diffraction energy produced by the geological body is about 0.1-0.01 times the reflected energy, which is a weak signal and is often submerged in the reflected wave imaging results. Therefore, it is necessary to remove the reflected wave. The commonly used methods include Harlan transform (Harlan et at .,1984), Radon transform (Zhang,2005), plane wave destruction filter (Taner et al.,2006; Fomel et al.,2006,2007), etc.

经专利检索及国内外文献调研，发现目前开展的绕射波成像方法，焦点多数放在反射波去除上，目的是通过压制反射波以突出绕射波。如偏移-倾角道集反射波去除法(Klokov and Fomel,2012)、聚焦-切除-反聚焦(Khaidukov et al.,2004)、共反射面元叠加(Dell and Gajewski,2011；Asgedomet al.,2011)、多聚焦(Berkovitch et al.,2009)等。After patent retrieval and domestic and foreign literature research, it is found that most of the current diffraction wave imaging methods focus on the removal of reflected waves, and the purpose is to highlight the diffracted waves by suppressing the reflected waves. Such as offset-dip gather reflection wave removal method (Klokov and Fomel, 2012), focus-cut-defocus (Khaidukov et al., 2004), common reflection surface element superposition (Dell and Gajewski, 2011; Asgedometal., 2011), multi-focus (Berkovitch et al., 2009), etc.

上述绕射波成像方法没有针对绕射波成像算子展开研究，实际上，绕射信息在空间分布上具有稀疏不连续性，完全可以利用数学上的L0半范数描述。The above-mentioned diffraction wave imaging method does not study the diffraction wave imaging operator. In fact, the diffraction information has sparse discontinuity in the spatial distribution, which can be described by the mathematical L0 semi-norm.

因此，本发明根据格林函数构建了基于L0半范数的绕射波成像模型，该模型通过约束绕射信息的稀疏性，可以达到高精度、高信噪比绕射波成像结果。与一般绕射波成像方法相比，基于L0半范数的绕射波成像方法，可以在提高地震信号分辨率的同时，进一步保证成像资料的信噪比，并减小多解性。Therefore, the present invention constructs a diffraction wave imaging model based on the L0 semi-norm based on the Green's function, and the model can achieve high precision and high signal-to-noise ratio diffraction wave imaging results by constraining the sparsity of diffraction information. Compared with the general diffracted wave imaging method, the diffracted wave imaging method based on the L0 half-norm can improve the resolution of seismic signals while further ensuring the signal-to-noise ratio of imaging data and reducing multi-solution.

发明内容Contents of the invention

本发明的目的是提供一种基于L0半范数的绕射波成像方法，解决现存绕射波成像技术中，在提高成像分辨率的同时无法兼顾地震资料信噪比的问题，该技术能够实现小断层、裂缝、溶洞等小尺度地质异常体精细刻画。The purpose of the present invention is to provide a diffraction wave imaging method based on the L0 semi-norm, which solves the problem that the existing diffraction wave imaging technology cannot take into account the signal-to-noise ratio of seismic data while improving the imaging resolution. Fine depiction of small-scale geological anomalies such as small faults, cracks, and caves.

本发明的另一目的是提供一种基于L0半范数的绕射波成像装置。Another object of the present invention is to provide a diffracted wave imaging device based on the L0 half-norm.

本发明所采用的技术方案是，一种基于L0半范数的绕射波成像方法，包括以下步骤：The technical solution adopted in the present invention is, a kind of diffracted wave imaging method based on L0 semi-norm, comprises the following steps:

步骤101：获取一个去除反射波的地震数据，作为输入数据；Step 101: Obtain a seismic data from which reflection waves have been removed as input data;

步骤102：对地下成像空间进行离散化，选取任意一个成像点作为当前绕射成像点，并依据输入地震数据中的炮点与检波点关系及给定的速度模型，计算上述任意绕射成像点的格林函数；Step 102: Discretize the underground imaging space, select any imaging point as the current diffraction imaging point, and calculate the above-mentioned arbitrary diffraction imaging point according to the relationship between the shot point and the receiver point in the input seismic data and the given velocity model Green's function;

步骤103：循环执行上述步骤102中地下成像空间中未选取的其它成像点，并计算其相应的格林函数，直至得出地下成像空间所有成像点格林函数；Step 103: cyclically execute other imaging points not selected in the underground imaging space in the above step 102, and calculate their corresponding Green's functions until the Green's functions of all imaging points in the underground imaging space are obtained;

步骤104：根据地下成像空间所有成像点的格林函数与初始成像模型，构建基于L0半范数的绕射波成像模型；Step 104: Construct a diffracted wave imaging model based on the L0 semi-norm according to the Green's function of all imaging points in the underground imaging space and the initial imaging model;

步骤105：由同伦分析迭代算法求解该模型，得出绕射波成像结果。Step 105: The model is solved by the homotopy analysis iterative algorithm, and the diffraction wave imaging result is obtained.

本发明所采用的另一技术方案是，一种基于L0半范数的绕射波成像装置，包括依次连接的地震数据获取模块，用于输入一个去除反射波的地震数据；绕射成像点选取模块，用于从离散化的地下成像空间选取成像点作为格林函数计算位置；格林函数计算模块，用于计算由炮点出发经绕射成像点到地震检波点的走时及振幅补偿项；L0半范数模型构建模块，依据地下空间所有成像点格林函数构建基于L0半范数的求解模型；模型求解器件，利用同伦分析算法求解由格林函数与绕射点模型构建的L0模型，得出绕射波成像结果。Another technical solution adopted in the present invention is, a kind of diffraction wave imaging device based on the L0 semi-norm, comprising sequentially connected seismic data acquisition modules, used to input a seismic data that removes reflected waves; diffraction imaging point selection The module is used to select the imaging point from the discretized underground imaging space as the Green’s function to calculate the position; the Green’s function calculation module is used to calculate the travel time and amplitude compensation item from the shot point to the seismic receiver point through the diffraction imaging point; L0 half The norm model building module constructs a solution model based on the L0 semi-norm based on the Green's function of all imaging points in the underground space; the model solving device uses the homotopy analysis algorithm to solve the L0 model constructed by the Green's function and the diffraction point model, and obtains the Radiographic imaging results.

本发明的有益效果是，对去除反射波的地震数据，进行格林函数计算，从而构建L0半范数反演模型。通常情况下，反演问题求解是非唯一的，除非加上限定条件以缩小搜索范围，考虑到绕射信息空间不连续性，在实施例中，限定了模型稀疏性，使得求解模型更为符合实际地质情况，并且通过迭代逼近可以在一定程度上压制噪声，提高绕射波成像资料信噪比。提高地震资料成像分辨率，增强绕射波成像信噪比，从而更容易识别与储层空间连通性有关的小尺度地质单元。The beneficial effect of the invention is that the Green's function calculation is performed on the seismic data from which reflected waves have been removed, so as to construct an L0 semi-norm inversion model. Usually, the solution to the inversion problem is non-unique, unless a limited condition is added to narrow the search range. Considering the discontinuity of the diffraction information space, in the embodiment, the model sparsity is limited to make the solution model more realistic Geological conditions, and through iterative approximation, the noise can be suppressed to a certain extent, and the signal-to-noise ratio of diffraction wave imaging data can be improved. Improving the imaging resolution of seismic data and enhancing the signal-to-noise ratio of diffraction wave imaging makes it easier to identify small-scale geological units related to reservoir spatial connectivity.

附图说明Description of drawings

图1是本发明方法的流程图。Figure 1 is a flow chart of the method of the present invention.

图2是本发明装置的结构框图。Fig. 2 is a structural block diagram of the device of the present invention.

其中，201.地震数据获取模块，202.绕射成像点选取模块，203.格林函数计算模块，204.L0半范数模型构建模块，205.模型求解器件。Among them, 201. Seismic data acquisition module, 202. Diffraction imaging point selection module, 203. Green's function calculation module, 204. L0 semi-norm model building module, 205. Model solving device.

具体实施方式Detailed ways

下面结合附图和具体实施方式对本发明进行详细说明。The present invention will be described in detail below in conjunction with the accompanying drawings and specific embodiments.

如图1所示，本发明提出的一种基于L0半范数的绕射波成像方法，包括以下步骤：As shown in Figure 1, a kind of diffracted wave imaging method based on L0 semi-norm that the present invention proposes comprises the following steps:

步骤101：获取一个去除反射波的地震数据，作为输入数据；相应的地震数据，可为叠前炮集地震数据或叠后地震数据，在反射波去除方法上，采用平面波破坏滤波方法；Step 101: Obtain a seismic data from which reflected waves have been removed as input data; the corresponding seismic data can be pre-stack shot seismic data or post-stack seismic data, and the method of removing reflected waves uses a plane wave damage filtering method;

步骤102：对地下成像空间进行离散化，并选取任意一个成像点作为当前绕射成像点，并依据输入地震数据中的炮点与检波点关系及给定的速度模型，计算上述任意绕射成像点的格林函数；Step 102: Discretize the underground imaging space, select any imaging point as the current diffraction imaging point, and calculate the above-mentioned arbitrary diffraction imaging according to the relationship between the shot point and the receiver point in the input seismic data and the given velocity model Green's function of the point;

在上述步骤102、103中的格林函数计算方法，依据射线追踪计算由炮点出发经成像点到检波点走时，并存储振幅加权项。In the Green's function calculation method in the above steps 102 and 103, the travel time from the shot point through the imaging point to the receiver point is calculated according to ray tracing, and the amplitude weighted items are stored.

在上述步骤104中，建立基于L0半范数的绕射波成像模型，设计如下：In the above step 104, a diffracted wave imaging model based on the L0 semi-norm is established, and the design is as follows:

$min min {J J}_{α α} ((m m)) : : = = \frac{11}{22} {| | | | Gm G m - - d d | | | |}_{{l l}_{22}}^{22} + + α α {| | | | m m | | | |}_{{l l}_{00}}$

其中，min表示最小化，J_α(m)为目标函数，m为求解的绕射模型，数学符号:＝表示定义为，G为格林函数，d为去除反射波的地震数据，α为正Wherein, min represents minimization, J _α (m) is objective function, m is the diffraction model of solution, mathematical symbol: = represents definition as, G is Green's function, d is the seismic data that removes reflected wave, α is positive

则化因子，表示l₂范数，为l₀半范数。regularization factor, Indicates the l ₂ norm, is the l ₀ semi-norm.

通过邻近点方法，可将上式写为：By the neighboring point method, the above formula can be written as:

$min min {H h}_{β β,, α α} (({m m}_{00},, m m)) : : = = (({G G}^{T T} {Gm G m}_{00} - - {G G}^{T T} d d + + α α \frac{d d}{dm dm} {| | | | m m | | | |}_{{l l}_{00},, m m = = {m m}_{00}},, m m - - {m m}_{00})) + + \frac{β β}{22} {| | | | m m - - {m m}_{00} | | | |}_{{l l}_{22}}^{22} + + α α {| | | | m m | | | |}_{{l l}_{00}}$

其中，H_β,α(m₀,m)为J_α(m)逼近式，m₀为初始模型，β为可调节的正则参数，(·,·)表示内积。Among them, H _β,α (m ₀ ,m) is the approximation formula of J _α (m), m ₀ is the initial model, β is the adjustable regular parameter, and (·,·) represent the inner product.

上式最小值问题，可通过硬阈值算子求得，如下所示：The minimum value problem of the above formula can be obtained through the hard threshold operator, as follows:

${V V}_{β β} (({m m}_{00})) = = \underset{m m}{arg arg min min} {H h}_{β β,, α α} (({m m}_{00},, m m))$

其中， $s_{β} (m) = m - \frac{1}{2 β} &dtri; {| | Gm - d | |}_{l_{2}}^{2},$ [·]_i表示矢量第_i个元素。in, ${the s}_{β} (m) = m - \frac{1}{2 β} &dtri; {| | G m - d | |}_{l_{2}}^{2},$ [·] _i represents _{the i} -th element of the vector.

上述步骤105中，对L0半范数模型构建模块，采取模型求解器件实现，其迭代过程由同伦分析算法完成，具体过程如下：In the above step 105, the construction module of the L0 semi-norm model is implemented by a model solving device, and the iterative process is completed by a homotopy analysis algorithm. The specific process is as follows:

步骤1：输入李普希茨参数β₀，正则化参数α₀，且要求Step 1: Input Lipschitz parameter β ₀ , regularization parameter α ₀ , and require

β₀∈[β_min,β_max]，β_min,β_max分别为普希茨参数β₀上下限，初始化参数，k＝0,ρ∈(0,1)，初始化模型m⁰；β ₀ ∈[β _min , β _max ], β _min , β _max are the upper and lower limits of Puschitz parameter β ₀ , initialization parameters, k=0, ρ∈(0,1), initialization model m ⁰ ;

步骤2:设定，i＝0,m^k,0＝m^k,β_k,0＝β_k；Step 2: setting, i=0, m ^{k, 0} = m ^k , β _{k, 0} = β _k ;

步骤3： $m^{k, i + 1} = V_{β_{k}, i} (m^{k, i + 1});$ Step 3: $m^{k, i + 1} = V_{β_{k}, i} (m^{k, i + 1});$

当 $J_{α_{k}} (m^{k, i}) - J_{α_{k}} (m^{k, i + 1}) < \frac{η}{2} | | m^{k, i} - m^{k, i + 1} | |,$ when $J_{α_{k}} (m^{k, i}) - J_{α_{k}} (m^{k, i + 1}) < \frac{η}{2} | | m^{k, i} - m^{k, i + 1} | |,$

执行β_k,i＝min{γβ_k,i,β_max}， Execute β _k,i =min{γβ _k,i ,β _max },

β_k,i+1＝β_k,i,i＝i+1；β _k,i+1 =β _k,i ,i=i+1;

返回步骤3，直至，||m^k,i-m^k,i+1||_∞≤ε₀ Return to step 3 until ||m ^k,i -m ^k,i+1 || _∞ ≤ε ₀

其中，终止参数ε₀一般为ε₀＝10^-1；Among them, the termination parameter ε ₀ is generally ε ₀ =10 ^-1 ;

步骤4：m^k+1＝m^k,i,β_k+1＝β_k,i,α_k+1＝ρα_k,k＝k+1；Step 4: m ^k+1 =m ^k,i ,β _k+1 =β _k,i ,α _k+1 =ρα _k ,k=k+1;

返回步骤2，直至，||m^k,i-m^k,i+1||_∞≤ε；Return to step 2 until ||m ^k,i -m ^k,i+1 || _∞ ≤ε;

其中，ε一般为ε＝10^-5；Among them, ε is generally ε=10 ^-5 ;

步骤5：输出最终迭代结果，m^*＝m^k。Step 5: output the final iteration result, m ^* = m ^k .

基于同一发明构思，本发明还提供了一种基于L0半范数的绕射波成像装置。由于一种基于L0半范数的绕射波成像装置解决问题的原理与一种基于L0半范数的绕射波成像方法相似，因此一种基于L0半范数的绕射波成像装置的实施可以参见一种基于L0半范数的绕射波成像方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。Based on the same inventive concept, the present invention also provides a diffracted wave imaging device based on the L0 half-norm. Since the problem-solving principle of a diffraction wave imaging device based on the L0 half-norm is similar to that of a diffraction wave imaging method based on the L0 half-norm, the implementation of a diffraction wave imaging device based on the L0 half-norm Reference may be made to the implementation of a diffracted wave imaging method based on the L0 half-norm, and repeated descriptions will not be repeated here. As used below, the term "unit" or "module" may be a combination of software and/or hardware that realizes a predetermined function. Although the devices described in the following embodiments are preferably implemented in software, implementations in hardware, or a combination of software and hardware are also possible and contemplated.

图2是本发明的一种基于L0半范数的绕射波成像装置的一种结构框图，包括依次连接的地震数据获取模块201、绕射成像点选取模块202、格林函数计算模块203、L0范数模型构建模块204与模型求解器件205，下面对该结构进行说明。Fig. 2 is a kind of structural block diagram of a kind of diffracted wave imaging device based on L0 semi-norm of the present invention, comprises sequentially connected seismic data acquisition module 201, diffraction imaging point selection module 202, Green's function calculation module 203, L0 The structure of the norm model building module 204 and the model solving device 205 will be described below.

地震数据获取模块201，用于输入一个去除反射波的地震数据；A seismic data acquisition module 201, configured to input seismic data from which reflected waves have been removed;

绕射成像点选取模块202，用于从离散化的地下成像空间选取成像点作为格林函数计算位置；The diffraction imaging point selection module 202 is used to select the imaging point from the discretized underground imaging space as the Green's function calculation position;

格林函数计算模块203，用于计算由炮点出发经绕射成像点到地震检波点的走时及振幅补偿项；The Green's function calculation module 203 is used to calculate the travel time and amplitude compensation items from the shot point to the seismic receiver point through the diffraction imaging point;

L0半范数模型构建模块204，依据地下空间所有成像点格林函数构建基于L0半范数的求解模型；The L0 semi-norm model construction module 204, constructs a solution model based on the L0 semi-norm according to the Green's function of all imaging points in the underground space;

模型求解器件205，利用同伦分析算法求解由格林函数与绕射点模型构建的L0模型，得出绕射波成像结果。The model solving device 205 uses the homotopy analysis algorithm to solve the L0 model constructed by the Green's function and the diffraction point model to obtain the diffraction wave imaging result.

其中，格林函数计算模块203，依据射线追踪计算由炮点出发经成像点到检波点走时，并存储振幅加权项。Wherein, the Green's function calculation module 203 calculates the travel time from the shot point through the imaging point to the receiver point according to ray tracing, and stores the amplitude weighted items.

L0半范数模型构建模块204，设计如下：The L0 semi-norm model building block 204 is designed as follows:

其中，min表示最小化，J_α(m)为目标函数，m为求解的绕射模型，数学符号:＝表示定义为，G为格林函数，d为去除反射波的地震数据，α为正则化因子，表示l₂范数，为l₀半范数。Among them, min represents the minimum, J _α (m) is the objective function, m is the diffraction model to be solved, the mathematical symbol: = represents the definition as, G is the Green's function, d is the seismic data that removes the reflected wave, and α is the regularization factor, Indicates the l ₂ norm, is the l ₀ semi-norm.

其中， $s_{β} (m) = m - \frac{1}{2 β} &dtri; {| | Gm - d | |}_{l_{2}}^{2},$ [·]_i表示矢量第i个元素。in, ${the s}_{β} (m) = m - \frac{1}{2 β} &dtri; {| | G m - d | |}_{l_{2}}^{2},$ [·] _i represents the i-th element of the vector.

可选的，对L0半范数模型构建模块，采取模型求解器件实现，其迭代Optionally, for the building block of the L0 semi-norm model, it is implemented by a model solving device, and its iteration

过程由同伦分析算法完成，过程如下：The process is completed by the homotopy analysis algorithm, and the process is as follows:

β_k,i+1＝β_k,i,i＝i+1；β _k,i+1 =β _k,i ,i=i+1;

其中，终止参数ε₀一般设定ε₀＝10^-1；Among them, the termination parameter ε ₀ is generally set to ε ₀ =10 ^-1 ;

其中，ε一般为ε＝10^-5；Among them, ε is generally ε=10 ^-5 ;

在另外一个实施例中，还提供了一种软件，该软件用于执行上述实施例及优选实施方式中描述的技术方案。In another embodiment, software is also provided, and the software is used to implement the technical solutions described in the above embodiments and preferred implementation manners.

在另外一个实施例中，还提供了一种存储介质，该存储介质中存储有上述软件，该存储介质包括但不限于：光盘、软盘、硬盘、可擦写存储器等。从以上的描述中，可以看出，本发明实施例实现了如下技术效果：一种基于L0半范数的绕射波成像方法及装置，通过限定模型稀疏性，使得求解模型更为符合实际地质情况，并且通过迭代逼近可以在一定程度上压制噪声，提高绕射波成像资料信噪比，该技术能够实现小断层、裂缝、溶洞等小尺度地质异常体精细刻画，在油气勘探储层研究中具有重要的应用价值。In another embodiment, there is also provided a storage medium, in which the software is stored, the storage medium includes but not limited to: optical discs, floppy disks, hard disks, rewritable memories, and the like. From the above description, it can be seen that the embodiment of the present invention achieves the following technical effects: a diffracted wave imaging method and device based on the L0 semi-norm, by limiting the sparsity of the model, the solution model is more in line with the actual geology and through iterative approximation, noise can be suppressed to a certain extent, and the signal-to-noise ratio of diffraction wave imaging data can be improved. This technology can achieve fine characterization of small-scale geological anomalies such as small faults, fractures, and caves. It has important application value.

最后需要说明的是：上述仅用以说明本发明而并非限制本发明所描述的技术方案；尽管本说明书对本发明已进行了详细的说明，但是，本领域的技术人员仍然可以对本发明进行修改或等同替换，一切不脱离本发明的精神和范围的技术方案及其改进，其均应涵盖在本发明的权利要求范围中。Finally, it should be noted that: the above is only used to illustrate the present invention rather than limit the technical solution described in the present invention; although this specification has described the present invention in detail, those skilled in the art can still modify or modify the present invention. Equivalent replacements, all technical solutions and their improvements that do not deviate from the spirit and scope of the present invention shall be included in the scope of the claims of the present invention.

Claims

1. a diffracted wave imaging method based on L0 half norm, is characterized in that, may further comprise the steps:

Step 101: Obtain a seismic data from which reflection waves have been removed as input data;

Step 102: Discretize the underground imaging space, select any imaging point as the current diffraction imaging point, and calculate the above-mentioned arbitrary diffraction imaging according to the relationship between the shot point and the receiver point in the input seismic data and the given velocity model Green's function of the point;

Step 103: cyclically execute other imaging points not selected in the underground imaging space in the above step 102, and calculate their corresponding Green's functions until the Green's functions of all imaging points in the underground imaging space are obtained;

Step 104: Construct a diffracted wave imaging model based on the L0 semi-norm according to the Green's function of all imaging points in the underground imaging space and the initial imaging model;

Step 105: The model is solved by the homotopy analysis iterative algorithm, and the diffraction wave imaging result is obtained.

2. a kind of diffracted wave imaging method based on L0 semi-norm according to claim 1, is characterized in that, the seismic data in the described step 101 is pre-stack shot seismic data or post-stack seismic data; reflection The wave removal method adopts the plane wave destruction filtering method.

3. a kind of diffracted wave imaging method based on L0 semi-norm according to claim 1, is characterized in that, the computing method of the Green's function in the described step 102,103, according to ray-tracing calculation, set out from shot point and pass through Travel time from the imaging point to the detector point, and store the amplitude weighted item.

4. a kind of diffraction wave imaging method based on L0 semi-norm according to claim 1, is characterized in that, in described step 104, constructs the diffraction wave imaging model based on L0 semi-norm, design is as follows:

min min {J J}_{α α} ((m m)) : : = = \frac{11}{22} {| | | | Gm G m - - d d | | | |}_{{l l}_{22}}^{22} + + α α {| | | | m m | | | |}_{{l l}_{00}}

Among them, min represents the minimum, J _α (m) is the objective function, m is the diffraction model to be solved, the mathematical symbol: = represents the definition as, G is the Green's function, d is the seismic data that removes the reflected wave, and α is the regularization factor, Indicates the l ₂ norm, is the l ₀ semi-norm.

By the neighboring point method, the above formula can be written as:

min min {H h}_{β β,, α α} (({m m}_{00},, m m)) : : = = (({G G}^{T T} G G {m m}_{00} - - {G G}^{T T} d d + + α α \frac{d d}{dm dm} {| | m m | | | |}_{{l l}_{00},, m m = = {m m}_{00}},, m m - - {m m}_{00})) + + \frac{β β}{22} {| | | | m m - - {m m}_{00} | | | |}_{{l l}_{22}}^{22} + + α α {| | | | m m | | | |}_{{l l}_{00}}

Among them, H _β,α (m ₀ ,m) is the approximation formula of J _α (m), m ₀ is the initial model, β is the adjustable regular parameter, (·,·) represents the inner product,

The minimum problem can be solved by the hard threshold operator:

{V V}_{β β} (({m m}_{00})) = = \underset{m m}{arg arg min min} {H h}_{β β,, α α} (({m m}_{00},, m m))

in,

{the s}_{β} (m) = m - \frac{1}{2 β} &dtri; {| | G m - d | |}_{l_{2}}^{2},

[·] _i represents _{the i} -th element of the vector.

5. a kind of diffracted wave imaging method based on L0 semi-norm according to claim 1, is characterized in that, in described step 105, the specific process of homotopy analysis iterative algorithm is as follows:

Step 1: Input Lipschitz parameter β ₀ , regularization parameter α ₀ , and require β ₀ ∈ [β _min , β _max ], β _min , β _max are the upper and lower limits of Puschitz parameter β ₀ respectively, initialization parameters, k= 0,ρ∈(0,1), initialize the model m ⁰ ;

Step 2: setting, i=0, m ^{k, 0} = m ^k , β _{k, 0} = β _k ;

Step 3: m ^k,i+1 =V _βk,i (m ^k,i+1 );

when

J_{α_{k}} (m^{k, i}) - J_{α_{k}} (m^{k, o = i + 1}) < \frac{η}{2} | | m^{k, i} - m^{k, i + 1} | |,

Execute β _k,i = min{γβ _k,i ,β _max }, m ^k,i+1 =T _βk,i (m ^k,i+1 );

β _k,i+1 =β _k,i ,i=i+1;

Return to step 3 until ||m ^k,i -m ^k,i+1 || _∞ ≤ε ₀ ;

Wherein, the termination parameter ε ₀ is ε ₀ =10 ^-1 ;

Step 4: m ^k+1 =m ^k,i ,β _k+1 =β _k,i ,α _k+1 =ρα _k ,k=k+1;

Return to step 2 until ||m ^k,i -m ^k,i+1 || _∞ ≤ε;

Wherein, ε is ε=10 ^-5 ;

Step 5: output the final iteration result, m ^* = m ^k .

6. A diffracted wave imaging device based on L0 half norm, is characterized in that, comprises the seismic data acquisition module (201) that is connected successively, is used for inputting a seismic data that removes reflected wave;

Diffraction imaging point selection module (202), used to select imaging points from the discretized underground imaging space as the Green's function calculation position;

Green's function calculation module (203), used for calculating travel time and amplitude compensation items from the shot point to the seismic receiver point through the diffraction imaging point;

The L0 semi-norm model building module (204), constructs a solution model based on the L0 semi-norm according to the Green's function of all imaging points in the underground space;

The model solving device (205) uses the homotopy analysis algorithm to solve the L0 model constructed by the Green's function and the diffraction point model, and obtains the diffraction wave imaging result.

7. A kind of diffraction wave imaging device based on L0 semi-norm according to claim 6, it is characterized in that, described Green's function calculation module (203), according to ray tracing calculation, departs from shot point and travels from imaging point to receiver point , and store the amplitude-weighted term.

8. a kind of diffracted wave imaging device based on L0 semi-norm according to claim 6, is characterized in that, described L0 semi-norm model construction module (204), is designed as follows:

min min {J J}_{α α} ((m m)) : : = = \frac{11}{22} {| | | | Gm G m - - d d | | | |}_{{l l}_{22}}^{22} + + α α {| | | | m m | | | |}_{{l l}_{00}}

Wherein, min represents minimization, J _α (m) is objective function, m is the diffraction model of solution, mathematical symbol: = represents definition as, G is Green's function, d is the seismic data that removes reflected wave, α is regularization factor, Indicates the l ₂ norm, is l ₀ semi-norm;

By the neighboring point method, the above formula can be written as:

min min {H h}_{β β,, α α} (({m m}_{00},, m m)) : : = = (({G G}^{T T} G G {m m}_{00} - - {G G}^{T T} d d + + α α \frac{d d}{dm dm} {| | m m | | | |}_{{l l}_{00},, m m = = {m m}_{00}},, m m - - {m m}_{00})) + + \frac{β β}{22} {| | | | m m - - {m m}_{00} | | | |}_{{l l}_{22}}^{22} + + α α {| | | | m m | | | |}_{{l l}_{00}}

The minimum problem is solved by the hard threshold operator:

{V V}_{β β} (({m m}_{00})) = = \underset{m m}{arg arg min min} {H h}_{β β,, α α} (({m m}_{00},, m m))

in,

{the s}_{β} (m) = m - \frac{1}{2 β} &dtri; {| | G m - d | |}_{l_{2}}^{2},

[·] _i represents the i-th element of the vector.

9. according to claim 6 a kind of diffracted wave imaging device based on L0 semi-norm, it is characterized in that, described model solving device (205) iterative process is finished by homotopy analysis algorithm, and process is as follows: