Lossy Compression Method for Seismic Data Based on Blind Inversion Theory

Abstract

The present invention provides a lossy compression method for seismic data based on the blind inversion theory. The lossy compression method for seismic data based on the blind inversion theory includes: Step 1, inputting seismic data corresponding to the current shot number; Step 2, extracting local seismic data within the spatial window of the current gather based on the preset spatial window ranges in the X and Y directions; Step 3, performing compression on the local seismic data based on the blind inversion theory; Step 4, determining whether the sliding window has traversed the current gather; Step 5, determining whether all shot gathers have been processed. The lossy compression method for seismic data based on the blind inversion theory utilizes the local time-space information of multi-channel seismic signals to achieve the compressed storage of local characteristic reflection event axes, and the data compression ratio is higher than that of general lossless compression technologies, which can meet the requirements of real-time data processing.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration data processing and interpretation, and particularly to a method for lossy compression of seismic data. Background Art

[0002] The combination of the efficient acquisition technology of vibrator sources and the high-density observation system has greatly improved the shot gather density and the seismic data acquisition efficiency. With the normalization of high shot gather density and efficient acquisition operations, it brings a huge amount of seismic data (PB level), posing a great challenge to the real-time quality analysis and monitoring of the seismic data collected in the field. The conventional lossless data storage method can no longer meet the real-time processing requirements of the massive seismic data. It is necessary to study the lossy compression technology for the massive seismic data, extract the main reflection event axes in the seismic data for real-time imaging of the subsurface structure, so as to realize the on-site data quality monitoring and the evaluation of the observation system.

[0003] Conventional lossless compression technologies, such as the "Method for Lossless Compression and Decompression of High-Efficiency Adaptive Seismic Data Streams" (Publication No.: CN 104378118 A), are characterized by performing lossless compression on the data stream in real time. The original 24-bit 3-byte form of the sampled data is adaptively compressed into 1 byte or 2 bytes or 3 bytes or a small amount of data is converted into 4 bytes by using an encoding method; in the patent "An Adaptive Real-Time Lossless Compression Method for Seismic Data Streams" (Publication No.: CN107135004 A), the seismic data of 3n bytes for n sampling points in a single channel is compressed. The data compression is characterized by two steps: (1) Nth-order differential prediction coding, (2) kth-order exponential Golomb coding. Although the above patents can compress the signal losslessly, the data compression ratio is not high.

[0004] Among many lossy compression technologies, compressive sensing is the most widely used guiding ideology. In the patent "Method and System for Reconstructing Seismic Data Using Compressive Sensing Algorithm" (Publication No.: CN 111428193 A), a method for reconstructing seismic data using compressive sensing algorithm is proposed, mainly including: using an overcomplete dictionary to replace the basis function, performing sparse transformation on the original data; selecting an initialized random unit matrix as the observation matrix; combining the sampling matrix with regularization and then using the orthogonal matching pursuit algorithm to realize the recovery and reconstruction of the missing seismic data.

[0005] In the patent "A Multi-Domain Sparse Seismic Data Reconstruction Method and System Based on Compressed Sensing" (Publication No.: CN 113109866 A), a multi-domain sparse seismic data reconstruction method and system based on compressed sensing are proposed. The method includes: performing data sampling on seismic observation data based on a seismic data observation system to obtain seismic sampling data; performing inversion on the seismic sampling data using an inversion algorithm with sparse constraints to obtain the solution of the sampling data in the sparse transform domain; and obtaining the seismic observation data on a reconstructed regular grid according to the solution of the sampling data in the sparse transform domain. The present invention can obtain seismic data with higher accuracy and higher resolution, and reduce the cost of seismic data acquisition.

[0006] In the patent "Real-Time Compression and High-Precision Reconstruction Method of Seismic Data in Wavelet Domain Based on Compressed Sensing" (Publication No.: CN 107045142 B), a real-time compression and high-precision reconstruction method of seismic data in wavelet domain based on compressed sensing is proposed, including the following steps: First, sparsely represent the microseismic signal in the wavelet domain; then construct a chaotic Bernoulli measurement matrix (CBMM) using a Logistic chaotic sequence, and compressively observe the sparsely represented microseismic signal with the measurement matrix; finally, use the Bayesian wavelet tree structure compressed sensing reconstruction algorithm (BTSWCS) to recover the complete original data.

[0007] In the patent "A Seismic Data Reconstruction Method Based on Spatial Constraint Compressed Sensing" (Publication No.: CN109490957 B), a seismic data reconstruction method based on spatial constraint compressed sensing is proposed, including: using a part of the data as training data, using the K-SVD dictionary learning to train an overcomplete dictionary to reconstruct the original seismic data; using the method of joint sparse decomposition to extract the common spatial information and transform the sensing matrix in the compressed sensing algorithm; improving the sparse degree adaptive matching pursuit algorithm, introducing the method of initial sparse degree estimation, and adopting a variable step size strategy to reconstruct the data. The reconstruction result not only has clearer details, but also the operation time is significantly reduced compared with IRLS and SAMP, and the horizontal transition is smoother, indicating that the algorithm designed in the present invention utilizes the relevant spatial information and the reconstruction result is more real.

[0008] In the patent "A Seismic Data Compression Method Based on Tensor Adaptive Rank Truncation" (Publication No.: CN106646595 A), a seismic data compression method based on tensor adaptive rank truncation is proposed. The truncation rank is set by given compression conditions and the singular values of different dimensions obtained immediately through high-order singular value decomposition, and at the same time, the size of the truncation rank is determined according to the distribution of the singular values for tensor decomposition; while ensuring the compression ratio, the peak signal-to-noise ratio of the compression is improved and the compression effect is improved.

[0009] In the patent "A Method and Device for Compressing Massive Seismic Data while Maintaining Spatial Attribute Information" (Publication Number: CN 103592684 A), a method and device for compressing massive seismic data while maintaining spatial attribute information are proposed. The main feature is that along the direction with the densest spatial distribution of the seismic wave field, the seismic data to be compressed within each grid is compressed. The original seismic data can be effectively compressed by 2 to 4 times, and various attribute information obtained from omnidirectional acquisition can be retained, providing convenience for subsequent series of processing.

[0010] The uniqueness of the technology of the present invention lies in considering the high-dimensional spatial characteristics of seismic signals, and regarding seismic signals as the linear superposition of local coherent axes. Through the blind inversion theory, the reflection wavelet waveform and its kinematic characteristics of the coherent axes are determined to achieve lossy data compression and reconstruction. Therefore, the above existing technologies are quite different from the present invention and fail to solve the technical problems we want to solve. There is an urgent need to form a lossy seismic data compression method with high computational efficiency and high compression ratio to meet the real-time data processing requirements. For this reason, we have invented a seismic data compression method based on the blind inversion theory. Summary of the Invention

[0011] The object of the present invention is to provide a method for lossy compression of seismic data, which extracts the main reflection information from the high-density acquired seismic data for real-time imaging of underground structures, thereby evaluating the quality of on-site data acquisition.

[0012] The present invention can be realized by the following technical measures: The method for lossy compression of seismic data based on the blind inversion theory includes:

[0013] Step 1: Input the seismic data corresponding to the current shot number.

[0014] Step 2: Based on the preset spatial window ranges in the X and Y directions, extract the local seismic data within the spatial window of the current gather.

[0015] Step 3: Achieve lossy compression of the local seismic data based on the blind inversion theory.

[0016] Step 4: Determine whether the sliding window has traversed the current gather. If the current gather has been traversed, output the compressed seismic data; otherwise, slide the spatial window and execute Step 2.

[0017] Step 5: Determine whether all shot gathers have been processed. If all shot gathers have been processed, end the data compression processing flow; otherwise, increment the shot number and execute Step 1.

[0018] The present invention can also be realized by the following technical measures:

[0019] In step 1, it is necessary to loop through the entire common-shot gather and read in the seismic data of one shot at a time. The input data structure is in the standard SEG-Y format, and the trace header keywords to be read are: the x, y, z coordinates of the source and geophones, the trace length, and the sampling interval.

[0020] In step 2, according to the parameters preset before the lossy compression processing of the seismic data, such as the spatial spread range of the sliding spatial window in the X and Y directions and the sliding distance each time (usually half of the sliding spatial window length), the seismic traces and their trace header keywords (i.e., local seismic data) within the local spatial window of the current input gather are extracted.

[0021] In step 3, according to the preset blind inversion parameters, the main reflection event in the local seismic data is compressed, and the lossy compressed local seismic data is output, specifically including steps 301 - 308 listed below:

[0022] In step 301, first, the input data and parameters of the blind inversion method are initialized, such as the maximum number of iterations N_iter, the relative energy ratio of the data residual ε (the ratio of the two-norm of the data residual to the two-norm of the input data), and the input local seismic data d obs = d(x,t) (where x is the spatial coordinate, t is the time, italic d is the seismic signal, and bold d is the vector representation of the seismic signal), and the data residual is initialized as e 0 = d obs (In the first step of the iterative inversion, the data residual vector is the input local seismic data, e 0 represents the residual vector, and the superscript 0 represents the initial value), and the iteration number variable k = 1 is set;

[0023] In step 302, the data residual is projected into the model space (the model space is the transform domain characterized by the intercept time (τ) and local slope (p)) using the following integral transform (i.e., the adjoint transform operator T H ):

[0024] ξ k = ξ k (τ,p) = T H e k-1 ≡ ∫e k-1 (τ + p·x)dx

[0025] In the above formula and all subsequent variables or vectors, the superscript k represents the k-th iteration.

[0026] In step 303, the discrete intercept time and slope are traversed to search for the discretized intercept time (τ k ) and local slope p k corresponding to the maximum value in the model space. The formula is expressed as:

[0027]

[0028] In step 304, extract the model space wavelet waveform w corresponding to (τ k , p k ), which is expressed by the formula: k as follows:

[0029] w k = w k (t) = H(t; τ k ) ξ k (τ k , p k )

[0030] where H(t; τ k ) is a one-dimensional square wave function H(t) with the center point at τ k time (the length of the square wave function can be automatically determined according to the data frequency band, and the two sides of the square wave function are usually attenuated by a windowing function to reduce the Gibbs effect)

[0031] In step 305, use the following prediction operator T to map the wavelet in the model space to the data space to obtain the prediction result of the local feature reflection event

[0032]

[0033] In step 306, adaptively subtract the current data residual from the prediction result to update the residual vector e k in the data space:

[0034]

[0035] where α is an adaptive amplitude correction factor that varies with the spatial coordinate and is expressed as

[0036]

[0037] In step 307, determine whether the iteration termination condition is satisfied. The preset iteration termination conditions include: (1) the number of iterations exceeds the maximum number of iterations (i.e., determine whether k + 1 is greater than N_iter); (2) the relative energy ratio of the data residual is less than the preset error energy ratio ε (i.e., is true, e 0 is the initial residual vector, and e k is the residual vector in the k-th iteration process). If either of the two termination conditions is satisfied, the iteration is terminated; otherwise, increment the iteration number variable k = k + 1.

[0038] In step 308, output the compressed data represented in the model space, denoted as W = [w1 , w 2 ,..., w k = [w 1 (t), w 2 (t), K, w k (t)].

[0039] In step 4, it is judged whether the sliding spatial window has traversed the current gather. If not, the sliding spatial window (the default sliding length is half of the spatial window length, which can be overridden by the user's input parameter) is slid and step 2 is continued to be executed.

[0040] In step 5, it is judged whether all shot gathers have been processed. If not, step 1 is continued to be executed. Otherwise, the data processing flow is ended and the lossy compression of seismic data is completed.

[0041] The lossy compression method of seismic data based on the blind inversion theory in the present invention realizes sparse representation of the main characteristic reflection event axes in the high-density acquired seismic signals, and simultaneously inverses the dictionary atoms and the corresponding coefficient vectors. Therefore, the method of the present invention belongs to the solution of non-convex optimization problems under the blind inversion framework in mathematics, and there is no unique solution in theory. Among the numerous non-unique solutions, the method of the present invention selects the characteristic solution that can best represent the prestack seismic signals, that is, uses the linear combination of local reflection event axes to approximate the original signal (i.e., the input local seismic data), thereby realizing the compressed storage of the original data. At the same time, since the data compression process is realized in the model space (representing information such as the travel time and slope of the signal), the kinematic attributes (the arrival time and slope corresponding to each reflection event axis) and dynamic information (the waveform characteristics of the reflection wavelet) of the reflection seismic signals are decoupled while the data is compressed, providing input data for fast prestack migration imaging (real-time imaging) and meeting the requirements of real-time data processing. Brief Description of the Drawings

[0042] Figure 1 is a flowchart of a specific embodiment of the lossy compression method of seismic data based on the blind inversion theory of the present invention;

[0043] Figure 2 is a flowchart of the specific implementation steps of the lossy compression of local seismic data in a specific embodiment of the present invention;

[0044] Figure 3 is a schematic diagram of the input local seismic data in a specific embodiment of the present invention;

[0045] Figure 4 is a schematic diagram of the change of the model space in each iteration process in a specific embodiment of the present invention;

[0046] Figure 5is a schematic diagram of wavelets in the model space extracted during each iteration in a specific embodiment of the present invention;

[0047] Figure 6 is a schematic diagram of the local reflection events predicted during each iteration in a specific embodiment of the present invention;

[0048] Figure 7 Schematic diagram of the change of data residual during each iteration in a specific embodiment of the present invention;

[0049] Figure 8 A comparison diagram of input data, reconstruction error, and reconstruction results in a specific embodiment of the present invention;

[0050] Figure 9 A schematic diagram of local high-dimensional seismic data (actual data from a certain three-dimensional work area) in a specific embodiment of the present invention;

[0051] Figure 10 Schematic diagram showing the comparison of data compression results (tracks 1-50), reconstructed seismic trace (track 51), reference seismic trace (track 52), and reconstruction error (track 53) in a specific embodiment of the present invention;

[0052] Figure 11 Schematic diagram showing the comparison of input data, reconstruction error and reconstruction result in a specific embodiment of the present invention. DETAILED DESCRIPTION

[0053] It should be noted that the following detailed descriptions are exemplary and intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present invention belongs.

[0054] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present invention. As used herein, unless the context clearly indicates otherwise, the singular form is intended to include the plural form. In addition, it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations and / or combinations thereof.

[0055] The present invention discloses a lossy compression method for seismic data based on blind inversion theory. By achieving sparse expression of characteristic reflection events in high-density acquired seismic signals under a blind inversion framework, the method decouples the kinematic properties and dynamic information of the reflected seismic signals while compressing the data. This provides input data for fast prestack migration imaging (i.e., real-time imaging), meeting the real-time processing requirements of seismic data.

[0056] likeFigure 1 As shown in FIG, the seismic data lossy compression method based on blind inversion theory includes the following data processing procedures:

[0057] Step 1: Input the seismic data corresponding to the current shot number;

[0058] Step 2: extracting local seismic data of the current gather within the spatial window based on the pre-set spatial window range in the X and Y directions;

[0059] Step 3: compress local seismic data based on blind inversion theory;

[0060] Step 4: Determine whether the sliding window has traversed the current gather. If it has, output the compressed seismic data; otherwise, slide the spatial window and execute step 2.

[0061] Step 5: Determine whether all shot gathers have been processed. If so, the data compression process ends; otherwise, the shot number is incremented and step 1 is executed.

[0062] The following are several specific embodiments of the present invention.

[0063] Example 1:

[0064] In a specific embodiment 1 of the present invention, the seismic data lossy compression method based on blind inversion theory includes the following steps:

[0065] In step 1, you need to loop through the entire common shot gather, reading in seismic data for one shot at a time. The input data is in the standard SEG-Y format, and the trace header keywords to be read are: the x, y, z coordinates of the source and receiver, the trace length, and the sampling interval.

[0066] In step 2, based on the parameters pre-set before the lossy compression processing of the seismic data, such as the spatial distribution range of the sliding spatial window in the X and Y directions and the sliding distance each time (usually half the length of the sliding spatial window), the seismic traces and their trace header keywords (i.e., local seismic data) of the current input data set within the local spatial window are extracted.

[0067] In step 3, the main reflection events in the local seismic data are compressed according to pre-set blind inversion parameters, and the lossy compressed local seismic data are output.

[0068] In step 4, it is determined whether the sliding spatial window has traversed the current gather. If not, the spatial window is slid (the default sliding length is half the spatial window length, which can be overridden by user input parameters) and the execution continues with step 2.

[0069] In step 5, it is judged whether all shot gathers have been processed. If not, step 1 is continued. Otherwise, the data processing flow is ended, and the lossy compression of seismic data is completed.

[0070] Embodiment 2:

[0071] In a specific embodiment 2 of applying the present invention, as Figure 2 shown, taking a two-dimensional synthetic seismic record as an example, the lossy compression processing flow and specific implementation steps of local seismic data are introduced in detail:

[0072] In step 301, first, the input data and parameters of the blind inversion method are initialized, such as the maximum number of iterations N_iter, the relative energy ratio ε of the data residual (the ratio of the two-norm of the data residual to the two-norm of the input data), the input local seismic data d obs = d(x,t) (as shown in the appendix Figure 3 shown, including 8 local linear event axes, a total of 21 traces, and 2000 samples per trace), and the data residual e 0 = d obs (in the first step of iterative inversion, the data residual vector is the input local seismic data), and the initial value of the iteration number variable k = 1 is set;

[0073] In step 302, the following integral transform (i.e., the adjoint transform operator T H ) is used to project the data residual into the model space (the model space is the transform domain characterized by the intercept time (τ) and local slope (p), as shown in the appendix Figure 4 shown):

[0074] ξ k = ξ k (τ,p) = T H e k-1 ≡ ∫ e k-1 (τ + p·x)dx

[0075] In step 303, the discrete intercept time and slope are traversed to search for the discretized intercept time (τ k ) and local slope p k corresponding to the maximum value in the model space, and the formula is expressed as:

[0076]

[0077] In step 304, the model space wavelet waveform w k ,p k ) corresponding to (τ k ) is extracted (as shown in the appendix Figure 5 shown), and the formula is expressed as:

[0078] w k = wk H(τ k )ξ k (τ k ,p k )

[0079] wherein, H(τ k ) is a one-dimensional square wave function (the center point is located at τ k time, the length of the square wave function can be automatically determined according to the data frequency band, and the two sides of the square wave function are usually attenuated by a windowing function to reduce the Gibbs effect)

[0080] In step 305, the wavelet in the model space is mapped to the data space by using the following prediction operator T to obtain the prediction result of the local feature reflection event (as shown in the appendix Figure 6 ):

[0081]

[0082] In step 306, the current data residual is adaptively subtracted from the prediction result to update the residual vector in the data space (as shown in the appendix Figure 7 ):

[0083]

[0084] wherein, α is an adaptive amplitude correction factor, expressed as

[0085] In step 307, it is judged whether the iteration termination condition is satisfied. The preset iteration termination conditions include: (1) the number of iterations exceeds the maximum number of iterations (that is, judge whether k + 1 is greater than N_iter); (2) the relative energy ratio of the data residual is less than the preset energy ratio (that is is established). If any of the two termination conditions is satisfied, the iteration is terminated; otherwise, the iteration number variable is incremented by k = k + 1.

[0086] In step 308, the compressed data represented in the model space is output, expressed as W = [w 1 ,w 2 ,...,w k = [w 1 (t),w 2 (t),K,w k (t)] (as shown in the appendix Figure 5 , the compressed data only contains 8 wavelets).

[0087] In order to verify the invention effect, the effectiveness of the present invention can be demonstrated by comparing the data reconstruction error. The input data, the reconstruction error, and the reconstruction result are combined together, as shown in the appendix Figure 8As shown, the amplitude of the reconstruction error is almost negligible compared to the input data, indicating that the compressed data of the present invention can reconstruct the input data well. In addition, the size of the original input data is: 21 traces * 2000 floating points / trace * 4 bytes / floating point = 168000 bytes (ignoring trace header keywords). The compressed data only contains 8 wavelets, and each wavelet only retains the values within the wavelet length (i.e., 200 floating points). Therefore, the size of the compressed data is: 8 traces * 200 floating points / trace * 4 bytes / floating point = 6400 bytes (ignoring trace header keywords). The lossy data compression ratio is: size before compression / size after compression = 168000 bytes / 6400 bytes = 26.25.

[0088] Example 3:

[0089] Next, taking local high-dimensional seismic data as an example, a data compression example of the lossy compression processing technology for local seismic data in a high-dimensional space will be introduced in detail. As Figure 9 shown, the local high-dimensional seismic data is expressed as: do bs = d(t, s, r), |s - s0| < W s && |r - r0| < W r (where s and r are the spatial coordinates of the source and the geophone respectively, s0 and r0 are the coordinates of the central source and the geophone respectively, and Ws and Wr are the local spatial window lengths at the source and geophone ends). By executing the lossy compression processing flow of local seismic data of the present invention and making the following generalizations to the integral transform operator and the prediction operator:

[0090] (1) Replace the integral transform in step 303 with a high-dimensional local linear integral transform, expressed as:

[0091] ξ k = ξ k (τ, p s , p r ) = T H e k-1 ≡ ∫∫ e k-1 (τ + p s ·(s - s0) + p r ·(r - s r )) ds dr

[0092] In the above integral transform, ps and pr are the slopes of the local plane waves at the source and geophone ends respectively.

[0093] (2) Replace the slope p in steps 303 and 304 with ps and pr;

[0094] (3) Replace the prediction operator in step 305 with a high-dimensional local linear prediction operator, expressed as:

[0095]

[0096] After the above promotion, the present invention can be used to compress local high-dimensional seismic data. As Figure 10 shown, seismic traces 1-50 represent Figure 9 50 compressed-domain wavelets retained after compression of the input data shown. Summing these 50 compressed wavelets can reconstruct the central trace (i.e., the central seismic trace corresponding to s0, r0), and the reconstruction result is as shown in Figure 10 trace 51 in. Compared with the center of the original data to ( Figure 10 trace 52 in), there is basically no visual difference between the two. The residual between the reconstruction result and the input data is as shown in Figure 10 trace 53 in. It can be seen that there is no obvious effective signal residue in the central trace reconstruction error. The reconstruction result and reconstruction error of the entire input data are as shown in Figure 11 shown, indicating that the present invention can be naturally extended to the compression and reconstruction of local high-dimensional data.

[0097] The present invention provides a lossy compression method for seismic data based on blind inversion theory. By using a sliding spatial window to extract local seismic data from the input seismic trace set and performing lossy compression on the local seismic data. The lossy compression process mainly includes: initializing blind inversion parameters; projecting the data residual to the model space using integral transformation (the model space is the transformed domain characterized by intercept time and local slope); traversing the discrete intercept time and slope, and searching for the discrete intercept time and local slope corresponding to the maximum value in the model space; extracting the wavelet waveform in the model space corresponding to the intercept time and local slope; using a prediction operator to map the wavelet in the model space to the data space to obtain the prediction result of the local characteristic reflection event; adaptively subtracting the current data residual from the prediction result to update the residual vector in the data space; determining whether the iteration termination condition is satisfied (if the iteration termination condition is satisfied, output the compressed seismic data characterized in the model space). This method utilizes the local time-space information of multi-channel seismic signals to achieve the compressed storage of local characteristic reflection events, and the data compression ratio is higher than that of general lossless compression techniques, which can meet the requirements of real-time data processing.

[0098] Finally, it should be noted that the above are only preferred embodiments of the present invention and are not used to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or perform equivalent replacements for some of the technical features. For example: (1) Although Example 2 provided by the present invention takes a common shot gather as an example to show the specific implementation process of lossy compression of seismic data. However, the input data of this method is not limited to common shot gathers, and any sorting method of seismic data, such as common geophone gathers, common midpoint gathers, common offset gathers, etc., can extract local seismic data by means of a sliding spatial window and use the present invention to achieve lossy compression of local seismic data; (2) Although the adjoint transform operator T H in step 302 and the prediction operator T in step 305 are both two-dimensional local linear transforms, other forms of mathematical transforms, including but not limited to high-dimensional local linear integral transforms (see Example 3 of the present invention), hyperbolic integral transforms, parabolic integral transforms, polynomial integral transforms, and even any family of basis functions as integral transform kernels, etc., all fall within the scope of this patent. Therefore, any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

[0099] Except for the technical features described in the specification, they are all well-known technologies to those skilled in the art.

Claims

1. A lossy compression method for seismic data based on the blind inversion theory, characterized in that The lossy compression method for seismic data based on the blind inversion theory includes: Step 1: Input the seismic data corresponding to the current shot number; Step 2: Extract the local seismic data within the spatial window of the current gather based on the preset spatial window ranges in the X and Y directions; Step 3: Compress the local seismic data based on the blind inversion theory; Step 4: Determine whether the sliding window has traversed the current gather; Step 5: Determine whether all shot gathers have been processed; In Step 3, based on the preset blind inversion parameters, compress the main reflection event axes in the local seismic data and output the lossy compressed local seismic data; Step 3 includes: Step 301: Initialize the input data and parameters of the blind inversion method; Step 302, using the following integral transform, i.e., the adjoint transform operator T H , project the data residual into the model space, where the model space is the transform domain characterized by the intercept time τ and the local slope p: ξ k = ξ k (τ, p) = T H e k-1 ≡ ∫ e k-1 (τ + p·x)dx Step 303: Traverse the discrete intercept times and slopes, and search for the discretized intercept time τ corresponding to the maximum value in the model space k and the local slope p k , which is expressed by the formula as: Step 304, extract the model space wavelet waveform w corresponding to (τ k , p k ), which is expressed by the formula as: k ​ w k = w k (t) = H(t; τ k )ξ k (τ k , p k ) Among them, H(t; τ k ) is a one-dimensional square wave function H(t) with its center point located at τ k at the time. The length of the square wave function is automatically determined according to the data frequency band, and the two sides of the square wave function are attenuated by a windowing function to reduce the Gibbs effect; Step 305: Use the prediction operator T to map the wavelet in the model space to the data space to obtain the prediction result of the local characteristic reflection event axis: Step 306: Subtract the current data residual from the prediction result adaptively to update the residual vector in the data space; Step 307: Determine whether the iteration termination condition is satisfied; Step 308, output the compressed data represented in the model space, denoted as W = [w 1 , w 2 ,..., w k = [w 1 (t), w 2 (t), …, w k (t)].

2. The lossy compression method for seismic data based on the blind inversion theory according to claim 1, wherein In Step 1, it is necessary to loop through the entire common shot gather and read in the seismic data of one shot at a time; the input data structure is in the standard SEG-Y format, and the trace header keywords to be read are: the x, y, z coordinates of the source and geophones, the trace length, and the sampling interval.

3. The lossy compression method of seismic data based on the blind inversion theory according to claim 1, wherein In Step 2, according to the parameters preset before the lossy compression processing of the seismic data, including the spatial spread ranges of the sliding spatial window in the X and Y directions and the sliding distance each time, extract the seismic traces and their trace header keywords within the local spatial window of the current gather.

4. The lossy compression method for seismic data based on the blind inversion theory according to claim 3, wherein, In Step 2, the sliding distance each time is half of the length of the sliding spatial window.

5. The lossy compression method for seismic data based on the blind inversion theory according to claim 1, wherein In step 301, the input data and parameters of the blind inversion method are first initialized, including the maximum number of iterations N_iter, the relative energy ratio ε of the data residual, that is, the ratio of the two-norm of the data residual to the two-norm of the input data, and the input local seismic data d obs = d (x,t) , where t is time, italic d is the seismic signal, bold d is the vector representation of the seismic signal, and the data residual is initialized as e 0 = d obs , in the first step of the iterative inversion, the data residual vector is the input local seismic data, e 0 represents the residual vector, and the superscript 0 represents the initial value.

6. The lossy compression method for seismic data based on the blind inversion theory according to claim 1, wherein In Step 306, the formula for updating the residual vector in the data space is: where α is the adaptive amplitude correction factor that varies with the spatial coordinates and is expressed as 7. The lossy compression method of seismic data based on the blind inversion theory according to claim 1, wherein In step 306, the preset iteration termination conditions include: (1) when the number of iterations exceeds the maximum number of iterations, it is judged whether k + 1 is greater than N_iter; (2) the relative energy ratio of the data residual is less than the preset error energy ratio ε, that is Whether it holds, e 0 Is the initial residual vector, e k Is the residual vector in the k-th iteration process; if any of the two termination conditions is satisfied, the iteration is terminated; otherwise, the iteration number variable is incremented by k = k + 1.

8. The lossy compression method for seismic data based on the blind inversion theory according to claim 1, wherein In Step 4, determine whether the sliding spatial window has traversed the current gather; if not, slide the spatial window and continue to execute Step 2.

9. The lossy compression method of seismic data based on the blind inversion theory according to claim 1, wherein, In Step 5, determine whether all shot gathers have been processed; if not, continue to execute Step 1; otherwise, end the data processing flow and complete the lossy compression of the seismic data.

Images