Seismic record calibration method and apparatus, and server and storage medium

Abstract

Provided in the present application are a seismic record calibration method and apparatus, and a server and a storage medium. The method comprises: receiving seismic records of a plurality of work areas that are collected by seismic sensors, and truncating the seismic records of all the work areas according to a preset time period, so as to obtain seismic data corresponding to all the work areas; determining a plurality of time points within the preset time period, and collecting corresponding target data from the seismic data of each work area on the basis of the plurality of time points; on the basis of the plurality of time points, calculating all item operators corresponding to all the target data; on the basis of the item operators, uniformly performing deconvolution processing on the seismic data of all the work areas, so as to obtain deconvolved seismic data of all the work areas; on the basis of the deconvolved seismic data of each work area, calibrating the seismic records, so as to obtain calibrated seismic records of the plurality of work areas; and outputting the calibrated seismic records of the plurality of work areas. Therefore, the consistency of seismic wavelets is improved, and the resolutions of images corresponding to seismic records are enhanced.

Description

Earthquake record calibration methods, devices, servers, and storage media

[0001] This application claims priority to Chinese Patent Application No. 202411941491.X, filed on December 26, 2024, entitled "Method, Apparatus, Server and Storage Medium for Seismic Record Calibration", the entire contents of which are incorporated herein by reference. Technical Field

[0002] This application relates to the field of seismic record processing, and more particularly to a seismic record calibration method, apparatus, server, and storage medium. Background Technology

[0003] Seismic data distortion refers to the distortion of seismic data caused by changes in the location of the seismic source and receiver, or changes in the type of seismic source, in areas with complex surface structures or drastic changes in surface structure, such as the transition zone between land and sea. This distortion hinders the interpretation of seismic information. Therefore, the accuracy of seismic records is extremely important.

[0004] In existing technologies, traditional seismic record calibration mainly involves setting up quality control guns in each of multiple seismic zones and processing seismic data in each zone using two combinations of excitation and reception methods to calibrate the seismic records.

[0005] However, in the existing technology, seismic data is processed by quality control guns in each work area according to two excitation and reception combinations to calibrate the seismic records. The quality control guns are located in a single position, which makes it difficult to reflect the overall differences between the work areas. This can easily lead to poor consistency of seismic wavelets, resulting in low resolution of the images corresponding to the seismic records. Summary of the Invention

[0006] This application provides a seismic record calibration method, apparatus, server, and storage medium to solve the problem that calibrating seismic records by statistically analyzing and calculating partial seismic data can easily lead to poor consistency of seismic wavelets, resulting in low resolution of the corresponding seismic record images.

[0007] In a first aspect, this application provides a seismic record calibration method, applied to a server, comprising:

[0008] It receives seismic records from multiple work areas collected by seismic sensors, and extracts seismic records from all work areas according to a preset time period to obtain seismic data for all work areas.

[0009] Multiple time points within the preset time period are determined, and corresponding target data are collected from the seismic data of all work areas based on the multiple time points;

[0010] Calculate the autocorrelation superposition data of each target data based on the multiple time points, and calculate all corresponding item operators based on their respective autocorrelation superposition data;

[0011] Based on the various operators, the seismic data of all work areas are uniformly deconvolved to obtain the deconvolved seismic data of all work areas;

[0012] The seismic records of the multiple work areas are calibrated based on the seismic data after deconvolution of each work area to obtain calibrated seismic records for multiple work areas;

[0013] Output the calibrated seismic records for the multiple work areas.

[0014] Based on the above technical content, multiple time points within a preset time period are determined, and corresponding target data are collected from seismic data of all work areas according to these multiple time points; all operators corresponding to each target data are calculated based on the multiple time points; the seismic data of all work areas are uniformly deconvolved according to each operator to obtain deconvolved seismic data of all work areas; the seismic records of multiple work areas are calibrated based on the deconvolved seismic data of each work area to obtain calibrated seismic records of multiple work areas; and the calibrated seismic records of multiple work areas are output. By considering all operators globally and performing deconvolution processing to calibrate the seismic records, the consistency of seismic wavelets is improved, and the resolution of the corresponding images of the seismic records is enhanced.

[0015] In one possible design, the step of calculating all-term operators corresponding to each target data based on the multiple time points includes: performing a Fourier transform on the respective related superimposed data to obtain the seismic frequency domain; calculating the corresponding logarithmic power spectrum based on the seismic frequency domain; decomposing the logarithmic power spectrum according to a preset Gauss-Seidel iteration method to obtain the response of the corresponding all terms; performing exponential extraction processing on the response of the all terms to obtain the corresponding all-term power spectrum; and calculating the corresponding all-term operators based on the all-term power spectrum.

[0016] Furthermore, by calculating all the operators corresponding to each target data at multiple time points, the wavelet differences caused by different excitation conditions, reception conditions, combination conditions and geological conditions in each work area are eliminated, thereby improving the consistency of subsequent wavelets.

[0017] In one possible design, the formula for calculating the seismic frequency domain by performing Fourier transforms on the target data based on the multiple time points is as follows:

[0018] X(f)=∫x(t)e -i2πft dt

[0019] In the formula, X(f) is the earthquake frequency domain; x(t) is the earthquake data; e -i2πft Let be a complex exponential function; i be the imaginary unit; t be the time points, where t = t1, t2, ..., tt3. n ,t1,t n These are the start and end times of the preset time period, respectively.

[0020] Accordingly, the formula for calculating the corresponding logarithmic power spectrum based on the earthquake frequency domain is:

[0021] C(f)=ln‖X(f)‖ 2

[0022] In the formula, C(f) is the logarithmic power spectrum; X(f) is the seismic frequency domain;

[0023] Accordingly, all of the terms are global terms, shot point terms, receiver point terms, shot-receiver midpoint terms, and shot-receiver distance terms; the calculation formula for decomposing the logarithmic power spectrum according to a preset consistency method to obtain the response of all corresponding terms is:

[0024] In equation (1), The response to the global item; k is the flag for each work zone; Ng is the maximum number of channels for all data; Logarithmic power spectrum; For the response of the shot point item; The response of the detector point term; The response to the midpoint item of the shot inspection; is the response of the shot-receiver distance term; (n) is the result of the nth iteration; (n-1) is the result of the (n-1)th iteration; i is the shot number in the shot set; j is the receiver point number;

[0025] In equation (2), The response of the shot point term; Ns is the largest path number in the shot set i;

[0026] In equation (3), The response of the detector point term; Nr is the largest channel number in the detector point set j;

[0027] In equation (4), The response of the shot-receiver midpoint term; Nm is the largest channel number in CMP(i+j) / 2;

[0028] In equation (5), The response to the shot-receiver distance term; Nh is the maximum track number in the shot-receiver distance (ij) / 2;

[0029] Accordingly, the formula for calculating the power spectrum of all terms by performing exponential extraction on the responses of all terms is as follows:

[0030] In the formula, G k (f) represents the global power spectrum; k is the identifier for each work zone; S i (f) represents the power spectrum of the shot point term; R j (f) is the power spectrum of the detector term; M (i+j) / 2 (f) is the power spectrum of the midpoint term of the shot-receiver; H (i-j) / 2 (f) represents the power spectrum of the shot-receiver distance term.

[0031] In one possible design, the step of calculating the corresponding term operator based on the power spectrum of all terms includes: determining the global term wavelet based on the global term power spectrum, and calculating the global term operator based on the global term wavelet; performing an inverse Fourier transform on the remaining term power spectra in the power spectrum of all terms to obtain the corresponding term autocorrelation; and calculating the corresponding term operator based on the term autocorrelation.

[0032] Furthermore, by calculating the corresponding operators based on the power spectrum of all terms, the consistency of subsequent wavelets is improved.

[0033] In one possible design, the formula for calculating the autocorrelation of the remaining power spectra in all the power spectra to obtain the corresponding autocorrelation is as follows:

[0034] In the formula, R s (t) represents the autocorrelation of the shot point term; R r (t) represents the autocorrelation of the detector term; R m (t) represents the autocorrelation of the midpoint term in the shot-receiver test; R h (t) Autocorrelation of the shot-receiver distance term; S(f) is the logarithmic power spectrum of the shot term; R(f) is the logarithmic power spectrum of the receiver term; M(f) is the logarithmic power spectrum of the shot-receiver midpoint term; H(f) is the logarithmic power spectrum of the shot-receiver distance term; e -i2πft Let be a complex exponential function; i be the imaginary unit; t be the time points, where t = t1, t2, ..., tt3. n ,t1,t n These are the start and end times of the preset time period, respectively.

[0035] In one possible design, the step of calculating the global term operator based on the global term wavelet includes: processing a preset expected wavelet according to a preset optimization method to obtain an optimal expected wavelet; performing correlation processing on the global term wavelet and the optimal expected wavelet to obtain the cross-correlation between the global term wavelet and the optimal expected wavelet; and calculating the global term operator based on the cross-correlation between the global term wavelet and the optimal expected wavelet.

[0036] Furthermore, by calculating the global term operator based on the global term wavelet, the consistency of subsequent wavelets is improved.

[0037] In one possible design, the formula for calculating the cross-correlation between the global term wavelet and the optimal expectation wavelet is as follows:

[0038] In the formula, the The cross-correlation between the optimal expected wavelet and the global term wavelet; w(i) is the global term wavelet; The optimal expected wavelet; i and j are the time point indices determined from each time point t;

[0039] Accordingly, the calculation formula for the global term operator based on the cross-correlation of the global term wavelet and the optimal expectation wavelet is as follows:

[0040] In the formula, f k (0) to f k (L-1) represents the elements of the global term operator; k is the identifier for each work zone; to These are the elements in the cross-correlation between the optimal expected wavelet and the global term wavelet.

[0041] In one possible design, the step of calculating the corresponding operators based on the autocorrelation of each item includes: calculating the shot point item operator based on the autocorrelation of the shot point item; calculating the receiver item operator based on the autocorrelation of the receiver item; calculating the shot-receiver midpoint item operator based on the autocorrelation of the shot-receiver midpoint item; and calculating the shot-receiver distance item operator based on the autocorrelation of the shot-receiver distance item.

[0042] In one possible design, the calculation formula for the operator that calculates the shot point term based on the autocorrelation of the shot point term is:

[0043] In the formula, h S (0) to h g (L-1) represents the elements of the shot point operator; R S (τ) to R S (τ+L-1) represents the elements of the autocorrelation of the shot point term.

[0044] Secondly, this application provides a seismic record calibration device, applied to a server, comprising:

[0045] The receiving module is used to receive seismic records from multiple work areas collected by seismic sensors, and to extract the seismic records of all work areas according to a preset time period to obtain the seismic data of all work areas.

[0046] The determination module is used to determine multiple time points within the preset time period, and to collect corresponding target data from the seismic data of all work areas based on the multiple time points;

[0047] The calculation module is used to calculate the autocorrelation superposition data of each target data according to the multiple time points, and to calculate the corresponding all-item operators according to their respective autocorrelation superposition data;

[0048] The processing module is used to perform deconvolution processing on the seismic data of all work areas according to various operators, so as to obtain the deconvolutioned seismic data of all work areas.

[0049] The calibration module is used to calibrate the seismic records of the multiple work areas based on the seismic data after deconvolution of each work area, so as to obtain the calibrated seismic records of the multiple work areas.

[0050] The output module is used to output the calibrated seismic records of the multiple work areas.

[0051] Thirdly, this application provides a server, including: at least one processor and a memory;

[0052] The memory stores computer-executed instructions;

[0053] The at least one processor executes computer execution instructions stored in the memory, causing the at least one processor to perform the seismic record calibration method as described in the first aspect and various possible designs of the first aspect.

[0054] Fourthly, this application provides a computer storage medium storing computer execution instructions, which, when executed by a processor, implement the seismic record calibration method described in the first aspect and various possible designs of the first aspect.

[0055] Fifthly, this application provides a computer program product, including a computer program that, when executed by a processor, implements the seismic record calibration method as described in the first aspect above.

[0056] The seismic record calibration method, apparatus, server, and storage medium provided in this application receive seismic records from multiple work areas acquired by seismic sensors, and extract seismic records from all work areas according to a preset time period to obtain seismic data corresponding to all work areas; determine multiple time points within the preset time period, and collect corresponding target data from the seismic data of all work areas based on the multiple time points; calculate all terms of operators corresponding to each target data based on the multiple time points; uniformly perform deconvolution processing on the seismic data of all work areas based on each operator to obtain deconvolutioned seismic data of all work areas; calibrate the seismic records of multiple work areas based on the deconvolutioned seismic data of each work area to obtain calibrated seismic records of multiple work areas; and output the calibrated seismic records of multiple work areas. By considering all terms of operators globally and performing deconvolution processing to calibrate the seismic records, the consistency of seismic wavelets is improved, and the resolution of the corresponding images of the seismic records is enhanced. Attached Figure Description

[0057] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0058] Figure 1 is a schematic diagram of the application scenario of the seismic record calibration method provided in the embodiment of this application;

[0059] Figure 2 is a schematic diagram of the seismic record calibration method provided in an embodiment of this application;

[0060] Figure 3 is a schematic diagram of the seismic record calibration method provided in the embodiment of this application (II).

[0061] Figure 4 is a schematic diagram of the structure of the seismic record calibration device provided in the embodiment of this application;

[0062] Figure 5 is a schematic diagram of the hardware structure of the server provided in an embodiment of this application. Detailed Implementation

[0063] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0064] In regions with undulating terrain, complex surface structures, or drastically changing surface structures such as the transition zone between land and sea, variations in excitation and reception conditions caused by changes in the location of the seismic source and receiver, or changes in the type of seismic source, can distort seismic data, hindering the interpretation of seismic information. Therefore, the accuracy of seismic records is crucial. Current technology primarily involves setting up quality control shots in each of multiple seismic monitoring areas and processing seismic data in each area using two different excitation and reception combinations to calibrate the seismic records. However, this method of using quality control shots to process seismic data in each area using two different excitation and reception combinations is limited by the single location of the quality control shots. It fails to reflect the overall differences between various monitoring areas, easily leading to poor consistency of seismic wavelets and resulting in low resolution of the corresponding seismic record images.

[0065] To address the aforementioned technical problems, this application proposes the following technical concept: Considering the seismic records collected from multiple work areas, the inventors collect corresponding target data from the seismic records of each work area based on multiple time points, calculate all terms of operators corresponding to each target data using multiple time points, and perform deconvolution processing on the seismic data of each work area according to each operator to calibrate the seismic records of multiple work areas, thereby improving the consistency of seismic wavelets and enhancing the resolution of the images corresponding to the seismic records.

[0066] Figure 1 is a schematic diagram of the application scenario of the seismic record calibration method provided in the embodiments of this application.

[0067] As shown in Figure 1, the scenario includes: earthquake sensor 101, server 102, and display terminal 103.

[0068] Among them, the earthquake sensor 101 can be a single sensor or a cluster of multiple sensors.

[0069] Server 102 can be a standalone server or a cluster of multiple servers.

[0070] Display terminal 103 can be a display screen or a terminal such as a personal computer.

[0071] Server 102 receives seismic records from multiple work areas collected by seismic sensor 101 via a wireless network, and collects corresponding target data from the seismic records of all work areas according to multiple predetermined time points; calculates the autocorrelation superposition data of each target data according to multiple time points, and calculates all corresponding operators based on the respective autocorrelation superposition data; performs deconvolution processing on the seismic data of all work areas according to each operator to obtain deconvolutioned seismic data of all work areas; calibrates the seismic records based on the deconvolutioned seismic data of each work area to obtain calibrated seismic records of multiple work areas; and outputs the calibrated seismic records of multiple work areas to display terminal 103 for display. A detailed embodiment is described below.

[0072] Figure 2 is a schematic flowchart of the seismic record calibration method provided in this embodiment. The execution entity in this embodiment can be the server shown in Figure 1, and no special limitation is made here. As shown in Figure 2, the method includes:

[0073] S201: Receives seismic records from multiple work areas collected by seismic sensors, and extracts seismic records from all work areas according to a preset time period to obtain seismic data for all work areas.

[0074] In this embodiment, the seismic record can be a pre-stack shot gather or a CMP gather.

[0075] Among them, pre-stack shot is a technique used in seismic exploration, mainly used to simulate the propagation process of seismic waves in order to obtain information about underground geological structures.

[0076] CMP gathers refer to the extraction of traces from different shot gathers that share a common center point, forming a new set in seismic data processing.

[0077] In this embodiment, the preset time period is a time window with a high signal-to-noise ratio.

[0078] In addition, various excitation and reception methods can be used to collect pre-stack shot gathers or CMP gathers for multiple work areas under different geological conditions, and the seismic data of different work areas can be marked to complete the subsequent calibration of seismic records for multiple work areas.

[0079] S202: Determine multiple time points within a preset time period, and collect corresponding target data from seismic data of all work areas based on multiple time points.

[0080] S203: Calculate the autocorrelation superposition data of each target data based on multiple time points, and calculate the corresponding all-term operators based on their respective autocorrelation superposition data.

[0081] S204: Perform deconvolution processing on the seismic data of all work areas according to the various operators to obtain the deconvolutioned seismic data of all work areas.

[0082] S205: Based on the seismic data after deconvolution of each work area, calibrate the seismic records of multiple work areas to obtain calibrated seismic records for multiple work areas.

[0083] For example, because the seismic source for ocean data is an air gun and the seismic source for land data is explosive, the difference in their sources results in the land data amplitude being much smaller than the source data amplitude on the profile before surface-consistent deconvolution, with discontinuous phase axes and a clear dividing line. However, after surface-consistent deconvolution, the amplitudes of the land data and the source data are basically consistent.

[0084] For example, before surface-consistent deconvolution, land data significantly lacks low frequencies compared to ocean data, indicating inconsistency in their wavelets. After surface-consistent deconvolution, the low frequencies of the land data are compensated, the bandwidth widens, and the spectrum becomes essentially consistent with that of the ocean data. This demonstrates that surface-consistent deconvolution improves wavelet consistency.

[0085] For example, in the autocorrelation data before surface-consistent deconvolution, the amplitudes of land and source data are much smaller. Near time zero, the autocorrelation phase axis of land data is thinner and the phase is different from that of ocean data, which also reflects the difference in wavelets between ocean and land data. However, in the autocorrelation data before surface-consistent deconvolution, the amplitude, phase, and continuity of autocorrelation in both ocean and land data are significantly improved.

[0086] S206: Outputs calibrated seismic records for multiple work areas.

[0087] Specifically, the calibrated seismic records from multiple work areas are output to the display terminal to display the corresponding seismic record images.

[0088] In summary, the seismic record calibration method provided in this embodiment receives seismic records from multiple work areas collected by seismic sensors, and extracts seismic records from all work areas according to a preset time period to obtain seismic data for all work areas; determines multiple time points within the preset time period, and collects corresponding target data from the seismic data of all work areas based on the multiple time points; calculates all terms operators corresponding to each target data based on the multiple time points; performs deconvolution processing on the seismic data of all work areas based on each operator to obtain deconvolutioned seismic data for all work areas; calibrates the seismic records of multiple work areas based on the deconvolutioned seismic data of each work area to obtain calibrated seismic records for multiple work areas; and outputs the calibrated seismic records for multiple work areas. By considering all terms operators globally and performing deconvolution processing to calibrate the seismic records, the consistency of seismic wavelets is improved, and the resolution of the corresponding images of the seismic records is enhanced.

[0089] In addition, the seismic record calibration method provided in this embodiment maintains the consistency of seismic wavelets in each work area by simultaneously calibrating the seismic records of multiple work areas.

[0090] In addition, the seismic record calibration method provided in this embodiment calculates all the operators corresponding to each target data based on multiple time points to eliminate the wavelet differences caused by different excitation conditions, reception conditions, combination conditions and geological conditions in each work area, thereby improving the consistency of subsequent wavelets.

[0091] Figure 3 is a schematic flowchart of the seismic record calibration method provided in this application embodiment. In this application embodiment, based on the embodiment provided in Figure 2, a detailed description is given of the specific implementation method for calculating all corresponding term operators based on their respective related stacking data in step S203. As shown in Figure 3, the method includes:

[0092] S301: Perform Fourier transform on the respective related superimposed data to obtain the seismic frequency domain.

[0093] In this embodiment, Fourier transforms are performed on the target data at multiple time points to obtain the seismic frequency domain. The calculation formula is as follows:

[0094] X(f)=∫x(t)e -i2πft dt

[0095] In the formula, X(f) represents the seismic frequency domain; x(t) represents the seismic data; e -i2πft Let be a complex exponential function; i be the imaginary unit; t be the time points, where t = t1, t2, ..., tt3. n ,t1,t n These are the start and end times of the preset time period, respectively.

[0096] S302: Calculate the corresponding logarithmic power spectrum based on the seismic frequency domain.

[0097] In this embodiment, the formula for calculating the corresponding logarithmic power spectrum based on the seismic frequency domain is as follows:

[0098] C(f)=ln‖X(f)‖ 2

[0099] In the formula, C(f) is the logarithmic power spectrum; X(f) is the seismic frequency domain.

[0100] S303: Decompose the logarithmic power spectrum according to the preset Gauss-Seidel iterative method to obtain the response of all corresponding terms.

[0101] In this embodiment, the default consistency method is the Gauss-Seidel method.

[0102] In this embodiment, all terms include global terms, shot point terms, receiver point terms, shot-receiver midpoint terms, and shot-receiver distance terms; the logarithmic power spectrum is decomposed according to a preset consistency method to obtain the response of all corresponding terms, and the calculation formula is as follows:

[0103] In equation (1), The response for the global item; k is the flag for each work area; Ng is the maximum number of channels for all data; Logarithmic power spectrum; For the response of the shot point item; The response of the detector point term; The response to the midpoint item of the shot inspection; is the response of the shot-receiver distance term; (n) is the result of the nth iteration; (n-1) is the result of the (n-1)th iteration; i is the shot number in the shot set; j is the receiver point number.

[0104] Furthermore, k in the above formula represents different work areas; its global term differs from the conventional surface uniform deconvolution method.

[0105] In equation (2), is the response of the gun point term; Ns is the largest path number in gun set i.

[0106] In equation (3), The response of the detector point term; Nr is the largest channel number in the detector point set j.

[0107] In equation (4), is the response of the shot detection midpoint term; Nm is the largest channel number in CMP(i+j) / 2.

[0108] In equation (5), is the response of the shot-receiver distance term; Nh is the maximum track number in the shot-receiver distance (ij) / 2.

[0109] S304: Perform exponential extraction processing on the responses of all terms to obtain the power spectrum of all terms.

[0110] In this embodiment, the response of all terms is subjected to exponential extraction processing to obtain the power spectrum of all terms. The calculation formula is as follows:

[0111] In the formula, G k (f) represents the global power spectrum; k is the identifier for each work zone; S i (f) represents the power spectrum of the shot point term; R j (f) is the power spectrum of the detector term; M (i+j) / 2 (f) is the power spectrum of the midpoint term of the shot-receiver; H (i-j) / 2(f) represents the power spectrum of the shot-receiver distance term.

[0112] S305: Calculate the corresponding operators for all terms based on the power spectrum of all terms.

[0113] Specifically, step S305 includes:

[0114] S3051: Determine the global term wavelet based on the global term power spectrum, and calculate the global term operator based on the global term wavelet.

[0115] Specifically, step S3051 includes:

[0116] S30511: Process the preset desired wavelet according to the preset optimization method to obtain the optimal desired wavelet.

[0117] In this embodiment, the preset optimization method is the optimization method that maximizes the cross-correlation coefficient (similarity coefficient) of multiple channels.

[0118] In this embodiment, the preset desired wavelet is a representative data segment from the first arrival, overlay profile, or well logging data of the seismic data, or the wavelet is directly provided by the user.

[0119] The criterion for the optimal expected wavelet is that the similarity coefficient of the actual global term wavelet reaches its maximum after applying the global term operator.

[0120] Furthermore, if the optimal desired wavelet is not ideal enough, iterative optimization can be performed on the desired wavelet.

[0121] In addition, the optimal expected wavelet can be appropriately subjected to frequency extension processing.

[0122] S30512: Perform correlation processing on the global term wavelet and the optimal expectation wavelet to obtain the cross-correlation between the global term wavelet and the optimal expectation wavelet.

[0123] In this embodiment, correlation processing is performed on the global term wavelet and the optimal expectation wavelet to obtain the cross-correlation between the global term wavelet and the optimal expectation wavelet. The calculation formula is as follows:

[0124] In the formula, The cross-correlation between the optimal expectation wavelet and the global term wavelet; w(i) is the global term wavelet; is the optimal expected wavelet; i and j are the time point numbers determined from each time point t.

[0125] S30513: Calculate the global term operator based on the cross-correlation between the global term wavelet and the optimal expectation wavelet.

[0126] In this embodiment, the global term operator is calculated based on the cross-correlation between the global term wavelet and the optimal expectation wavelet. The calculation formula is as follows:

[0127] In the formula, f k (0) to f k (L-1) represents the elements of the global term operator; k is the identifier for each work zone; to These are the elements in the cross-correlation between the optimal expectation wavelet and the global term wavelet.

[0128] S3052: Perform an inverse Fourier transform on the power spectra of the remaining terms in all terms to obtain the corresponding autocorrelation.

[0129] In this embodiment, an inverse Fourier transform is performed on the power spectra of the remaining terms in all terms to obtain the corresponding autocorrelation of each term. The calculation formula is as follows:

[0130] In the formula, R s (t) represents the autocorrelation of the shot point term; R r (t) represents the autocorrelation of the detector term; R m (t) represents the autocorrelation of the midpoint term in the shot-receiver test; R h (t) Autocorrelation of the shot-receiver distance term; S(f) is the logarithmic power spectrum of the shot term; R(f) is the logarithmic power spectrum of the receiver term; M(f) is the logarithmic power spectrum of the shot-receiver midpoint term; H(f) is the logarithmic power spectrum of the shot-receiver distance term; e -i2πft Let be a complex exponential function; i be the imaginary unit; t be the time points, where t = t1, t2, ..., tt3. n ,t1,t n These are the start and end times of the preset time period, respectively.

[0131] S3053: Calculate the corresponding operators based on the autocorrelation of each item.

[0132] Specifically, step S3053 includes:

[0133] S30531: Calculate the shot point term operator based on the autocorrelation of the shot point term.

[0134] S30532: Calculate the detector term operator based on the autocorrelation of the detector term.

[0135] S30533: Calculate the midpoint term operator based on the autocorrelation of the midpoint term.

[0136] S30534: Calculate the operator for the shot-receiver distance term based on the autocorrelation of the shot-receiver distance term.

[0137] Specifically, step S30531 includes: calculating the shot point term operator based on the autocorrelation of the shot point term and according to the Torbritz equations; step S30532 includes: calculating the receiver point term operator based on the autocorrelation of the receiver point term and according to the Torbritz equations; step S30533 includes: calculating the shot-receiver midpoint term operator based on the autocorrelation of the shot-receiver midpoint term and according to the Torbritz equations; and step S30534 includes: calculating the shot-receiver distance term operator based on the autocorrelation of the shot-receiver distance term and according to the Torbritz equations.

[0138] In this embodiment, the calculation formula for the shot point term operator based on the autocorrelation of shot point terms is as follows:

[0139] In the formula, h S (0) to h g (L-1) represents the elements of the shot point operator; R S (τ) to R S (τ+L-1) represents the elements of the autocorrelation of the shot point term.

[0140] In summary, the seismic record calibration method provided in this embodiment obtains the seismic frequency domain by performing Fourier transforms on each target data at multiple time points; calculates the corresponding logarithmic power spectrum based on the seismic frequency domain; decomposes the logarithmic power spectrum according to a preset consistency method to obtain the response of all terms; performs exponential extraction processing on the response of all terms to obtain the corresponding power spectrum of all terms; and calculates the corresponding all-term operators based on the power spectrum of all terms. By calculating the global term operator, shot point term operator, receiver point term operator, shot-receiver midpoint term operator, and shot-receiver distance term operator, the seismic record is calibrated, thereby eliminating seismic data distortion caused by source and receiver coupling, improving seismic wavelet consistency, and enhancing the resolution of the image corresponding to the seismic record.

[0141] Figure 4 is a schematic diagram of the structure of the seismic record calibration device provided in an embodiment of this application. As shown in Figure 4, the seismic record calibration device includes: a receiving module 401, a determining module 402, a calculating module 403, a processing module 404, a calibration module 405, and an output module 406.

[0142] The receiving module 401 is used to receive seismic records from multiple work areas collected by the seismic sensor, and to extract the seismic records of all work areas according to a preset time period to obtain the seismic data of all work areas.

[0143] The determination module 402 is used to determine multiple time points within a preset time period, and to collect corresponding target data from the seismic data of all work areas based on the multiple time points;

[0144] The calculation module 403 is used to calculate the autocorrelation superimposed data of each target data based on multiple time points, and to calculate the corresponding all-item operators based on their respective autocorrelation superimposed data;

[0145] Processing module 404 is used to perform deconvolution processing on the seismic data of all work areas according to various operators, so as to obtain the deconvolutioned seismic data of all work areas.

[0146] The calibration module 405 is used to calibrate the seismic records of multiple work areas based on the seismic data after deconvolution of each work area, so as to obtain the calibrated seismic records of multiple work areas.

[0147] Output module 406 is used to output calibrated seismic records from multiple work areas.

[0148] In one possible implementation, the computing module 403 specifically includes:

[0149] Transformation unit 4031 is used to perform Fourier transform on the respective related superimposed data to obtain the seismic frequency domain;

[0150] The first calculation unit 4032 is used to calculate the corresponding logarithmic power spectrum based on the earthquake frequency domain.

[0151] Decomposition unit 4033 is used to decompose the logarithmic power spectrum according to the preset Gauss-Seidel iterative method to obtain the response of all corresponding terms;

[0152] The processing unit 4034 is used to perform exponential extraction processing on the responses of all terms to obtain the power spectrum of all terms.

[0153] The second calculation unit 4035 is used to calculate the corresponding all-term operators based on the power spectrum of all terms.

[0154] In one possible implementation, Fourier transforms are performed on the target data at multiple time points to obtain the seismic frequency domain, and the calculation formula is as follows:

[0155] X(f)=∫x(t)e -i2πft dt

[0156] In the formula, X(f) represents the seismic frequency domain; x(t) represents the seismic data; e -i2πft Let be a complex exponential function; i be the imaginary unit; t be the time points, where t = t1, t2, ..., tt3. n ,t1,t n These are the start and end times of the preset time period, respectively.

[0157] Accordingly, the formula for calculating the corresponding logarithmic power spectrum based on the seismic frequency domain is:

[0158] C(f)=ln‖X(f)‖ 2

[0159] In the formula, C(f) is the logarithmic power spectrum; X(f) is the seismic frequency domain;

[0160] Accordingly, all terms are global terms, shot point terms, receiver point terms, shot-receiver midpoint terms, and shot-receiver distance terms; the logarithmic power spectrum is decomposed according to a preset consistency method to obtain the response of all corresponding terms, and the calculation formula is as follows:

[0161] In equation (1), The response for the global item; k is the flag for each work area; Ng is the maximum number of channels for all data; Logarithmic power spectrum; For the response of the shot point item; The response of the detector point term; The response to the midpoint item of the shot inspection; is the response of the shot-receiver distance term; (n) is the result of the nth iteration; (n-1) is the result of the (n-1)th iteration; i is the shot number in the shot set; j is the receiver point number;

[0162] In equation (2), The response of the shot point term; Ns is the largest path number in the shot set i;

[0163] In equation (3), The response of the detector point term; Nr is the largest channel number in the detector point set j;

[0164] In equation (4), The response of the shot-receiver midpoint term; Nm is the largest channel number in CMP(i+j) / 2;

[0165] In equation (5), The response to the shot-receiver distance term; Nh is the maximum track number in the shot-receiver distance (ij) / 2;

[0166] Accordingly, the responses of all terms are subjected to exponential extraction processing to obtain the power spectrum of all terms, and the calculation formula is as follows:

[0167] In the formula, G k (f) represents the global power spectrum; k is the identifier for each work zone; S i (f) represents the power spectrum of the shot point term; R j (f) is the power spectrum of the detector term; M (i+j) / 2 (f) is the power spectrum of the midpoint term of the shot-receiver; H (i-j) / 2 (f) represents the power spectrum of the shot-receiver distance term.

[0168] In one possible implementation, the second computing unit 4035 specifically includes:

[0169] The first calculation unit 40351 is used to determine the global term wavelet based on the global term power spectrum and to calculate the global term operator based on the global term wavelet.

[0170] The inverse transform unit 40352 is used to perform an inverse Fourier transform on the power spectra of the remaining terms in all terms to obtain the corresponding autocorrelation terms.

[0171] The second calculation unit 40353 is used to calculate the corresponding operators based on the autocorrelation of each item.

[0172] In one possible implementation, an inverse Fourier transform is performed on the power spectra of the remaining terms in all terms to obtain the corresponding autocorrelation, and the formula for calculating the autocorrelation is as follows:

[0173] In the formula, R s (t) represents the autocorrelation of the shot point term; R r (t) represents the autocorrelation of the detector term; R m (t) represents the autocorrelation of the midpoint term in the shot-receiver test; R h (t) Autocorrelation of the shot-receiver distance term; S(f) is the logarithmic power spectrum of the shot term; R(f) is the logarithmic power spectrum of the receiver term; M(f) is the logarithmic power spectrum of the shot-receiver midpoint term; H(f) is the logarithmic power spectrum of the shot-receiver distance term; e -i2πft Let be a complex exponential function; i be the imaginary unit; t be the time points, where t = t1, t2, ..., tt3. n ,t1,t n These are the start and end times of the preset time period, respectively.

[0174] In one possible implementation, the first computing unit 40351 specifically includes:

[0175] The first processing unit 403511 is used to process the preset desired wavelet according to the preset optimization method to obtain the optimal desired wavelet.

[0176] The second processing unit 403512 is used to perform correlation processing on the global term wavelet and the optimal expectation wavelet to obtain the cross-correlation between the global term wavelet and the optimal expectation wavelet;

[0177] The computation unit 403513 is used to calculate the global term operator based on the cross-correlation between the global term wavelet and the optimal expectation wavelet.

[0178] In one possible implementation, the global term wavelet and the optimal expectation wavelet are correlated to obtain their cross-correlation, and the formula for calculating this cross-correlation is as follows:

[0179] In the formula, The cross-correlation between the optimal expectation wavelet and the global term wavelet; w(i) is the global term wavelet; The optimal expected wavelet; i and j are the time point indices determined from each time point t;

[0180] Accordingly, the global term operator is calculated based on the cross-correlation between the global term wavelet and the optimal expectation wavelet, and the calculation formula is as follows:

[0181] In the formula, f k (0) to f k (L-1) represents the elements of the global term operator; k is the identifier for each work zone; to These are the elements in the cross-correlation between the optimal expectation wavelet and the global term wavelet.

[0182] In one possible implementation, the second computing unit 40353 specifically includes:

[0183] The first calculation unit 403531 is used to calculate the shot point term operator based on the autocorrelation of the shot point term;

[0184] The first calculation unit 403532 is used to calculate the detector point term operator based on the autocorrelation of the detector point term;

[0185] The first calculation unit 403533 is used to calculate the midpoint term operator based on the autocorrelation of the midpoint term;

[0186] The first calculation unit 403534 is used to calculate the shot-receiver distance operator based on the autocorrelation of the shot-receiver distance term.

[0187] In one possible implementation, the shot point term operator is calculated based on the autocorrelation of the shot point term, and the calculation formula is as follows:

[0188] In the formula, h S (0) to h g (L-1) represents the elements of the shot point operator; R S (τ) to R S (τ+L-1) represents the elements of the autocorrelation of the shot point term.

[0189] The apparatus provided in this embodiment can be used to execute the technical solutions of the above method embodiments. Its implementation principle and technical effects are similar, and will not be described again here.

[0190] Figure 5 is a schematic diagram of the hardware structure of the server provided in this embodiment. As shown in Figure 5, the server in this embodiment includes: a processor 501 and a memory 502; the memory stores computer execution instructions; at least one processor executes the computer execution instructions stored in the memory, causing at least one processor to execute the above-described seismic record calibration method.

[0191] Alternatively, the memory 502 can be either standalone or integrated with the processor 501.

[0192] When the memory 502 is configured independently, the server also includes a bus 503 for connecting the memory 502 and the processor 501.

[0193] This application also provides a computer storage medium storing computer execution instructions. When the processor executes the computer execution instructions, the above-mentioned seismic record calibration method is implemented.

[0194] This application also provides a computer program product, including a computer program that, when executed by a processor, implements the above-described seismic record calibration method.

[0195] In the several embodiments provided in this application, it should be understood that the disclosed devices and methods can be implemented in other ways. For example, the device embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be indirect coupling or communication connection through some interfaces, devices, or modules, and may be electrical, mechanical, or other forms.

[0196] The modules described as separate components may or may not be physically separate. The components shown as modules may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to implement the solution of this embodiment according to actual needs.

[0197] Furthermore, the functional modules in the various embodiments of this application can be integrated into one processing unit, or each module can exist physically separately, or two or more modules can be integrated into one unit. The unit composed of the above modules can be implemented in hardware or in the form of hardware plus software functional units.

[0198] The integrated modules described above, implemented as software functional modules, can be stored in a computer-readable storage medium. These software functional modules, stored in a storage medium, include several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) or processor to execute some steps of the methods of the various embodiments of this application.

[0199] It should be understood that the aforementioned processor can be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), etc. A general-purpose processor can be a microprocessor or any conventional processor. The steps of the method disclosed in this invention can be directly manifested as being executed by a hardware processor, or executed by a combination of hardware and software modules within the processor.

[0200] The memory may include high-speed RAM, and may also include non-volatile storage (NVM), such as at least one disk storage device, and may also be a USB flash drive, external hard drive, read-only memory, disk or optical disc, etc.

[0201] The bus can be an Industry Standard Architecture (ISA) bus, a Peripheral Component Interconnect (PCI) bus, or an Extended Industry Standard Architecture (EISA) bus, etc. Buses can be categorized as address buses, data buses, control buses, etc. For ease of illustration, the buses shown in the accompanying drawings are not limited to a single bus or a single type of bus.

[0202] The aforementioned storage media can be implemented from any type of volatile or non-volatile storage device or a combination thereof, such as Static Random-Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read-Only Memory (EPROM), Programmable Read-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk. The storage media can be any available medium accessible to general-purpose or special-purpose computers.

[0203] An exemplary storage medium is coupled to a processor, enabling the processor to read information from and write information to the storage medium. Alternatively, the storage medium can be an integral part of the processor. Both the processor and the storage medium can reside in an Application Specific Integrated Circuit (ASIC). Alternatively, the processor and storage medium can exist as discrete components in an electronic device or host device.

[0204] Those skilled in the art will understand that all or part of the steps of the above-described method embodiments can be implemented by hardware related to program instructions. The aforementioned program can be stored in a computer-readable storage medium. When executed, the program performs the steps of the above-described method embodiments; and the aforementioned storage medium includes various media capable of storing program code, such as ROM, RAM, magnetic disks, or optical disks.

[0205] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of this application.

Claims

1. A method for calibrating seismic records, characterized in that, Applied to servers, including: It receives seismic records from multiple work areas collected by seismic sensors, and extracts seismic records from all work areas according to a preset time period to obtain seismic data for all work areas. Multiple time points within the preset time period are determined, and corresponding target data are collected from the seismic data of all work areas based on the multiple time points; Calculate the autocorrelation superposition data of each target data based on the multiple time points, and calculate all corresponding item operators based on their respective autocorrelation superposition data; Based on the various operators, the seismic data of all work areas are uniformly deconvolved to obtain the deconvolved seismic data of all work areas; The seismic records of the multiple work areas are calibrated based on the seismic data after deconvolution of each work area to obtain calibrated seismic records for multiple work areas; Output the calibrated seismic records for the multiple work areas.

2. The method according to claim 1, characterized in that, The calculation of all corresponding terms based on their respective overlay data includes: Fourier transform is performed on the respective related superimposed data to obtain the earthquake frequency domain; Calculate the corresponding logarithmic power spectrum based on the earthquake frequency domain; The logarithmic power spectrum is decomposed according to a preset Gauss-Seidel iterative method to obtain the responses of all corresponding terms; The responses to all the items are subjected to exponential extraction processing to obtain the power spectrum of all the corresponding items; Calculate all corresponding operators based on the power spectra of all terms.

3. The method according to claim 1 or 2, characterized in that, The formula for calculating the seismic frequency domain by performing Fourier transform on the target data based on the multiple time points is as follows: In the formula, X(f) represents the earthquake frequency domain; x(t) represents the earthquake data; Let be a complex exponential function; i be the imaginary unit; t be the time points, where t = t1, t2, ..., tt3. n ,t1,t n These are the start and end times of the preset time period, respectively. Accordingly, the formula for calculating the corresponding logarithmic power spectrum based on the earthquake frequency domain is: C(f)=ln‖X(f)‖ 2 In the formula, C(f) is the logarithmic power spectrum; X(f) is the seismic frequency domain; Accordingly, all of the terms are global terms, shot point terms, receiver point terms, shot-receiver midpoint terms, and shot-receiver distance terms; the calculation formula for decomposing the logarithmic power spectrum according to a preset consistency method to obtain the response of all corresponding terms is: In equation (1), The response to the global item; k is the flag for each work zone; Ng is the maximum number of channels for all data; Logarithmic power spectrum; response for the shot point term; response of the detection point term; response for the shot-inside-point term; is the response of the shot-receiver distance term; (n) is the result of the nth iteration; (n-1) is the result of the (n-1)th iteration; i is the shot number in the shot set; j is the receiver point number; in formula (2), The response of the shot point term; Ns is the largest path number in the shot set i; in formula (3), The response of the detector point term; Nr is the largest channel number in the detector point set j; in formula (4), The response of the shot-receiver midpoint term; Nm is the largest channel number in CMP(i+j) / 2; in formula (5), The response to the shot-receiver distance term; Nh is the maximum track number in the shot-receiver distance (ij) / 2; Accordingly, the response of the all terms is subjected to an exponential extraction process to obtain a corresponding all term power spectrum, and the calculation formula is: In the formula, G k (f) represents the global power spectrum; k is the identifier for each work zone; S i (f) represents the power spectrum of the shot point term; R j (f) is the power spectrum of the detector term; M (i+j) / 2 (f) is the power spectrum of the midpoint term of the shot-receiver; H (i-j) / 2 (f) represents the power spectrum of the shot-receiver distance term.

4. The method according to any one of claims 1 to 3, characterized in that, The step of calculating the corresponding operators based on the power spectrum of all terms includes: The global term wavelet is determined based on the global term power spectrum, and the global term operator is calculated based on the global term wavelet. Perform an inverse Fourier transform on the power spectra of the remaining terms in all the power spectra to obtain the corresponding autocorrelation terms; The corresponding operators are calculated based on the autocorrelation of each item.

5. The method according to any one of claims 1 to 4, characterized in that, The formula for calculating the autocorrelation of the remaining power spectra in all the power spectra is as follows: In the formula, R s (t) represents the autocorrelation of the shot point term; R r (t) represents the autocorrelation of the detector term; R m (t) represents the autocorrelation of the midpoint term in the shot-receiver test; R h (t) Autocorrelation of the shot-receiver distance term; S(f) is the logarithmic power spectrum of the shot point term; R(f) is the logarithmic power spectrum of the receiver point term; M(f) is the logarithmic power spectrum of the shot-receiver midpoint term; H(f) is the logarithmic power spectrum of the shot-receiver distance term; Let i be a complex exponential function; i be the imaginary unit; t be the time points, where t = t1, t2, ..., tt3. n ,t1,t n These are the start and end times of the preset time period, respectively.

6. The method according to any one of claims 1 to 4, characterized in that, The step of calculating the global term operator based on the global term wavelet includes: The preset desired wavelet is processed according to the preset optimization method to obtain the optimal desired wavelet; The global term wavelet and the optimal expectation wavelet are subjected to correlation processing to obtain the cross-correlation between the global term wavelet and the optimal expectation wavelet; The global term operator is calculated based on the cross-correlation between the global term wavelet and the optimal expectation wavelet.

7. The method according to claim 6, characterized in that, The formula for calculating the cross-correlation between the global term wavelet and the optimal expectation wavelet, obtained by performing correlation processing on the global term wavelet and the optimal expectation wavelet, is as follows: In the formula, the The cross-correlation between the optimal expected wavelet and the global term wavelet; w(i) is the global term wavelet; The optimal expected wavelet; i and j are the time point indices determined from each time point t; Accordingly, the calculation formula for the global term operator based on the cross-correlation of the global term wavelet and the optimal expectation wavelet is as follows: In the formula, f k (0) to f k (L-1) represents the elements of the global term operator; k is the identifier for each work zone; to These are the elements in the cross-correlation between the optimal expected wavelet and the global term wavelet.

8. The method according to any one of claims 5 to 7, characterized in that, The step of calculating the corresponding operators based on the autocorrelation of each term includes: Calculate the shot point term operator based on the autocorrelation of the shot point term; calculate the receiver term operator based on the autocorrelation of the receiver term; calculate the shot-receiver midpoint term operator based on the autocorrelation of the shot-receiver midpoint term; calculate the shot-receiver distance term operator based on the autocorrelation of the shot-receiver distance term.

9. The method according to any one of claims 5 to 8, characterized in that, The calculation formula for the operator that calculates the shot point term based on the autocorrelation of the shot point term is as follows: In the formula, h S (0) to h g (L-1) represents the elements of the shot point operator; R S (τ) to R S (τ+L-1) represents the elements of the autocorrelation of the shot point term.

10. A seismic record calibration device, characterized in that, Applied to servers, including: The receiving module is used to receive seismic records from multiple work areas collected by seismic sensors, and to extract the seismic records of all work areas according to a preset time period to obtain the seismic data of all work areas. The determination module is used to determine multiple time points within the preset time period, and to collect corresponding target data from the seismic data of all work areas based on the multiple time points; The calculation module is used to calculate the autocorrelation superposition data of each target data according to the multiple time points, and to calculate the corresponding all-item operators according to their respective autocorrelation superposition data; The processing module is used to perform deconvolution processing on the seismic data of all work areas according to various operators, so as to obtain the deconvolutioned seismic data of all work areas. The calibration module is used to calibrate the seismic records of the multiple work areas based on the seismic data after deconvolution of each work area, so as to obtain the calibrated seismic records of the multiple work areas. The output module is used to output the calibrated seismic records of the multiple work areas.

11. A server, characterized in that, include: At least one processor and memory; The memory stores computer-executed instructions; The at least one processor executes computer execution instructions stored in the memory, causing the at least one processor to perform the seismic record calibration method as described in any one of claims 1 to 9.

12. A computer storage medium, characterized in that, The computer storage medium stores computer execution instructions, and when the processor executes the computer execution instructions, it implements the seismic record calibration method as described in any one of claims 1 to 9.

13. A computer program product, characterized in that, It includes a computer program that, when executed by a processor, implements the seismic record calibration method according to any one of claims 1 to 9.

Images