A least squares integral method prestack depth migration method and device

Abstract

The application discloses a least square integral method pre-stack depth migration method and device. The method comprises the following steps: the application performs migration calculation, reverse migration imaging and the modulus square of the migration imaging result and the reverse migration imaging result by iteration, calculates a first step length, updates the residual error and the migration imaging result by using the first step length, and uses the cross-correlation coefficient of the residual error and the residual error update direction when updating the residual error; if the cross-correlation coefficient is less than 0, the residual error and the residual error update direction of the channel are both set to 0; then the migration calculation is performed by using the updated residual error, the conjugate gradient direction is calculated, the reverse migration calculation is repeatedly performed again by using the conjugate gradient direction, and the preset convergence standard is reached. The application effectively improves the imaging iteration speed of the model in the middle and deep layer position, and improves the imaging quality of the least square migration result by suppressing the residual error data with the negative cross-correlation coefficient of the update direction in the iteration process.

Description

Technical Field

[0001] This invention relates to the field of petroleum geophysical exploration data processing technology, and in particular to a method and apparatus for pre-stack depth migration using the least squares integral method. Background Technology

[0002] Conventional pre-stack depth migration is affected by factors such as the observation system and missing data, which often leads to poor imaging quality, such as uneven illumination and low resolution. In order to solve these problems and further improve the quality of migration imaging, the academic community has proposed the least squares migration method, which can effectively improve the imaging quality.

[0003] In the application of real-world 3D data, the pre-stack depth migration technique using the least squares integral method faces challenges such as managing multiple sets of data (travel time tables, migration results, model gradients, simulation data, residual update directions, etc.) within limited memory, parallel implementation of inverse migration, slow convergence of 3D data inversion, and significant noise issues in the inversion results. Summary of the Invention

[0004] In view of the above problems, the present invention is proposed to provide a least-squares integral pre-stack depth migration method and apparatus that overcomes or at least partially solves the above problems.

[0005] In a first aspect, embodiments of the present invention provide a pre-stack depth migration method using least squares integration, comprising:

[0006] Migration calculations are performed on the acquired seismic data to obtain migration imaging results;

[0007] Using the offset imaging result as input, reverse offset calculation is performed to obtain the reverse offset imaging result;

[0008] The first step length is calculated based on the square of the modulus of the offset imaging result and the square of the modulus of the reverse offset imaging result.

[0009] Based on the first step length, update the residual of this iteration and the offset imaging result obtained in this iteration; when updating the residual, calculate the cross-correlation coefficient between each residual and the residual update direction. If the cross-correlation coefficient is less than 0, set both the residual and the residual update direction of that trace to 0; the residual update direction is the reverse offset result of that trace.

[0010] Determine whether the current iteration meets the preset convergence criteria;

[0011] If the conditions are met, the final offset imaging result is generated based on the offset imaging result obtained in this iteration;

[0012] If the condition is not met, the offset calculation is performed again using the residual updated in this iteration to obtain the offset imaging result; and the conjugate gradient direction is determined based on the offset imaging result obtained in the current calculation and the offset imaging result obtained in the previous calculation.

[0013] Return to the conjugate gradient direction and repeat the inverse offset calculation step, incrementing the iteration count by 1, until the preset convergence criterion is reached.

[0014] In one embodiment, migration calculations are performed on the acquired seismic data to obtain migration imaging results, including:

[0015] Based on memory size and offset parameters, calculate part or all of the offset group data and part or all of the imaging space of the offset group in the seismic data that a single offset unit task needs to process;

[0016] Calculate the travel time table based on the offset velocity field;

[0017] Multiple migration unit tasks are executed through the job node, the seismic data corresponding to the task and the travel time data corresponding to the travel time table are read, the migration imaging is completed and saved;

[0018] Imaging results from the same offset group are merged to obtain offset imaging results from different offset groups.

[0019] In one embodiment, the offset calculation process further includes: preconditioning, wherein the preconditioning process involves adding weight coefficients to the offset result. Where a ranges from 0.5 to 2, b ranges from 0.00025 to 0.01, and offset represents the middle offset length of the offset group. The pre-conditioning imaging formula is as follows:

[0020]

[0021] in, This represents the spatial imaging position; in two dimensions, it is a two-dimensional vector, and in three dimensions, it is a three-dimensional vector, where the depth dimension is represented by z. Indicates spatial imaging position Offset grouping The imaging results The weighting coefficients mentioned above are: , Indicates the weighting coefficient. Indicates gun and gun inspection configuration , The spectrum of earthquake data, Indicates offset grouping, In two dimensions, it is a one-dimensional vector; in three dimensions, it is a two-dimensional vector. This represents the filter factor; in the two-dimensional case, it is a semi-derivative filter, i.e.: In the three-dimensional case, it is derivative filtering, that is: , Represents the spectrum of the seismic wavelet. , These represent the distances from the shot point and receiver point to the imaging point, respectively. When they left, In the two-dimensional case, it represents a single integral, and in the three-dimensional case, it represents a double integral.

[0022] In one embodiment, the offset imaging result is used as input to perform inverse offset calculation to obtain the inverse offset imaging result, including:

[0023] The task node executes multiple inverse offset unit tasks, reads the results of the offset imaging, the travel time data in the corresponding travel time table, and the track head information of the data in the input offset distance group, performs the corresponding inverse offset calculation and saves it.

[0024] The inverse offset imaging results of the same offset group are merged to obtain inverse offset imaging results of different offset groups.

[0025] In one embodiment, it also includes:

[0026] In the inverse migration calculation, the imaging results are pre-conditioned along the depth direction. The pre-conditioned inverse migration formula is as follows:

[0027]

[0028] in, Indicates gun and gun inspection configuration , Inverse offset synthetic data, Indicates offset grouping, In two dimensions, it is a one-dimensional vector; in three dimensions, it is a two-dimensional vector. This represents the filter factor; in the two-dimensional case, it is a semi-derivative filter, i.e.: In the three-dimensional case, it is derivative filtering, that is: , Represents the spectrum of the seismic wavelet. This represents the spatial imaging position; in two dimensions, it is a two-dimensional vector, and in three dimensions, it is a three-dimensional vector, where the depth dimension is represented by z. Indicates the weighting coefficient. Indicates spatial imaging position Offset grouping The imaging results , These represent the distances from the shot point and receiver point to the imaging point, respectively. When they left, The weighting coefficients mentioned above are: , In the two-dimensional case, it represents a double integral, and in the three-dimensional case, it represents a triple integral.

[0029] In one embodiment, the first step length is obtained by dividing the square of the modulus of the offset imaging result by the inverse offset imaging result.

[0030] In one embodiment, updating the residual of the current iteration and the offset imaging result obtained in the current iteration based on the first step length includes:

[0031] When the cross-correlation coefficient between the input data at the same location and the reverse migration imaging result is greater than or equal to zero, the input data of that channel is subtracted from the product of the first step length and the reverse migration imaging result of this iteration to obtain the residuals of multiple offset groups in this iteration. When the cross-correlation coefficient is less than zero, the residual of that channel is set to 0.

[0032] The offset imaging result of the previous iteration is obtained by adding the first step length to the offset imaging result of the current iteration;

[0033] The input data is seismic data in the first iteration and residuals obtained from the previous iteration in subsequent iterations.

[0034] The offset imaging result of the previous iteration is data 0 in the first iteration, and the offset imaging result of the previous iteration in subsequent iterations.

[0035] In one embodiment, if a preset convergence criterion is met, a final migration imaging result is generated based on the migration imaging result obtained in this iteration, including:

[0036] If the preset convergence criteria are met, the migration imaging result obtained in this iteration is multiplied by the precondition factor to generate the final migration imaging result.

[0037] In one embodiment, if the condition is not met, the offset calculation is performed again using the residual updated in this iteration to obtain the offset imaging result; and the conjugate gradient direction is determined based on the currently calculated offset imaging result and the offset imaging result obtained in the previous offset calculation, including:

[0038] The second step length is obtained by dividing the squared modulus of the offset imaging result obtained from the previous offset calculation by the squared modulus of the offset imaging result obtained from the previous offset calculation.

[0039] Based on the second step size, the offset imaging result obtained from this offset calculation, and the offset imaging result obtained from the previous offset calculation, the conjugate gradient direction is obtained.

[0040] Secondly, embodiments of the present invention provide a pre-stack depth migration device using the least squares integral method, comprising: a migration calculation module, an inverse migration calculation module, an update module, a judgment module, a result output module, and a conjugate calculation module;

[0041] The migration calculation module is used to perform migration calculations on the acquired seismic data to obtain migration imaging results; when the judgment module determines that the current iteration does not meet the preset convergence criteria, it uses the residual updated after the current iteration to perform migration calculations again to obtain migration imaging results.

[0042] The reverse offset calculation module is used to take the offset imaging result as input, perform reverse offset calculation, and obtain the reverse offset imaging result;

[0043] The update module is used to calculate the first step length based on the square of the modulus of the offset imaging result and the square of the modulus of the reverse offset imaging result; update the residual of this iteration and the offset imaging result obtained in this iteration based on the first step length; when updating the residual, calculate the cross-correlation coefficient between the residual and the residual update direction for each trace; if the cross-correlation coefficient is less than 0, set both the residual and the residual update direction of that trace to 0; the residual update direction is the reverse offset result of that trace.

[0044] The judgment module is used to determine whether the current iteration meets the preset convergence criteria;

[0045] The result output module is used to generate the final offset imaging result based on the offset imaging result obtained in the current iteration when the judgment module determines that the current iteration meets the preset convergence criterion.

[0046] The conjugate calculation module is used to determine the conjugate gradient direction based on the offset imaging result currently calculated by the offset calculation module and the offset imaging result obtained in the previous offset calculation when the judgment module determines that the current iteration does not meet the preset convergence criterion.

[0047] Thirdly, embodiments of the present invention provide a computing device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the pre-stack depth migration method of least squares integration as described above.

[0048] Fourthly, embodiments of the present invention provide a computer-readable storage medium storing a computer program that, when executed by a processor, implements the aforementioned least-squares integral pre-stack depth migration method.

[0049] The beneficial effects of the above-described technical solutions provided in the embodiments of the present invention include at least the following:

[0050] This invention provides a method and apparatus for pre-stack depth migration using the least squares integral method. It iteratively calculates the migration and back-migration imaging results, and uses the squared modulus of the migration and back-migration imaging results to calculate the first-step length. Using this first-step length, the residuals and migration imaging results are updated. When updating the residuals, the cross-correlation coefficient between the residuals and the residual update direction is used. If the cross-correlation coefficient is less than 0, both the residuals and the residual update direction for that trace are set to 0. The process then returns to the conjugate gradient direction and repeats the back-migration calculation until a preset convergence criterion is reached. This invention effectively improves the imaging iteration speed of the model in mid-to-deep layers. Furthermore, by suppressing residual data with negative cross-correlation coefficients with the update direction during iteration, it improves the signal-to-noise ratio of the inversion results, overcoming the slow inversion convergence and severe noise problems present in existing technologies.

[0051] Furthermore, by dividing the data into multiple unit tasks during the offset calculation process and implementing the least squares offset process using a parallel computing framework, computational efficiency can be improved while effectively enhancing the imaging effect of actual data.

[0052] Other features and advantages of the invention will be set forth in the following description, and will be apparent in part from the description, or may be learned by practicing the invention. The objects and other advantages of the invention may be realized and obtained by means of the structures particularly pointed out in the written description and the accompanying drawings.

[0053] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0054] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings:

[0055] Figure 1 This is a flowchart of the pre-stack depth migration method using the least squares integral method in an embodiment of the present invention;

[0056] Figure 2A-2B This is a schematic diagram comparing the least squares offset results before and after preconditioning in an embodiment of the present invention;

[0057] Figure 3A This is an example diagram showing the cross-correlation results between residuals and residual update directions in an embodiment of the present invention;

[0058] Figure 3B This is an example diagram showing the energy of each channel of the residual data in an embodiment of the present invention;

[0059] Figure 3CThis is an example diagram showing the energy of each path in the residual update direction in an embodiment of the present invention;

[0060] Figure 3D This is an example diagram showing the maximum energy of each channel of the residual data in an embodiment of the present invention;

[0061] Figure 3E This is an example diagram showing the maximum energy of each path in the residual update direction in an embodiment of the present invention;

[0062] Figure 4A and 4B This is an example diagram comparing the least squares offset results before and after residual data processing based on the cross-correlation results between the residuals and the residual update direction in an embodiment of the present invention.

[0063] Figure 5 This is a structural block diagram of the pre-stack depth migration device using the least squares integral method in an embodiment of the present invention. Detailed Implementation

[0064] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0065] The following detailed description, with reference to the accompanying drawings, illustrates the specific implementation of the least squares integral pre-stack depth migration method and related apparatus provided in the embodiments of the present invention.

[0066] This invention provides a pre-stack depth migration method using the least squares integral method, referring to... Figure 1 As shown, it includes the following steps:

[0067] S11. Perform migration calculations on the acquired seismic data to obtain the migration imaging results;

[0068] S12. Using the offset imaging result as input, perform reverse offset calculation to obtain the reverse offset imaging result;

[0069] S13. Calculate the first step length based on the square of the modulus of the offset imaging result and the square of the modulus of the reverse offset imaging result.

[0070] S14. Based on the first step length, update the residual of this iteration and the offset imaging result obtained in this iteration; when updating the residual, calculate the cross-correlation coefficient between each residual and the residual update direction. If the cross-correlation coefficient is less than 0, set both the residual and the residual update direction of that trace to 0; the residual update direction is the reverse offset result of that trace.

[0071] S15. Determine whether the current iteration meets the preset convergence criterion; if it does, proceed to step S16; if it does not, proceed to step S17.

[0072] S16. Based on the migration imaging results obtained in this iteration, generate the final migration imaging results;

[0073] S17. Using the residual updated in this iteration, perform the offset calculation again to obtain the offset imaging result;

[0074] S18. Based on the offset imaging results obtained from the current calculation and the offset imaging results obtained from the previous offset calculation, determine the conjugate gradient direction; then return to execute S12, increment the iteration count by 1, until the preset convergence criterion is reached.

[0075] Furthermore, step S11 above performs migration calculations on the acquired seismic data to obtain migration imaging results, which can be achieved, for example, through the following steps:

[0076] Based on memory size and offset parameters, calculate part or all of the offset group data and part or all of the imaging space of the offset group in the seismic data that a single offset unit task needs to process;

[0077] Calculate the travel time table based on the offset velocity field;

[0078] Multiple migration unit tasks are executed through the job node, the seismic data corresponding to the task and the travel time data corresponding to the travel time table are read, the migration imaging is completed and saved;

[0079] Imaging results from the same offset group are merged to obtain offset imaging results from different offset groups.

[0080] Furthermore, in step S12 above, the offset imaging result is used as input to perform inverse offset calculation to obtain the reverse offset imaging result. This can be implemented using the following process:

[0081] The task node executes multiple inverse offset unit tasks, reads the results of the offset imaging, the travel time data in the corresponding travel time table, and the track head information of the data in the input offset distance group, performs the corresponding inverse offset calculation and saves it.

[0082] The inverse offset imaging results of the same offset group are merged to obtain inverse offset imaging results of different offset groups.

[0083] Furthermore, in step S14 above, the residual of this iteration and the offset imaging result obtained in this iteration are updated according to the length of the first step. This can be achieved, for example, through the following process:

[0084] When the cross-correlation coefficient between the input data at the same location and the reverse migration imaging result is greater than or equal to zero, the input data of that channel is subtracted from the product of the first step length and the reverse migration imaging result of this iteration to obtain the residuals of multiple offset groups in this iteration. When the cross-correlation coefficient is less than zero, the residual of that channel is set to 0.

[0085] The offset imaging result of the previous iteration is obtained by adding the first step length to the offset imaging result of the current iteration;

[0086] The input data is seismic data in the first iteration and residuals obtained from the previous iteration in subsequent iterations.

[0087] The offset imaging result of the previous iteration is data 0 in the first iteration, and the offset imaging result of the previous iteration in subsequent iterations.

[0088] Furthermore, in step S18 above, the conjugate gradient direction can be calculated in the following manner:

[0089] The second step length is obtained by dividing the squared modulus of the offset imaging result obtained from the previous offset calculation by the squared modulus of the offset imaging result obtained from the previous offset calculation.

[0090] Based on the second step size, the offset imaging result obtained from this offset calculation, and the offset imaging result obtained from the previous offset calculation, the conjugate gradient direction is obtained.

[0091] To better illustrate the pre-stack depth migration method using the least squares integral method provided in this embodiment of the invention, the steps described above will be explained in detail below with reference to a specific embodiment:

[0092] The specific implementation steps of this embodiment are as follows:

[0093] 1) Acquire seismic data and migration velocity field.

[0094] 2) Calculate the size of the unit task based on the memory size and offset parameters; each unit task corresponds to a part or all of the data of a certain offset group, or a part or all of the imaging space of a certain offset group.

[0095] This step is used to divide the size of the unit task, divide the data and imaging space of the offset group, and thus divide the range of the corresponding travel time table.

[0096] The action nodes are used to execute unit tasks separately. The execution process of each unit task can be independent of each other. Dividing the data to be calculated into several unit tasks can avoid the problem of excessive memory resource consumption caused by excessive computation.

[0097] 3) Use the offset velocity field to calculate the travel time table.

[0098] Each shot point has a series of receiver points, and the receiver points have three-dimensional spatial data. Based on the velocity data of the offset velocity field, the travel time can be calculated to form a travel time table.

[0099] 4) The master node distributes tasks. After the job node obtains each migration unit task, it reads the seismic data and travel time table corresponding to the task. After reading, it retrieves the next unit task from the master node. When the job node executes the migration unit task, multiple migration processes can complete the migration imaging. After processing, it continues to retrieve the next unit task from the master node, and so on in a loop.

[0100] 5) During the offset calculation process, a pre-conditioning process is added, that is, weighting coefficients are added to the offset results. Where a ranges from 0.5 to 2, b ranges from 0.00025 to 0.01, and offset represents the middle offset length of the offset group. The pre-conditioning imaging formula is as follows:

[0101]

[0102] in, This represents the spatial imaging position; in two dimensions, it is a two-dimensional vector, and in three dimensions, it is a three-dimensional vector, where the depth dimension is represented by z. Indicates spatial imaging position Offset grouping The imaging results The weighting coefficients mentioned above are: , Indicates the weighting coefficient. Indicates gun and gun inspection configuration , The spectrum of earthquake data, Indicates offset grouping, In two dimensions, it is a one-dimensional vector; in three dimensions, it is a two-dimensional vector. This represents the filter factor; in the two-dimensional case, it is a semi-derivative filter, i.e.: In the three-dimensional case, it is derivative filtering, that is: , Represents the spectrum of the seismic wavelet. , These represent the distances from the shot point and receiver point to the imaging point, respectively. When they left, In the two-dimensional case, it represents a single integral, and in the three-dimensional case, it represents a double integral.

[0103] The preconditioning process described above in the offset calculation can improve the imaging iteration speed of the model in mid-to-deep layers.

[0104] Figure 2A-2BThe image shows a comparison of the least squares shift results before and after preconditioning. Figure 2A The result is the least squares offset before preconditioning. Figure 2B This is the least squares offset result after preconditioning.

[0105] 6) After the offset unit task is completed, the imaging results are saved to a disk file, and the unit task completion status and result file information are returned to the master node. The master node saves this information into memory.

[0106] 7) After all the offset tasks are completed, the master node merges the imaging results of the offset distance groups according to the completed unit task information stored in memory, and obtains the final offset imaging results of different offset distance groups.

[0107] Simultaneously, the imaging results are placed in the designated directory and named according to the specified rules, serving as the model input data for the next stage of inverse migration calculation, i.e., the gradient direction of this iteration (the migration imaging results are the gradient direction); at the same time, the square of the modulus of the migration results of each migration distance group (i.e., the two-modulus data, or the 2-norm data) is calculated, and the results corresponding to all migration distance groups are saved to a file.

[0108] The process of calculating the square of the modulus of the offset result for each offset group is as follows:

[0109] Since each offset group contains data from several channels, and each channel contains data from several sampling time points, for each channel, the sum of squares of the data from each sampling time point in that channel is first calculated to obtain the value of the sum of squares for each channel. Then, the data from each channel are summed again to obtain the square of the modulus of the offset result of each offset group (i.e., L2 modulus data, or L2 norm data).

[0110] 8) After the above stages are completed, the unit tasks for inverse offset are distributed. After the job node obtains each inverse offset unit task, it reads the model input data corresponding to the task (the input data is the offset imaging result of the offset distance grouping of the above offset, which is the output of the offset step), the corresponding pre-calculated travel time table, and the track head information of the data in the input offset distance group. After reading, it takes the next unit task from the master node, and the inverse offset process completes the inverse offset calculation.

[0111] 9) In the inverse migration calculation, the imaging results are pre-conditioned along the depth direction. The pre-conditioned inverse migration formula is as follows:

[0112]

[0113] in, Indicates gun and gun inspection configuration , Inverse offset synthetic data, Indicates offset grouping, In two dimensions, it is a one-dimensional vector; in three dimensions, it is a two-dimensional vector. This represents the filter factor; in the two-dimensional case, it is a semi-derivative filter, i.e.: In the three-dimensional case, it is derivative filtering, that is: , Represents the spectrum of the seismic wavelet. This represents the spatial imaging position; in two dimensions, it is a two-dimensional vector, and in three dimensions, it is a three-dimensional vector, where the depth dimension is represented by z. Indicates the weighting coefficient. Indicates spatial imaging position Offset grouping The imaging results , These represent the distances from the shot point and receiver point to the imaging point, respectively. When they left, The weighting coefficients mentioned above are: , In the two-dimensional case, it represents a double integral, and in the three-dimensional case, it represents a triple integral.

[0114] Similar to step 5, the preconditioning process in the above-mentioned inverse migration calculation can also improve the imaging iteration speed of the model in the middle and deep layers.

[0115] 10) After the inverse offset unit task is completed, if the output position of the current inverse offset unit task is different from that of the next inverse offset unit task, or if the next task is empty, the output result of the current inverse offset unit task is written to the disk file, and the inverse offset unit task completion information and inverse offset imaging result file information are returned to the master node. The master node saves this information in memory and then continues to the next inverse offset unit task. The job node continues to execute the next inverse offset unit task until all inverse offset unit tasks are completed.

[0116] 11) After all reverse offset tasks are completed, the reverse offset results are merged. That is, the master node sums the results of the same output position completed by different job nodes according to the recorded unit task completion status and result file information. This process is because different nodes may be processing the same offset group. Here, it is necessary to merge the same offset group to obtain the reverse offset imaging result of each offset group. The result file is named according to a certain naming rule and saved to the specified folder. Then, the previous temporary files are deleted to obtain the reverse offset imaging result of each offset group.

[0117] Based on the inverse migration imaging results of each offset group, the square of the modulus of the inverse migration result data for each offset group can be further calculated (i.e., binary modulus data, L2 norm data):

[0118] Squaring the sampled values ​​at each time point of each channel within each offset group, then summing the squared values ​​of all sampled time points for the same channel to obtain the data corresponding to that channel, and finally summing the data of all channels to obtain the squared modulus corresponding to each offset group, and saving the results to a file. The calculation process of the squared modulus is similar to the implementation process of step 7 above, and will not be repeated here.

[0119] 12) Calculate the first step length based on the offset imaging result obtained in step 7 (i.e., the gradient direction of the objective function) and the square of the magnitude of the inverse offset imaging result obtained in step 11.

[0120] The length of the first step is equal to the square of the modulus of the offset result divided by the square of the modulus of the reverse offset result.

[0121] 13) Update the residual and offset results based on the length of the first step.

[0122] The residual is calculated using the following formula:

[0123] Input data - (first step length multiplied by inverse offset result) = residuals of multiple offset groups in this iteration;

[0124] Update the offset model results using the following formula:

[0125] The offset model result of the previous iteration + (the length of the first step multiplied by the offset result of this iteration) = the offset model result of this iteration;

[0126] The offset model result from the previous iteration is 0 when this step is executed for the first time. The offset model result obtained after this execution can be saved as the "offset model result from the previous iteration" for the next calculation.

[0127] The input data mentioned above refers to the actual earthquake data input during the first calculation; during the second iteration, the input data is the residual obtained from the previous iteration.

[0128] The residual data is the data after each iteration using the first step length and the inverse offset result.

[0129] During the residual update process, the cross-correlation coefficient between each residual and the residual update direction (i.e., the inverse migration result in step 11 above) is first calculated. When the cross-correlation coefficient is less than 0, it indicates that the two are inversely correlated. In this case, if the data of this trace is updated, it will have the opposite effect on improving the imaging effect of the pre-stack depth migration result. Therefore, the residual and residual update direction of this trace are both set to zero and do not participate in the iterative update.

[0130] The imaging result after this update is obtained by multiplying the migration imaging result updated with residuals by a preconditioning factor. At the same time, the residuals, the cross-correlation coefficients between the residual update directions, their respective maximum energies, and the overall energy of the gather are output to facilitate quality control and analysis of the data during the iteration process.

[0131] An example of the cross-correlation result between a residual and its update direction can be found in [reference]. Figure 3A As shown. An example plot of the energy for each channel of the residual data can be found in [reference needed]. Figure 3B As shown, an example of the energy for each path in the residual update direction can be found in [reference needed]. Figure 3C As shown, an example of the maximum energy per channel of the residual data can be found in [reference needed]. Figure 3D As shown, the maximum energy for each path in the residual update direction can be referenced. Figure 3E As shown.

[0132] If the convergence criterion for iteration is met at this point, meaning the iteration can end, then the imaging result after this update iteration can be used as the final output.

[0133] If the convergence criterion for iteration is met at this point, meaning the iteration has not ended, the imaging result after this update iteration can only be considered an intermediate result.

[0134] The convergence criteria can be preset, such as the number of iterations reaching a preset number, or the objective function reaching a preset convergence condition.

[0135] The objective function of this invention embodiment can be, for example:

[0136]

[0137] Where L represents the forward modeling operator, m represents the reflection coefficient model, d represents the acquired seismic data, x, y, z represent the three-dimensional spatial coordinates, and h represents the offset grouping parameter. Represents the regularization parameter. The L2 norm represents the sum of squares of the vector elements. This represents the first norm, which is the sum of the absolute values ​​of the vector elements. It can also be the second norm.

[0138] Figure 4A and Figure 4B To compare the least squares offset results before and after residual data processing based on the cross-correlation results of the residuals and residual update directions, the following is performed: Figure 4A For the residual data before processing, Figure 4B This is the processed residual data.

[0139] 14) Read in the updated residual, and then perform the offset calculation in steps (4) to (7) to obtain the offset imaging result, i.e. the new gradient direction;

[0140] 15) Calculate the second step size based on the squared magnitude of the gradient directions from the two steps (step 7 and step 14). Then, based on the second step size and the two gradient directions, obtain the conjugate gradient direction, which will be used as the model input data for the next stage of inverse offset calculation.

[0141] Given the limited memory, only a specified size is read at a time during the data reading process, and this process is repeated iteratively.

[0142] The formula for calculating the second step length is as follows:

[0143] Second step length = squared magnitude of current gradient direction divided by squared magnitude of previous gradient direction;

[0144] 16) Continue the iterative process described in steps (8) to (15) to finally complete the iterative solution.

[0145] Based on the same inventive concept, this embodiment of the invention also provides a pre-stack depth migration device using the least squares integral method. Since the principle of solving the problem by these devices is similar to that of the aforementioned pre-stack depth migration method using the least squares integral method, the implementation of this device can refer to the implementation of the aforementioned method, and the repeated parts will not be described again.

[0146] This invention provides a pre-stack depth migration device using the least squares integral method, with reference to... Figure 5 As shown, it includes: offset calculation module 51, inverse offset calculation module 52, update module 53, judgment module 54, result output module 55, and conjugate calculation module 56;

[0147] The migration calculation module 51 is used to perform migration calculation on the acquired seismic data to obtain the migration imaging result; when the judgment module determines that the current iteration does not meet the preset convergence criterion, it uses the residual updated after the current iteration to perform migration calculation again to obtain the migration imaging result.

[0148] The inverse offset calculation module 52 is used to take the offset imaging result as input, perform inverse offset calculation, and obtain the inverse offset imaging result.

[0149] The update module 53 is used to calculate the first step length based on the square of the modulus of the offset imaging result and the square of the modulus of the reverse offset imaging result; update the residual of this iteration and the offset imaging result obtained in this iteration based on the first step length; when calculating the residual, calculate the cross-correlation coefficient between the residual and the residual update direction for each trace; if the cross-correlation coefficient is less than 0, set both the residual and the residual update direction of that trace to 0; the residual update direction is the reverse offset result of that trace.

[0150] The judgment module 54 is used to determine whether the current iteration meets the preset convergence criteria;

[0151] The result output module 55 is used to generate the final offset imaging result based on the offset imaging result obtained in the current iteration when the judgment module determines that the current iteration meets the preset convergence criterion.

[0152] The conjugate calculation module 56 is used to determine the conjugate gradient direction based on the offset imaging result currently calculated by the offset calculation module and the offset imaging result obtained in the previous offset calculation when the judgment module determines that the current iteration does not meet the preset convergence criterion.

[0153] This invention also provides a computing device, comprising:

[0154] processor;

[0155] Memory used to store processor-executable commands;

[0156] The processor is configured to implement the pre-stack depth migration method using least squares integration as described above when executing the program.

[0157] This invention also provides a non-transitory computer-readable storage medium, which, when the instructions in the storage medium are executed by the processor of a mobile terminal, enables the mobile terminal to execute the above-described least squares integral pre-stack depth offset method.

[0158] Regarding the least squares integral pre-stack depth migration device in the above embodiments, the specific methods by which each module performs its operations have been described in detail in the embodiments related to the method, and will not be elaborated here.

[0159] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage and optical storage) containing computer-usable program code.

[0160] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1A device that provides the functions specified in one or more boxes.

[0161] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0162] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0163] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A pre-stack depth migration method using least squares integration, characterized in that, include: Migration calculations are performed on the acquired seismic data to obtain migration imaging results; Using the offset imaging result as input, reverse offset calculation is performed to obtain the reverse offset imaging result; The first step length is calculated based on the square of the modulus of the offset imaging result and the square of the modulus of the reverse offset imaging result. Based on the first step length, update the residual of this iteration and the offset imaging result obtained in this iteration; when updating the residual, calculate the cross-correlation coefficient between each residual and the residual update direction. If the cross-correlation coefficient is less than 0, set both the residual and the residual update direction of that trace to 0; the residual update direction is the reverse offset result of that trace. Determine whether the current iteration meets the preset convergence criteria; If the conditions are met, the final offset imaging result is generated based on the offset imaging result obtained in this iteration; If the conditions are not met, the offset calculation is performed again using the residual updated in this iteration to obtain the offset imaging result; The conjugate gradient direction is determined based on the current offset imaging result and the offset imaging result obtained from the previous offset calculation. Return to the conjugate gradient direction and repeat the inverse offset calculation step, incrementing the iteration count by 1, until the preset convergence criterion is reached; The step of updating the residual of the current iteration and the offset imaging result obtained in the current iteration based on the first step length includes: When the cross-correlation coefficient between the input data at the same location and the reverse migration imaging result is greater than or equal to zero, the input data of that channel is subtracted from the product of the first step length and the reverse migration imaging result of this iteration to obtain the residuals of multiple offset groups in this iteration. When the cross-correlation coefficient is less than zero, the residual of that channel is set to 0. The offset imaging result of the previous iteration is obtained by adding the first step length to the offset imaging result of the current iteration; The input data is seismic data in the first iteration and residuals obtained from the previous iteration in subsequent iterations. The offset imaging result of the previous iteration is data 0 in the first iteration, and the offset imaging result of the previous iteration in subsequent iterations.

2. The method as described in claim 1, characterized in that, Migration calculations are performed on the acquired seismic data to obtain migration imaging results, including: Based on memory size and offset parameters, calculate part or all of the offset group data and part or all of the imaging space of the offset group in the seismic data that a single offset unit task needs to process; Calculate the travel time table based on the offset velocity field; Multiple migration unit tasks are executed through the job node, the seismic data corresponding to the task and the travel time data corresponding to the travel time table are read, the migration imaging is completed and saved; Imaging results from the same offset group are merged to obtain offset imaging results from different offset groups.

3. The method as described in claim 2, characterized in that, The offset calculation process also includes: preconditioning, which involves adding weight coefficients to the offset result. Where a ranges from 0.5 to 2, b ranges from 0.00025 to 0.01, and offset represents the middle offset length of the offset group; the pre-conditioning imaging formula is as follows: ； in, The spatial imaging position is represented by a two-dimensional vector in two-dimensional cases and a three-dimensional vector in three-dimensional cases, where the depth dimension is represented by z. Indicates spatial imaging position Offset grouping The imaging results The weighting coefficients mentioned above are: , Indicates the weighting coefficient. Indicates gun and gun inspection configuration , The spectrum of earthquake data, Indicates offset grouping, In two dimensions, it is a one-dimensional vector; in three dimensions, it is a two-dimensional vector. This represents the filter factor; in the two-dimensional case, it is a semi-derivative filter, i.e.: In the three-dimensional case, it is derivative filtering, that is: , Represents the spectrum of the seismic wavelet. , These represent the distances from the shot point and receiver point to the imaging point, respectively. When they left, In the two-dimensional case, it represents a single integral, and in the three-dimensional case, it represents a double integral.

4. The method as described in claim 1, characterized in that, Using the offset imaging result as input, inverse offset calculation is performed to obtain the reverse offset imaging result, including: The task node executes multiple inverse offset unit tasks, reads the results of the offset imaging, the travel time data in the corresponding travel time table, and the track head information of the data in the input offset distance group, performs the corresponding inverse offset calculation and saves it. The inverse offset imaging results of the same offset group are merged to obtain inverse offset imaging results of different offset groups.

5. The method as described in claim 4, characterized in that, Also includes: In the inverse migration calculation, the imaging results are pre-conditioned along the depth direction. The pre-conditioned inverse migration formula is as follows: ； in, Indicates gun and gun inspection configuration , Inverse offset synthetic data, Indicates offset grouping, In two dimensions, it is a one-dimensional vector; in three dimensions, it is a two-dimensional vector. This represents the filter factor; in the two-dimensional case, it is a semi-derivative filter, i.e.: In the three-dimensional case, it is derivative filtering, that is: , Represents the spectrum of the seismic wavelet. The spatial imaging position is represented by a two-dimensional vector in two-dimensional cases and a three-dimensional vector in three-dimensional cases, where the depth dimension is represented by z. Indicates the weighting coefficient. Indicates spatial imaging position Offset grouping The imaging results , These represent the distances from the shot point and receiver point to the imaging point, respectively. When they left, The weighting coefficients mentioned above are: , In the two-dimensional case, it represents a double integral, and in the three-dimensional case, it represents a triple integral.

6. The method as described in claim 1, characterized in that, The first step length is obtained by dividing the square of the modulus of the offset imaging result by the inverse offset imaging result.

7. The method as described in claim 1, characterized in that, If the preset convergence criteria are met, the final migration imaging result is generated based on the migration imaging results obtained in this iteration, including: If the preset convergence criteria are met, the migration imaging result obtained in this iteration is multiplied by the precondition factor to generate the final migration imaging result.

8. The method as described in claim 1, characterized in that, If the conditions are not met, the offset calculation is performed again using the residual updated in this iteration to obtain the offset imaging result; Based on the currently calculated migration imaging results and the migration imaging results obtained from the previous migration calculation, the conjugate gradient direction is determined, including: The second step length is obtained by dividing the squared modulus of the offset imaging result obtained from the previous offset calculation by the squared modulus of the offset imaging result obtained from the previous offset calculation. Based on the second step size, the offset imaging result obtained from this offset calculation, and the offset imaging result obtained from the previous offset calculation, the conjugate gradient direction is obtained.

9. A pre-stack depth migration device using least squares integration, characterized in that, include: Offset calculation module, inverse offset calculation module, update module, judgment module, result output module, and conjugate calculation module; The migration calculation module is used to perform migration calculations on the acquired seismic data to obtain migration imaging results; when the judgment module determines that the current iteration does not meet the preset convergence criteria, it uses the residual updated after the current iteration to perform migration calculations again to obtain migration imaging results. The reverse offset calculation module is used to take the offset imaging result as input, perform reverse offset calculation, and obtain the reverse offset imaging result. The update module is used to calculate the first step length based on the square of the modulus of the offset imaging result and the square of the modulus of the reverse offset imaging result. Based on the first step length, update the residual of this iteration and the offset imaging result obtained in this iteration; when updating the residual, calculate the cross-correlation coefficient between each residual and the residual update direction. If the cross-correlation coefficient is less than 0, set both the residual and the residual update direction of that trace to 0; the residual update direction is the reverse offset result of that trace. The step of updating the residual of the current iteration and the offset imaging result obtained in the current iteration based on the first step length includes: When the cross-correlation coefficient between the input data at the same location and the reverse migration imaging result is greater than or equal to zero, the input data of that channel is subtracted from the product of the first step length and the reverse migration imaging result of this iteration to obtain the residuals of multiple offset groups in this iteration. When the cross-correlation coefficient is less than zero, the residual of that channel is set to 0. The offset imaging result of the previous iteration is obtained by adding the first step length to the offset imaging result of the current iteration; The input data is seismic data in the first iteration and residuals obtained from the previous iteration in subsequent iterations. The offset imaging result of the previous iteration is data 0 in the first iteration, and the offset imaging result of the previous iteration in subsequent iterations; The judgment module is used to determine whether the current iteration meets the preset convergence criteria; The result output module is used to generate the final offset imaging result based on the offset imaging result obtained in the current iteration when the judgment module determines that the current iteration meets the preset convergence criterion. The conjugate calculation module is used to determine the conjugate gradient direction based on the offset imaging result currently calculated by the offset calculation module and the offset imaging result obtained in the previous offset calculation when the judgment module determines that the current iteration does not meet the preset convergence criterion.

10. A computing device, characterized in that, include: The memory, the processor, and the computer program stored in the memory and executable on the processor, wherein the processor, when executing the program, implements the least squares integral pre-stack depth migration method according to any one of claims 1-9.

11. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the least squares integral pre-stack depth migration method according to any one of claims 1-9.

Images