CN106405648B - The imaging method and device of diffracted wave - Google Patents

Abstract

The present invention provides a diffracted wave imaging method and device, which relate to the technical field of seismic wave imaging, including: acquiring initial shot data, wherein the initial shot data carries geological information in the target area; Data preprocessing is performed on the shot set data to obtain the common offset diffraction wave data, wherein the common offset diffraction wave data have the same offset; based on the offset of the diffraction wave in the common offset diffraction wave data The re-weighted imaging model is constructed using the moving velocity and common offset diffraction wave data; the re-weighted imaging model is calculated using a preset algorithm, and the calculation result is used as the target imaging result of the diffracted wave, which solves the problem of using in the prior art In the process of determining faults and collapsed column areas by diffraction wave imaging technology, poor imaging results lead to technical problems of inaccurate determination.

Description

Imaging method and device for diffracted waves

technical field

The invention relates to the technical field of seismic wave imaging, in particular to a diffraction wave imaging method and device.

Background technique

In the process of coal mining, in order to prevent geological disasters, it is necessary to detect small-scale discontinuous and heterogeneous geological structures before mining, such as faults and collapse columns. Among them, faults and collapse columns are efficient and safe for mechanized coal mines. Production is critical. Faults can cause coal and rock formations to connect with strong aquifers, and then induce water seepage accidents and even well flooding; if a large amount of gas accumulates in the fault fracture zone, it will cause gas outburst accidents, which will cause irreparable damage to people. In areas where collapse pillars are developed, coal seams and surrounding rocks in coal-measure strata are often destroyed, resulting in a reduction in coal reserves and even loss of mining value in this area; moreover, it is difficult to arrange longwall mining face in this area, This has also seriously hindered mechanized coal mining.

In order to reduce property losses, relevant scholars in the industry have conducted a lot of research on the identification of faults and collapsed columns. At present, the methods used in the industry to identify faults and collapsed columns are mainly based on the theory of seismic reflection waves, and on this basis, the analysis of seismic coherent volumes is carried out; spectral decomposition algorithms can also be used. However, the reflection wave theory is based on the assumption that the underground reflection interface is smooth and infinite, and due to the limited resolution, the reflection wave theory is difficult to meet the detection of small-scale discontinuous and heterogeneous geological bodies.

In addition to the theory based on reflected waves, faults and collapse columns can also be determined using seismic diffraction waves. In the process of using diffraction waves to identify faults and collapse columns, diffraction wave separation and diffraction wave imaging are two major issues. The key issue. Relevant scholars have made many attempts on the separation of diffracted waves, for example, the method of diffracted wave separation based on the combination of local dip filtering and prediction and inversion, the method of diffracted wave field separation based on plane wave destruction filter (PWD), etc. .

At present, most of the core content of diffraction wave technology focuses on the separation of diffraction waves, and there is no research on the difference between the physical inversion models of diffraction points and reflection interfaces. Although the sparse discontinuity of the spatial distribution of diffraction information is considered in the prior art, the constraint on the solution model is too strong in the prior art, and it is difficult to achieve the optimal imaging effect during the inversion process.

Contents of the invention

The purpose of the present invention is to provide a diffracted wave imaging method and device to alleviate the technical problem of inaccurate determination due to poor imaging effect in the process of using diffractive wave imaging technology to determine faults and collapsed column areas in the prior art .

According to an aspect of an embodiment of the present invention, a diffracted wave imaging method is provided, including: acquiring initial shot data, wherein the initial shot data carries geological information in the target area, and the geological information includes At least one of the following: geological information of stratum layer structure, geological information of fault form, geological information of karst caves; performing data preprocessing on the acquired initial shot data to obtain common offset diffraction wave data, Wherein, the common offset diffraction wave data has the same offset; based on the offset velocity of the diffraction wave in the common offset diffraction wave data and the common offset diffraction wave data to construct A re-weighted imaging model: using a preset algorithm to calculate the re-weighted imaging model, and use the calculation result as the target imaging result of the diffracted wave.

Further, constructing a reweighted imaging model based on the migration velocity of the diffraction wave in the common offset diffraction wave data and the common offset diffraction wave data includes: according to the migration velocity of the diffraction wave Calculate the target Green's function, wherein the target Green's function represents the propagation time and amplitude compensation factor of the diffracted wave from the shot point position to the receiver point position through any imaging point position in the underground imaging space; based on the target Green's function and constructing the reweighted imaging model with the co-offset diffraction wave data.

Further, constructing the reweighted imaging model based on the target Green's function and the common offset diffraction wave data includes: using the formula Constructing the re-weighted imaging model, w _i is the re-weighting coefficient, G is the matrix form of the target Green's function, _ri is the diffraction wave imaging result r( _xi ) of the imaging point x _i in the underground imaging space The scalar form of , d _obs is the common offset diffraction wave data, i takes the value from 1 to N in turn, and N represents the number of imaging points in the underground imaging space.

Further, the preset algorithm includes an adaptive homotopy algorithm, using the preset algorithm to calculate the reweighted imaging model, and using the calculation result as the target imaging result of the diffracted wave includes: The adaptive homotopy algorithm performs a superposition operation on the reweighted imaging model, and uses the superposition result as the target imaging result.

Further, performing a superposition operation on the reweighted imaging model by using an adaptive homotopy algorithm, and using the superimposed result as the target imaging result includes: taking the preset initial parameter value of the target parameter as the current parameter value, perform the following steps until the parameter value of the target parameter satisfies the preset condition, wherein the target parameter includes: the reweighting coefficient, the diffraction wave imaging result of the imaging point x _i in the underground imaging space, Iteration termination parameter; the first calculation step, according to the formula Computes the scalar value of the parameter value for the current reweighting factor and computes the current update direction vector Wherein, G _Γ is the matrix that is made up of the target column vector in the matrix G of described target Green's function, and the sequence number of described target column vector in the matrix of described Green's function is corresponding to the index number in the current set Γ, The index number in the current set Γ is composed of the serial number corresponding to the non-zero value in the current inversion solution vector r( _xi ), and the current inversion solution vector r( _xi ) is composed of r _i , a diagonal matrix W and diagonal The diagonal elements of are respectively composed of w _i and composition; the second calculation step, according to the formula Calculate the current update step size Δr _i , where s _i is the ith element of the current update direction vector s; the third calculation step, according to the formula r _i := r _i +(Δr _i )s _i and the formula Calculate the current iteration result; the first update step is used to delete the element corresponding to i in the current set Γ when it is judged that Δr _i < 1, or, when it is judged that Δr _i ≥ 1, in Add a new index number in the current set Γ; the second update step, according to the formula Updating the parameter value of the current reweighting coefficient; judging step, judging whether the updated current reweighting coefficient and the current iteration termination parameter satisfy the preset condition, wherein the preset condition is max(w _i ) ≤τ holds true, or the current iteration termination parameter is greater than or equal to the target threshold, i=1, 2,...,N; wherein, if it is judged that the preset condition is satisfied, the current inversion solution vector r is output ( _xi ), if it is judged that the preset condition is not met, then control the parameter value of the current iteration termination parameter to increase the preset value, and add the parameter value of _ri after iteration in the third calculation step and The updated parameter value of the current reweighting coefficient in the second updating step is used as the current parameter value, and the execution of the first calculation step is returned.

Further, performing data preprocessing on the acquired initial shot data to obtain the co-offset diffraction wave data includes: filtering the initial shot data to obtain the co-offset shot data, wherein the The common offset shot data has the same offset distance; according to the sparse Radon hyperbolic transformation method, the common offset shot data is transformed to obtain the Radon domain after transformation; the distance between the reflected wave and the Radon domain is cut The part corresponding to the frequency spectrum; the inverse sparse Radon hyperbolic transformation is performed on the cut-off Radon domain to obtain the common offset diffraction wave data.

According to an aspect of an embodiment of the present invention, there is also provided an imaging device for diffracted waves, including: an acquisition unit, configured to acquire initial shot data, wherein the initial shot data carries geological information in the target area , the geological information includes at least one of the following: geological information of stratum layer structure, geological information of fault morphology, geological information of karst caves; a processing unit, configured to perform data preprocessing on the acquired initial shot data , to obtain the common offset diffraction wave data, wherein the common offset diffraction wave data has the same offset; the construction unit is used for diffracting waves based on the common offset diffraction wave data The migration velocity and the common offset diffraction wave data construct a reweighted imaging model; the calculation unit is used to calculate the reweighted imaging model using a preset algorithm, and use the calculation result as the diffraction wave target imaging results.

Further, the construction unit includes: a first calculation subunit, which is used to calculate the target Green's function according to the migration velocity of the diffracted wave, wherein the target Green's function represents that the diffracted wave passes from the shot point position to The propagation time and amplitude compensation factor from any imaging point position in the underground imaging space to the receiver point position; constructing a subunit for constructing the reweighted imaging model based on the Green's function and the common offset diffraction wave data .

Further, the construction subunit includes: a construction module, which is used to pass the formula Constructing the re-weighted imaging model, _wi is the re-weighting coefficient, G is the matrix form of the target Green's function, ri is the scalar form of the diffraction wave imaging result of the imaging point _xi in the underground imaging space, d _obs is the common-offset diffraction wave data, i is sequentially taken from 1 to N, and N represents the number of imaging points in the underground imaging space.

Further, the preset algorithm includes an adaptive homotopy algorithm, and the calculation unit includes: a second calculation subunit, which is used to perform superposition operations on the reweighted imaging model by using the adaptive homotopy algorithm to obtain the superposition The result after adding is used as the target imaging result.

In the diffraction wave imaging method provided in the embodiment of the present invention, the initial shot data carrying geological information is first obtained, and then the acquired data is subjected to data preprocessing to obtain the common offset diffraction wave data, and then Next, a reweighted imaging model is constructed according to the common offset diffraction wave data and migration velocity obtained after processing. Finally, the preset algorithm is used to calculate the reweighted imaging model to obtain the target imaging result of the diffraction wave. In the embodiment of the present invention, the method of determining the diffraction wave imaging results through the re-weighted imaging model achieves the purpose of accurately detecting faults and collapse columns, and alleviates the problem of using diffraction wave imaging technology to determine faults and collapses in the prior art. In the process of the column area, poor imaging results lead to inaccurate technical problems, thus achieving the technical effect of improving the detection accuracy of faults and collapsed columns.

Description of drawings

In order to more clearly illustrate the specific implementation of the present invention or the technical solutions in the prior art, the following will briefly introduce the accompanying drawings that need to be used in the specific implementation or description of the prior art. Obviously, the accompanying drawings in the following description The drawings show some implementations of the present invention, and those skilled in the art can obtain other drawings based on these drawings without any creative work.

FIG. 1 is a flow chart of a diffracted wave imaging method according to an embodiment of the present invention;

2 is a flowchart of a method for constructing a reweighted imaging model according to an embodiment of the present invention;

FIG. 3 is a flow chart of calculating a reweighted imaging model by an adaptive homotopy algorithm according to an embodiment of the present invention;

4 is a flow chart of a method for processing initial shot data according to an embodiment of the present invention;

Fig. 5 is a schematic diagram of an imaging device for diffracted waves according to an embodiment of the present invention.

Detailed ways

The technical solutions of the present invention will be clearly and completely described below in conjunction with the accompanying drawings. Apparently, the described embodiments are some of the embodiments of the present invention, but not all of them. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without making creative efforts belong to the protection scope of the present invention.

In the description of the present invention, it should be noted that the terms "center", "upper", "lower", "left", "right", "vertical", "horizontal", "inner", "outer" etc. The indicated orientation or positional relationship is based on the orientation or positional relationship shown in the drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the referred device or element must have a specific orientation, or in a specific orientation. construction and operation, therefore, should not be construed as limiting the invention. In addition, the terms "first", "second", and "third" are used for descriptive purposes only, and should not be construed as indicating or implying relative importance.

In the description of the present invention, it should be noted that unless otherwise specified and limited, the terms "installation", "connection" and "connection" should be understood in a broad sense, for example, it can be a fixed connection or a detachable connection. Connected, or integrally connected; it can be mechanically connected or electrically connected; it can be directly connected or indirectly connected through an intermediary, and it can be the internal communication of two components. Those of ordinary skill in the art can understand the specific meanings of the above terms in the present invention in specific situations.

Fig. 1 is a flow chart of a method for imaging a diffracted wave according to an embodiment of the present invention. As shown in Fig. 1, the method includes the following steps:

Step S102, acquiring initial shot data, wherein the initial shot data carries geological information in the target area, and the geological information includes at least one of the following: geological information of stratum layer structure, geological information of fault morphology, geological information of karst caves information.

In the embodiment of the present invention, the shot data can also become seismic data, which is the seismic wave data detected by the geophone at the detection point, wherein the seismic wave data includes reflected wave data and diffracted wave data, except for reflected waves and diffracted waves In addition, the seismic wave also includes other waveforms, but in the embodiment of the present invention, the reflected wave and the diffracted wave are mainly processed, and then the target imaging result is obtained. Therefore, in the embodiment of the present invention, in addition to the reflected wave Waveforms other than waves and diffracted waves are not described in detail.

Assuming that relevant technicians set up a shot point in the target area, when the shot point explodes, seismic waves will be generated. At this time, a plurality of geophones may be set on the ground in the target area, that is, a plurality of geophone points may be set, and then seismic waves at each geophone point may be detected by the multiple geophones. It should be noted that the data described above can also become single-shot data, and multiple single-shot data constitute shot set data.

Step S104, performing data preprocessing on the acquired initial shot data to obtain common-offset diffraction wave data, wherein the common-offset diffraction wave data have the same offset.

In the embodiment of the present invention, after the shot set data is acquired, data preprocessing needs to be performed on the shot set data, and the common offset diffraction wave data is obtained after the processing. In the embodiment of the present invention, the co-offset diffracted wave data obtained after processing is data with the same offset, where the offset is the horizontal distance between the shot point position and the receiver point position. That is to say, in the common-offset diffraction wave data, the distance between the shot point and the receiver point is equal.

It should be noted that the initial shot data includes diffraction wave data and reflected wave data, and the process of processing the initial shot data includes the process of extracting diffraction wave data from the initial shot data, specifically, extracting The process will be described in detail in the following examples.

Step S106, constructing a reweighted imaging model based on the migration velocity of the diffraction wave in the common offset diffraction wave data and the common offset diffraction wave data.

In the embodiment of the present invention, after obtaining the common offset diffraction wave data in step S104, the migration velocity file can be loaded to obtain the migration velocity of the diffraction wave stored in the migration velocity file; furthermore, according to the offset The reweighted imaging model was constructed using the moving velocity and co-offset diffraction wave data. It should be noted that, in the embodiment of the present invention, the above-mentioned reweighted imaging model is preferably a Kirchhoff high-resolution imaging model, and reweighting is used as a constraint value of the Kirchhoff high-resolution imaging model.

It should be noted that, in the embodiment of the present invention, the migration velocity file is a file obtained in advance by the relevant technical personnel, which includes the propagation velocity of seismic waves (for example, diffracted waves and reflected waves) in the underground imaging space . Specifically, relevant technicians can collect relevant data of migration velocity in the field, and then load the relevant data of migration velocity collected through the observation system, and then perform processing processes such as denoising processing and migration velocity analysis on the data, Finally, the processed data is used as an offset velocity file.

Step S108, using a preset algorithm to calculate the reweighted imaging model, and use the calculation result as the target imaging result of the diffracted wave.

In the embodiment of the present invention, after the reweighted imaging model is constructed in step S106, the reweighted imaging model can be calculated according to a preset algorithm, and the result obtained after the calculation is used as the imaging result of the diffracted wave (that is, the target imaging result ).

In the diffraction wave imaging method provided in the embodiment of the present invention, the initial shot data carrying geological information is first obtained, and then the acquired data is subjected to data preprocessing to obtain the common offset diffraction wave data, and then Next, a reweighted imaging model is constructed according to the common offset diffraction wave data and migration velocity obtained after processing. Finally, the preset algorithm is used to calculate the reweighted imaging model to obtain the target imaging result of the diffraction wave. In the embodiment of the present invention, the method of determining the diffraction wave imaging results through the re-weighted imaging model achieves the purpose of accurately detecting faults and collapse columns, and alleviates the problem of using diffraction wave imaging technology to determine faults and collapses in the prior art. In the process of the column area, poor imaging results lead to inaccurate technical problems, thus achieving the technical effect of improving the detection accuracy of faults and collapsed columns.

Fig. 2 is a flow chart of a method for constructing a reweighted imaging model according to an embodiment of the present invention, as shown in Fig. 2, based on the migration velocity and the common offset distance The construction of a reweighted imaging model from diffraction wave data includes the following steps:

Step S201, calculate the target Green's function according to the migration velocity of the diffracted wave, wherein the target Green's function represents the propagation time and amplitude compensation factor of the diffracted wave from the shot point position to the receiver point position via any imaging point position in the underground imaging space ;

Step S202, constructing a reweighted imaging model based on the target Green's function and the co-offset diffraction wave data.

In the embodiment of the present invention, after the common-offset diffraction wave data is acquired, a reweighted imaging model can be constructed according to the common-offset diffraction wave data and the pre-acquired migration velocity file. When constructing the reweighted imaging model, on the premise that the migration velocity of the diffraction wave is known, the diffraction wave can be calculated according to the migration velocity of the diffraction wave. The diffraction wave starts from the shot point and passes through any imaging point in the underground imaging space. The travel time (that is, the propagation time) to each geophone location, and then, the target Green's function is established according to the calculated travel time.

After the target Green's function is calculated, a reweighted imaging model can be constructed based on the calculated target Green's function and the co-offset diffraction wave data. Preferably, in the embodiment of the present invention, a reweighted imaging model may be constructed using a Kirchhoff imaging algorithm. Among them, the reweighted imaging model obtained by using the Kirchhoff imaging algorithm is a Kirchhoff high-resolution imaging model based on reweighted sparse constraints.

In an optional embodiment of the present invention, a reweighted imaging model is constructed based on the target Green's function and the common offset diffraction wave data, specifically:

by formula Construct a reweighted imaging model, where w _i is the reweighting coefficient, G is the matrix form of the target Green's function, _{ri is the scalar form of the diffraction wave imaging result of the imaging point x i in} the underground imaging space, and d _obs is the common bias For distance-shifted diffraction wave data, i takes the value from 1 to N in sequence, and N represents the number of imaging points in the underground imaging space.

In the embodiment of the present invention, the reweighted imaging model can be constructed according to the above formula. In the above formula, w _i is the reweighted coefficient, wherein, w _i >0; G is the matrix form of the target Green's function; _ri is the underground imaging Diffraction wave imaging result of imaging point _xi in space, r( _xi ) is the vector representation of _ri , where i=1,2,...,N, N represents the number of discrete sample points in the underground imaging space, That is, the number of imaging points; the vector d _obs is the above-mentioned common offset diffraction wave data.

using the formula After the reweighted imaging model is constructed, the model needs to be calculated, and then the calculation result is obtained, wherein the calculation result is used to determine the target imaging result of the diffracted wave. It should be noted that, in the above formula middle, for the constraint part.

Preferably, in the embodiment of the present invention, the preset algorithm may be selected as an adaptive homotopy algorithm (also called a homotopy adaptive inversion algorithm). Specifically, an adaptive homotopy algorithm may be used to perform a superposition operation on the reweighted imaging model, and the superposition result may be used as a target imaging result. Among them, the homotopy adaptive inversion algorithm can invert the oblique formation model with a relatively complex structure, and can converge at a relatively fast speed. The specific calculation process of the adaptive homotopy algorithm will be described in detail below.

Fig. 3 is a flow chart of calculating a reweighted imaging model according to an adaptive homotopy algorithm according to an embodiment of the present invention. As shown in Fig. 3, the reweighted imaging model is superposed using the adaptive homotopy algorithm to obtain The subsequent result is used as the target imaging result, including the following steps S301 to S307:

Before executing the calculation method described in the following steps S301 to S307, first, obtain the initial parameter value of the preset target parameter, and use the obtained initial parameter value as the current parameter value to perform the following steps S301 to S307, Until the parameter value of the target parameter satisfies the preset conditions in the following description, wherein the target parameter includes: the reweighting coefficient w _i , the diffraction wave imaging result r _i of the imaging point x _i in the underground imaging space, and the iteration termination parameter τ.

In the embodiment of the present invention, the initial parameter values of the above-mentioned target parameters can be selected in the following manner: Among them, max means to obtain the maximum value, the scalar τ is the iteration termination parameter, τ is determined by the user, and the vector g _i is the i-th column of the above-mentioned target Green's function matrix G.

S301, the first calculation step, according to the formula Computes the scalar value of the parameter value for the current reweighting factor and computes the current update direction vector Among them, G _Γ is a matrix composed of target column vectors in the matrix G of the target Green's function, the serial number of the target column vector in the matrix of Green's functions corresponds to the index number in the current set Γ, and the index in the current set Γ The number is composed of the serial number corresponding to the non-zero value in the current inversion solution vector r( _xi ), the current inversion solution vector r( _xi ) is composed of r _i , the diagonal matrix W and the diagonal matrix The diagonal elements of are respectively composed of w _i and composition;

In the embodiment of the present invention, after the system obtains the current parameter value of the target parameter, it can follow the formula Calculate the scalar value of the current reweighting coefficient, and then calculate the current update direction vector S according to the matrix of the scalar value and the matrix of the current reweighting coefficient, wherein the current update direction vector is used to determine the direction of change of r _i and w _i , that is to say, it is used to determine the change direction of the current update step Δr _i in the second calculation step.

S302, the second calculation step, according to the formula Calculate the current update step size Δr _i , where s _i is the ith element of the current update direction vector s;

After the current update direction vector is calculated in the first calculation step, it can be calculated according to the formula Calculate the current update step size Δr _i .

S303, the third calculation step, according to the formula r _i :=r _i +(Δr _i )s _i and the formula Calculate the result of the current iteration;

After the current update step size Δr _i is calculated in the second calculation step, the calculation result r _i after iteration can be calculated according to the formula _ri :=ri+(Δr _i ) _{s i ;} and according to the formula Calculate the calculation result w _i after iteration.

S304, the first update step is used to delete the element corresponding to i in the current set Γ when it is judged that Δr _i <1, or, if it is judged that Δr _i ≥ 1, add new index number;

In the embodiment of the present invention, after the superposition r _i and w _i are calculated in the above third calculation step, the elements in the set Γ need to be updated for subsequent superposition. When updating the elements in the set Γ, first judge whether Δr _i <1 is true, and if it is judged that Δr _i <1 is true, then remove the element corresponding to i in the set Γ, if it is judged that Δr _i <1 is not true, then Add a new element in the set Γ (wherein, the new element is the above-mentioned new index number), and the selection method of the added element is: Among them, Γ ^c is composed of the serial number corresponding to the zero value in the inversion solution vector r( _xi ), and argmax means to obtain the maximum value.

S305, the second updating step, according to the formula Update the parameter value of the current reweighting coefficient;

After updating the elements in the set Γ according to the above-mentioned first update step, it is also necessary to update the parameter value of the current reweighting coefficient, and the specific update can be according to the formula

S306, judging step, judging whether the updated current reweighting coefficient and the current iteration termination parameter meet the preset condition, wherein the preset condition is that max(w _i )≤τ holds true, or the current iteration termination parameter is greater than or equal to the target threshold , i=1,2,...,N; wherein, if it is judged that the preset condition is satisfied, then execute step S307 to output the current inversion solution vector r( _xi ), if it is judged that the preset condition is not satisfied, control the current The parameter value of the iteration termination parameter is increased by a preset value, and the parameter value of r _i after iteration in the third calculation step and the parameter value of the current reweighting coefficient after updating in the second update step are used as the current parameter value, and return to execute the first A calculation step.

After the parameter value of the current reweighting coefficient is updated in the above-mentioned second update step, it is judged whether max(w _i )≤τ holds true, or whether the current iteration termination parameter is greater than or equal to the target threshold, wherein the target threshold represents the above steps S301 to The maximum number of iterations required by step S306. Among them, if it is judged that max(w _i )≤τ holds true, or it is judged that the current iteration termination parameter is greater than or equal to the target threshold, then stop and output the vector r( _xi ); otherwise, the current iteration termination parameter increases the preset value (for example, increase by 1), then, the parameter value of r _i after iteration in the third calculation step and the parameter value of the current reweighting coefficient after updating in the second update step are used as the current parameter value, return to continue to execute the first calculation Steps until the calculated result satisfies the preset condition.

Finally, after the output vector r( _xi ) is obtained in the judgment step, the output vector is substituted into the above formula , and take the calculation result of the above formula as the target imaging result.

Fig. 4 is a flow chart of a method for processing initial shot data according to an embodiment of the present invention. As shown in Fig. 4, data preprocessing is performed on the acquired initial shot data to obtain common offset diffraction wave data Including the following steps:

Step S401, screening the initial shot data to obtain co-offset shot data, wherein the co-offset shot data have the same offset distance;

Step S402, transforming the co-offset shot data according to the sparse Radon hyperbolic transformation method to obtain the transformed Radon field;

Step S403, cutting off the part corresponding to the frequency spectrum of the reflected wave in the Radon domain;

Step S404, perform inverse-sparse Radon hyperbolic transformation on the cut-out Radon field to obtain common-offset diffraction wave data.

In the embodiment of the present invention, since the obtained initial shot data contains multiple types of data, in order to obtain the diffraction wave seismic data used in the embodiment of the present invention, it is necessary to analyze the obtained initial shot data Data are preprocessed accordingly. First, the acquired seismic initial shot data is loaded in the observation system; then, denoising and other processing are performed on the initial shot data. After processing in the above manner, the seismic initial shot data is screened according to the keywords (for example, offset and track number) in the processed seismic initial shot file to obtain the common offset seismic data, That is, the shot data with the same offset distance is screened, where the offset distance is the horizontal distance between the shot point position and the receiver point position on the acquisition ground, and the trace number is the number of the geophone in the seismic shot data; the screening operation is based on Corresponding seismic data are extracted with the same offset.

After screening the common-offset shot data, according to the sparse Radon hyperbolic transformation method, the common-offset shot data is separated into the common-offset seismic data containing only diffracted waves. The specific separation methods include: The Radon hyperbolic transform transforms the data of each common-offset shot into the Radon domain, and then cuts off the spectral component corresponding to the reflected wave in the Radon domain, and finally performs an inverse-sparse Radon hyperbolic transform on the cut-off Radon domain spectrum to obtain Separated co-offset diffracted wave data.

To sum up, the diffraction wave imaging method provided by the embodiment of the present invention includes: first sorting the preprocessed common-offset seismic data; then, according to the above-mentioned sorted common-offset seismic data, according to Sparse Radon hyperbolic transformation method to separate the common-offset seismic data containing only diffracted waves; next, according to the input migration velocity file and the separated common-offset diffracted wave data, the Kirchhoff imaging method is used to construct A high-resolution imaging model based on reweighted sparse constraints; finally, an adaptive homotopy algorithm is used to solve the high-resolution imaging model based on reweighted sparse constraints, and the diffraction wave imaging result is obtained. On the basis of in-depth analysis of the limitations of the conventional diffraction physical model, the present invention proposes a diffraction wave adaptive sparse imaging method based on a reweighted model. Compared with the conventional diffraction wave imaging method, the model can be self-adaptive Adjust the weighting coefficient, that is, increase the weight coefficient of the position where the model solution value is small and reduce the weight coefficient of the position where the model solution value is large, so as to ensure the stability and convergence of the iterative inversion process of diffraction wave imaging, and achieve optimal diffraction wave inversion Therefore, it can accurately detect fractures and small-scale collapse columns, and reduce safety hazards such as induced water inrush and gas leakage in coal mining.

The embodiment of the present invention also provides a diffracted wave imaging device, the diffracted wave imaging device is mainly used to implement the diffracted wave imaging method provided in the above content of the embodiment of the present invention, the following provides the embodiment of the present invention The imaging device of the diffracted wave will be introduced in detail.

Fig. 5 is a schematic diagram of a diffracted wave imaging device according to an embodiment of the present invention. As shown in Fig. 5, the diffracted wave imaging device mainly includes an acquisition unit 51, a processing unit 53, a construction unit 55 and a computing unit 57 ,in:

The acquiring unit 51 is configured to acquire initial shot data, wherein the initial shot data carries geological information in the target area, and the geological information includes at least one of the following: geological information of stratum layer structure, geological information of fault form, karst Geological information of caves;

In the embodiment of the present invention, the shot data can also become seismic data, which is the seismic wave data detected by the geophone at the detection point, wherein the seismic wave data includes reflected wave data and diffracted wave data, except for reflected waves and diffracted waves In addition, the seismic wave also includes other waveforms, but in the embodiment of the present invention, the reflected wave and the diffracted wave are mainly processed, and then the target imaging result is obtained. Therefore, in the embodiment of the present invention, in addition to the reflected wave Waveforms other than waves and diffracted waves are not described in detail.

Assuming that relevant technicians set up a shot point in the target area, when the shot point explodes, seismic waves will be generated. At this time, a plurality of geophones may be set on the ground in the target area, that is, a plurality of geophone points may be set, and then seismic waves at each geophone point may be detected by the multiple geophones. It should be noted that the data described above can also become single-shot data, and multiple single-shot data constitute shot set data.

The processing unit 53 is configured to perform data preprocessing on the acquired initial shot data to obtain common-offset diffraction wave data, wherein the common-offset diffraction wave data have the same offset;

In the embodiment of the present invention, after the shot set data is acquired, data preprocessing needs to be performed on the shot set data, and the common offset diffraction wave data is obtained after the processing. In the embodiment of the present invention, the co-offset diffracted wave data obtained after processing is data with the same offset, where the offset is the horizontal distance between the shot point position and the receiver point position. That is to say, in the common-offset diffraction wave data, the distance between the shot point and the receiver point is equal.

It should be noted that the initial shot data includes diffraction wave data and reflected wave data. The process of processing the initial shot data includes the process of extracting diffraction wave data from the initial shot data. Specifically, the extraction process A detailed introduction will be made in the following examples.

A construction unit 55, configured to construct a reweighted imaging model based on the migration velocity of the diffraction wave in the common offset diffraction wave data and the common offset diffraction wave data;

In the embodiment of the present invention, after obtaining the common offset diffraction wave data in step S104, the migration velocity file can be loaded to obtain the migration velocity of the diffraction wave stored in the migration velocity file; furthermore, according to the offset The reweighted imaging model was constructed using the moving velocity and co-offset diffraction wave data. It should be noted that, in the embodiment of the present invention, the above-mentioned reweighted imaging model is preferably a Kirchhoff high-resolution imaging model, and reweighting is used as a constraint value of the Kirchhoff high-resolution imaging model.

It should be noted that, in the embodiment of the present invention, the migration velocity file is a file obtained in advance by the relevant technical personnel, which includes the propagation velocity of seismic waves (for example, diffracted waves and reflected waves) in the underground imaging space . Specifically, relevant technicians can collect relevant data of migration velocity in the field, and then load the relevant data of migration velocity collected through the observation system, and then perform processing processes such as denoising processing and migration velocity analysis on the data, Finally, the processed data is used as an offset velocity file.

The calculation unit 57 is configured to use a preset algorithm to calculate the reweighted imaging model, and use the calculation result as the target imaging result of the diffracted wave.

In the embodiment of the present invention, after the reweighted imaging model is constructed in step S106, the reweighted imaging model can be calculated according to a preset algorithm, and the result obtained after the calculation is used as the imaging result of the diffracted wave (that is, the target imaging result ).

In the diffraction wave imaging method provided in the embodiment of the present invention, the initial shot data carrying geological information is first obtained, and then the acquired data is subjected to data preprocessing to obtain the common offset diffraction wave data, and then Next, a reweighted imaging model is constructed according to the common offset diffraction wave data and migration velocity obtained after processing. Finally, the preset algorithm is used to calculate the reweighted imaging model to obtain the target imaging result of the diffraction wave. In the embodiment of the present invention, the method of determining the diffraction wave imaging results through the re-weighted imaging model achieves the purpose of accurately detecting faults and collapse columns, and alleviates the problem of using diffraction wave imaging technology to determine faults and collapses in the prior art. In the process of the column area, poor imaging results lead to inaccurate technical problems, thus achieving the technical effect of improving the detection accuracy of faults and collapsed columns.

Optionally, the construction unit includes: a first calculation subunit, which is used to calculate the target Green's function according to the migration velocity of the diffracted wave, wherein the target Green's function indicates that the diffracted wave is imaged by any one of the underground imaging spaces from the shot point position The travel time and amplitude compensation factor from the point position to the receiver point position; construct a subunit for constructing a reweighted imaging model based on Green's function and common offset diffraction wave data.

Optionally, the building subunits include: building blocks for passing the formula Construct the reweighted imaging model, w _i is the reweighting coefficient, G is the matrix form of the target Green's function, ri is the scalar form of the diffraction wave imaging result of the imaging point _xi in the underground imaging space, d _obs is the common offset around For radio wave data, i takes the value from 1 to N in sequence, and N represents the number of imaging points in the underground imaging space.

Optionally, the preset algorithm includes an adaptive homotopy algorithm, and the computing unit includes: a second computing subunit, configured to perform superposition operations on the reweighted imaging model by using the self-adaptive homotopy algorithm, and obtain a result after superposition as Target imaging results.

Optionally, the second calculation subunit includes: taking the preset initial parameter value of the target parameter as the current parameter value, and performing the following steps until the parameter value of the target parameter satisfies a preset condition, wherein the target parameter includes: a reweighting coefficient , the diffraction wave imaging result of the imaging point xi in the underground imaging space, the iteration termination parameter; the first calculation module is used to follow the formula Computes the scalar value of the parameter value for the current reweighting factor and computes the current update direction vector Among them, G _Γ is a matrix composed of target column vectors in the matrix G of the target Green's function, the serial number of the target column vector in the matrix of Green's functions corresponds to the index number in the current set Γ, and the index in the current set Γ The number is composed of the sequence number corresponding to the non-zero value in the inversion solution vector r( _xi ), the inversion solution vector r( _xi ) is composed of r _i , the diagonal matrix W and the diagonal matrix The diagonal elements of are respectively composed of w _i and Composition; the second calculation module is used to follow the formula Calculate the current update step size Δr _i , where s _i is the ith element of the current update direction vector s; the third calculation module is used to follow the formula r _i := r _i +(Δr _i )s _i and the formula Calculate the current iteration result; the first update module is used to delete the element corresponding to i in the current set Γ when it is judged that Δr _i < 1 is true, or when it is judged that Δr _i ≥ 1 is true, in the current Add a new index number in the set Γ; the second update module is used to follow the formula Updating the parameter value of the current reweighting coefficient; a judging module for judging whether the updated current reweighting coefficient and the current iteration termination parameter meet a preset condition, wherein the preset condition is that max(w _i )≤τ holds true, or, The current iteration termination parameter is greater than or equal to the target threshold, i=1,2,...,N; where, if it is judged that the preset condition is satisfied, the current inversion solution vector r( _xi ) is output, and if it is judged that the preset condition is not satisfied condition, then control the parameter value of the current iteration termination parameter to increase the preset value, and use the parameter value of r _i after iteration in the third calculation step and the parameter value of the current reweighting coefficient after updating in the second update step as the current parameter value, returns to perform the first calculation step.

Optionally, the processing unit includes: a screening module, configured to filter the initial shot data to obtain co-offset shot data, wherein the co-offset shot data have the same offset; the first transformation module uses According to the sparse Radon hyperbolic transformation device, the common offset shot data is transformed to obtain the transformed Radon domain; the cutting module is used to cut the part corresponding to the frequency spectrum of the reflected wave in the Radon domain; the second transformation module uses The inverse-sparse Radon hyperbolic transformation is performed on the cut-off Radon domain to obtain the common-offset diffraction wave data.

Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present invention, rather than limiting them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that: It is still possible to modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements for some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the technical solutions of the various embodiments of the present invention. scope.

Claims (3)

1. a kind of imaging method of diffracted wave, it is characterized in that, comprising: Acquiring initial shot data, wherein the initial shot data carries geological information in the target area, and the geological information includes at least one of the following: geological information of stratum layer structure, geological information of fault form, karst cave geological information; performing data preprocessing on the acquired initial shot data to obtain common-offset diffraction wave data, wherein the common-offset diffraction wave data have the same offset; Constructing a reweighted imaging model based on the migration velocity of the diffraction wave in the common offset diffraction wave data and the common offset diffraction wave data; Using a preset algorithm to calculate the reweighted imaging model, and use the calculation result as the target imaging result of the diffracted wave; Constructing a reweighted imaging model based on the migration velocity of the diffraction wave in the common offset diffraction wave data and the common offset diffraction wave data includes: Calculate the target Green's function according to the migration velocity of the diffracted wave, wherein the target Green's function represents the sum of the propagation time and Amplitude compensation factor; Constructing the reweighted imaging model based on the target Green's function and the common offset diffraction wave data; Constructing the reweighted imaging model based on the target Green's function and the common offset diffraction wave data includes: by formula Construct the re-weighted imaging model, wherein, w _i is a re-weighting coefficient, G is the matrix form of the target Green's function, r _i is the diffraction wave imaging result r(x of imaging point x _i in the underground imaging space _i ) scalar form, d _obs is the common offset diffraction wave data, i takes 1 to N in turn, and N represents the number of imaging points in the underground imaging space; The preset algorithm includes an adaptive homotopy algorithm, using a preset algorithm to calculate the reweighted imaging model, and using the calculation result as the target imaging result of the diffracted wave includes: performing a superposition operation on the reweighted imaging model by using the adaptive homotopy algorithm, and using the superposition result as the target imaging result; The reweighted imaging model is superposed by using an adaptive homotopy algorithm, and the superimposed result is used as the target imaging result including: Taking the preset initial parameter value of the target parameter as the current parameter value, perform the following steps until the parameter value of the target parameter meets the preset condition, wherein the target parameter includes: the reweighting coefficient, the underground imaging Diffraction wave imaging result of spatial imaging point x _i , iteration termination parameter; The first calculation step, according to the formula Computes the scalar value of the parameter value for the current reweighting factor and computes the current update direction vector Wherein, G _Γ is the matrix that is made up of the target column vector in the matrix G of described target Green's function, and the sequence number of described target column vector in the matrix of described Green's function is corresponding to the index number in the current set Γ, The index number in the current set Γ is composed of the sequence number corresponding to the non-zero value in the current inversion solution vector r( _xi ), the current inversion solution vector r( _xi ) is composed of r _i , and the diagonal matrix W and diagonal The diagonal elements of are respectively composed of w _i and composition; The second calculation step, according to the formula Calculate the current update step size Δr _i , where s _i is the ith element of the current update direction vector s; The third calculation step, according to the formula r _i := _ri +(Δr _i )s _i and the formula Calculate the result of the current iteration; The first update step is used to delete the element corresponding to i in the current set Γ when it is judged that Δr _i <1, or, when it is judged that Δr _i ≥ 1, delete the element corresponding to i in the current set Γ Add a new index number in; The second update step, according to the formula i∈Γ ^c updates the parameter value of the current reweighting coefficient; Judging step, judging whether the updated current reweighting coefficient and the current iteration termination parameter satisfy the preset condition, wherein the preset condition is max(w _i )≤τi=1,2,...,N holds true , or, the current iteration termination parameter is greater than or equal to the target threshold; Wherein, if it is judged that the preset condition is satisfied, the current inversion solution vector r( _xi ) is output, and if it is judged that the preset condition is not satisfied, the parameter value controlling the termination parameter of the current iteration is increased Preset values, and use the parameter value of r _i after iteration in the third calculation step and the parameter value of the current reweighting coefficient after updating in the second update step as the current parameter value, and return to execute The first calculation step. 2. The method according to claim 1, characterized in that, performing data preprocessing on the acquired initial shot data, obtaining common offset diffraction wave data comprises: Filtering the initial shot data to obtain co-offset shot data, wherein the co-offset shot data have the same offset distance; According to the sparse Radon hyperbolic transformation method, the common offset shot set data is transformed to obtain the transformed Radon domain; cutting off the part corresponding to the frequency spectrum of the reflected wave in the Radon domain; Inverse-sparse Radon hyperbolic transformation is performed on the cut-off Radon field to obtain the common-offset diffraction wave data. 3. An imaging device for diffracted waves, comprising: An acquisition unit, configured to acquire initial shot data, wherein the initial shot data carries geological information in the target area, and the geological information includes at least one of the following: geological information of stratum layer structure, geological information of fault form information, geological information of karst caves; A processing unit, configured to perform data preprocessing on the acquired initial shot data to obtain common-offset diffraction wave data, wherein the common-offset diffraction wave data have the same offset; A construction unit, configured to construct a reweighted imaging model based on the migration velocity of the diffraction wave in the common offset diffraction wave data and the common offset diffraction wave data; a calculation unit, configured to use a preset algorithm to calculate the reweighted imaging model, and use the calculation result as the target imaging result of the diffracted wave; The first calculation subunit is used to calculate the target Green's function according to the migration velocity of the diffracted wave, wherein the target Green's function represents the position of the diffracted wave from the shot point position to any imaging point in the underground imaging space The travel time and amplitude compensation factor to the geophone location; Constructing a subunit for constructing the reweighted imaging model based on the target Green's function and the co-offset diffraction wave data; Building blocks for passing formulas Constructing the re-weighted imaging model, w _i is the re-weighting coefficient, G is the matrix form of the target Green's function, _ri is the diffraction wave imaging result r( _xi ) of the imaging point x _i in the underground imaging space scalar form, d _obs is the common offset diffraction wave data, i takes 1 to N in turn, and N represents the number of imaging points in the underground imaging space; The second calculation subunit is used to perform superposition operation on the reweighted imaging model by using an adaptive homotopy algorithm, and obtain a result after superposition as the target imaging result; The second calculation subunit includes: taking the preset initial parameter value of the target parameter as the current parameter value, and performing the following steps until the parameter value of the target parameter meets the preset condition, wherein the target parameter includes: reweighting coefficient, underground imaging space The diffraction wave imaging result of the imaging point xi, the iteration termination parameter; the first calculation module is used to follow the formula Computes the scalar value of the parameter value for the current reweighting factor and computes the current update direction vector Among them, G _Γ is a matrix composed of target column vectors in the matrix G of the target Green's function, the serial number of the target column vector in the matrix of Green's functions corresponds to the index number in the current set Γ, and the index in the current set Γ The number is composed of the sequence number corresponding to the non-zero value in the inversion solution vector r( _xi ), the inversion solution vector r( _xi ) is composed of r _i , the diagonal matrix W and the diagonal matrix The diagonal elements of are respectively composed of w _i and Composition; the second calculation module is used to follow the formula Calculate the current update step size Δr _i , where s _i is the ith element of the current update direction vector s; the third calculation module is used to follow the formula r _i := r _i +(Δr _i )s _i and the formula Calculate the current iteration result; the first update module is used to delete the element corresponding to i in the current set Γ when it is judged that Δr _i <1 is true, or when it is judged that Δr _i ≥ 1 is true, in the current Add a new index number in the set Γ; the second update module is used to follow the formula ^i∈Γc updates the parameter value of the current reweighting coefficient; the judging module is used to judge whether the updated current reweighting coefficient and the current iteration termination parameter meet the preset condition, wherein the preset condition is max(w _i )≤τ is established, or the current iteration termination parameter is greater than or equal to the target threshold, i=1,2,...,N; wherein, if it is judged that the preset condition is satisfied, the current inversion solution vector r( _xi ) is output, and if it is judged that If the preset condition is not satisfied, the parameter value of the current iteration termination parameter is controlled to increase the preset value, and the parameter value of r _i after iteration in the third calculation step and the parameter of the current reweighting coefficient after updating in the second update step value as the current parameter value, return to perform the first calculation step.