An adaptive offset-aperture marchenko imaging method for large-scale subsurface regions

Abstract

The application belongs to the technical field of exploring the earth and is an adaptive offset aperture Marchenko imaging method for large-scale underground areas, comprising: calculating the travel time field information of preprocessed seismic data; performing longitudinal difference calculation according to the calculated travel time field information to obtain a longitudinal travel time difference field; completing the offset aperture size corresponding to all underground imaging points; calculating the background velocity model gradient to obtain the normal direction of the underground structure interface, making an extension line along the normal method to the surface to obtain an intersection point as the offset aperture center position corresponding to the underground imaging point; calculating the final offset aperture range; screening a local reflection data body; calculating the local focusing function and the local Green function corresponding to the underground imaging point; and performing imaging operation on the underground imaging point according to the imaging condition and the local focusing function and the local Green function.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of exploring the earth, and particularly relates to a self-adaptive offset aperture Marchenko imaging method for a large-scale underground area. BACKGROUND

[0002] The Marchenko imaging method is a new method developed in the field of exploration geophysics in the past decade, and has obvious advantages over the traditional migration imaging method in terms of interlayer multiple waves in the underground. The traditional seismic migration imaging method has multiple wave artifacts in the imaging results when facing the interlayer multiple waves in the seismic data. The focusing function obtained by solving the Marchenko equation correctly extends the surface excitation surface acquisition seismic data to the underground to obtain the Green function of the underground excitation surface acquisition, and the Green function is applied to imaging to effectively suppress the multiple wave artifacts in the imaging.

[0003] As an integral method, the Marchenko method has high requirements for the sampling density of the excitation shot and the receiving acquisition, which leads to the need to process a large amount of data. When facing a large-scale area, the data volume increases exponentially with the increase of the area, the computing resources consumed by imaging also increase exponentially, and the consumption of computer memory also increases exponentially. At the same time, the cost of acquisition according to the standard requirements of the Marchenko method is also huge. Moreover, for seismic data of a large-scale area, the imaging result will be affected by the far offset data interference without considering the offset aperture, resulting in a decrease in the resolution of the imaging result. The offset aperture problem exists in many applications of the conventional seismic imaging method, such as the Kichhoff imaging method. However, the Marchenko method has high requirements for the spatial sampling rate of the seismic data, and the setting of the offset aperture involves specific operations on the reflection data volume, which is more complex than the traditional method. Therefore, the application introduces the concept of self-adaptive offset aperture into the Marchenko method, limits the data by using the offset aperture, optimizes the computing efficiency, and improves the imaging resolution and adaptability of the method in practical application. SUMMARY

[0004] The embodiment of the application provides a self-adaptive offset aperture Marchenko imaging method for a large-scale underground area, which solves the problem of excessive consumption of computing resources by the conventional Marchenko method when facing a large-scale area.

[0005] According to the self-adaptive offset aperture Marchenko imaging method for a large-scale underground area, the method comprises the following steps:

[0006] After removing direct waves from the original seismic data, deconvolution calculation is performed to obtain preprocessed seismic data including reflected waves and multiple waves.

[0007] Based on the background velocity model, the travel time field information of the preprocessed seismic data is calculated;

[0008] Based on the calculated travel time field information, longitudinal subtraction is performed to obtain the longitudinal travel time difference field;

[0009] For underground imaging points, by comparing with the longitudinal travel time field, it is determined whether the shot point coordinates are within the offset aperture. The size of the offset aperture is obtained based on the difference between the farthest shot point and the nearest shot point within the offset aperture. By traversing all underground imaging points, the size of the offset aperture corresponding to all underground imaging points is obtained.

[0010] Calculate the gradient of the background velocity model to obtain the normal direction of the underground structural interface. Extend the normal direction to the surface to obtain the intersection point, which is used as the offset aperture center position corresponding to the underground imaging point.

[0011] Set the maximum offset distance and, for the entire seismic data volume, filter seismic data gathers where the distance between the shot point coordinates and the receiver point coordinates is less than or equal to the maximum offset distance.

[0012] Take the smaller value between the maximum offset distance and the offset aperture size as the final offset aperture size, and calculate the final offset aperture range based on the final offset aperture size and the center position of the offset aperture.

[0013] Subsurface imaging points are selected, and data within the final migration aperture range are extracted from the filtered seismic data trace set as local reflection data volumes;

[0014] The local Marchenko equation is solved using local reflection data volumes to calculate the local focusing function and local Green's function corresponding to the underground imaging points.

[0015] Based on the local focusing function and the local Green's function, imaging calculations are performed on the underground imaging points using imaging conditions;

[0016] Complete the calculation of all underground imaging points to obtain the final imaging result.

[0017] Furthermore, determining whether the shot point coordinates are within the offset aperture includes: matching the longitudinal travel time difference fields at different offset distances with the vertical cell grid travel time at the underground imaging point. Make a judgment based on the multiple.

[0018] Furthermore, the travel time of the vertical direction cell grid is represented as follows: , The vertical depth distance The horizontal distance. a background velocity model.

[0019] Further, the gradient of the background velocity model is calculated, including: calculating the partial derivative of the background velocity model along the horizontal direction and the vertical direction.

[0020] Further, the final migration aperture range is calculated from the final migration aperture size and the migration aperture center position, including:

[0021] the difference between the migration aperture center position and half of the final migration aperture size is the minimum value of the final migration aperture range;

[0022] the sum of the migration aperture center position and half of the final migration aperture size is the maximum value of the final migration aperture range.

[0023] Further, the local focusing function and the local Green's function corresponding to the imaging point are calculated by solving the Marchenko equation using the local reflection data volume, including:

[0024] the local direct wave Green's function at the imaging point is obtained by convolving the travel time field information corresponding to the imaging point in the subsurface with the wavelet at the surface;

[0025] the local direct wave Green's function is inverted along the time axis to obtain the initial local downgoing focusing function, and the initial local upgoing focusing function is directly set to 0;

[0026] a local time window function is created, and the corresponding surface receiver of the local time window function is located within the final migration aperture range;

[0027] the upgoing local focusing function and the downgoing local focusing function are obtained by iteratively calculating the local Marchenko equation using the local reflection data volume, the initial local downgoing focusing function and the local upgoing focusing function, and the local time window function;

[0028] the downgoing local Green's function is calculated using the upgoing local focusing function and the downgoing local focusing function;

[0029] the downgoing local Green's function and the downgoing local focusing function are used for imaging.

[0030] Compared with the prior art, the application has the beneficial effects that: by using the adaptive migration aperture, the local seismic data volume is extracted for imaging, which on the one hand reduces the calculation amount and memory consumption of the Marchenko imaging algorithm, and on the other hand reduces the imaging interference caused by far migration data and improves the imaging resolution by setting a variable migration aperture. BRIEF DESCRIPTION OF DRAWINGS

[0031] Figure 1 A flow chart of the adaptive offset-aperture Marchenko imaging method for large-scale subsurface regions provided by the embodiment of the present application is shown in FIG. 1.

[0032] Figure 2 A diagram of the velocity model and the density model provided by the embodiment of the present application is shown in FIG. 2, (a) is the velocity model, and (b) is the density model.

[0033] Figure 3 A diagram of the used seismic data and wavelet provided by the embodiment of the present application is shown in FIG. 3.

[0034] Figure 4 A diagram of the background velocity model provided by the embodiment of the present application is shown in FIG. 4.

[0035] Figure 5 A diagram of the travel-time field information provided by the embodiment of the present application is shown in FIG. 5.

[0036] Figure 6 A diagram of the vertical travel-time difference field provided by the embodiment of the present application is shown in FIG. 6.

[0037] Figure 7 A diagram of the background velocity gradient provided by the embodiment of the present application is shown in FIG. 7, (a) is the gradient horizontal component, (b) is the gradient vertical component, and (c) is the gradient vector field directly drawn according to the two components.

[0038] Figure 8 A diagram of the offset-aperture center position provided by the embodiment of the present application is shown in FIG. 8.

[0039] Figure 9 A diagram of the reflection data volume and the offset-aperture data volume provided by the embodiment of the present application is shown in FIG. 9, (a) is the full-coverage reflection data volume, and (b) is the offset-aperture data volume.

[0040] Figure 10 The final imaging result provided by the embodiment of the present application is shown in FIG. 10, (a) is the imaging result diagram of the entire model, (b), (c), (d) and (e) are the enlarged display diagrams of the selected local part of the large diagram. DETAILED DESCRIPTION

[0041] In order to make the purpose, technical scheme and advantages of the present application more clear, the present application is further described in detail below with reference to the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application, and are not used to limit the present application.

[0042] Referring to Figure 1 shown, an adaptive offset-aperture Marchenko imaging method for large-scale subsurface regions includes:

[0043] S1 performs deconvolution on the raw seismic data after removing direct waves, resulting in preprocessed seismic data including reflected waves and multiples. First, the seismic data that has undergone direct wave removal and deconvolution is prepared as follows: The seismic data includes reflected wave and multiple wave data, formula The center represents the coordinates of the receiving point. Here are the coordinates of the gun point. For time, its array dimension is , This refers to the number of time samples of earthquake data. It is the total number of gun points. This represents the total number of receiving points. Unlike the conventional Marchenko method, the seismic data here is not required to be fully sampled.

[0044] S2 calculates the travel time field information of the preprocessed seismic data based on the background velocity model;

[0045] The background velocity model here is represented as: The first arrival travel time field was calculated by solving the equation using the finite difference method, yielding the travel time field information for the first arrival of seismic waves at all shot points. The travel time field information is expressed as follows: Coordinates in the formula The vertical depth distance The horizontal distance. These are the coordinates of the gun point;

[0046] S3 performs longitudinal subtraction calculation based on the calculated travel time field information to obtain the longitudinal travel time difference field; the formula is expressed as: , This indicates the longitudinal travel time difference field. The vertical grid spacing of the background velocity model. The vertical depth distance is indicated by the vertical grid spacing. From the travel field information below, we can see the coordinates of each shot point. A two-dimensional travel time field is obtained by subtracting the two-dimensional longitudinal travel time difference field. This step requires calculating the longitudinal travel time difference fields for all shot points.

[0047] For underground imaging points, S4 compares the shot point coordinates with the longitudinal travel time field to determine whether they are within the offset aperture. The offset aperture size is obtained by subtracting the farthest and nearest shot points within the offset aperture. This process is repeated for all underground imaging points to complete the offset aperture size for all underground imaging points.

[0048] For underground imaging points, the following methods are used: This indicates the determination of the shot point coordinates. Whether it is within the offset aperture requires the following condition to be met: Parameters in the formula This parameter adjusts the size of the offset aperture; the larger the parameter, the smaller the overall offset aperture. This determination is made by comparing the longitudinal travel time difference field at different offset distances with the travel time of the vertical cell grid at the underground imaging point. The judgment is made by multiplying the distance; theoretically, the smaller the offset, the shorter the longitudinal travel time difference. Travel time with vertical cell grid The closer they are. For underground imaging points, the coordinates of all shot points that satisfy this criterion are calculated. This leads to the determination of the offset aperture size. ,in These are the coordinates of all shot points that meet the conditions. The maximum value in, and These are the coordinates of all shot points that meet the conditions. The minimum value in, This indicates the offset aperture size. This step calculates the offset aperture size for all underground imaging points following the process described above.

[0049] S5 calculates the gradient of the background velocity model to obtain the normal direction of the underground structural interface. The intersection point is obtained by extending the normal method to the surface, and it is used as the offset aperture center position corresponding to the underground imaging point.

[0050] The gradient of the background velocity model is calculated and expressed as: That is, calculating the background velocity model Partial derivatives along the horizontal and vertical directions and The calculation is performed using spatial difference: , , It is the horizontal grid spacing. (This refers to the vertical grid spacing). Indicates that it is located at Background speed at that location, Indicates that it is located at Background speed at that location, Indicates that it is located at Background speed at that location.

[0051] The normal direction of the subsurface structural interface is the direction of the steepest change in wave velocity, and therefore perpendicular to the gradient vector. To obtain the normal direction of the subsurface structural interface, the gradient vector can be normalized to obtain a unit normal vector, which represents the normal direction of the subsurface structural interface. The intersection point of the line extended from the subsurface imaging point along the gradient direction to the surface is taken as the center position of the offset aperture corresponding to the subsurface imaging point, denoted as... .

[0052] S6 sets the maximum offset distance, and selects the seismic data gathers with the distance between the shot coordinates and the receiver coordinates less than or equal to the maximum offset distance for the whole seismic data volume;

[0053] The maximum offset distance is set for the whole imaging process, and is expressed as The reflection data volume is selected according to the maximum offset distance. For the whole seismic data volume The distance between the shot coordinates and the receiver coordinates less than or equal to the maximum offset distance is selected, that is The selected seismic data gathers are expressed as and are stored;

[0054] S7 takes the smaller value between the maximum offset distance and the offset aperture size as the final offset aperture size, and calculates the final offset aperture range according to the final offset aperture size and the offset aperture center position;

[0055] The final offset aperture range of all imaging points is calculated. For the underground imaging point , the calculated offset aperture size and the set maximum offset distance are compared, and the smaller value is taken as the final offset aperture range, that is expressed as: Then, the final offset aperture range is calculated according to the offset aperture center position , that is, the offset aperture range of the underground imaging point is .

[0056] The final offset aperture range is calculated according to the final offset aperture size and the offset aperture center position, including:

[0057] The difference between the offset aperture center position and half of the final offset aperture size is taken as the minimum value of the final offset aperture range;

[0058] The sum of the offset aperture center position and half of the final offset aperture size is taken as the maximum value of the final offset aperture range.

[0059] is expressed as: ;

[0060] ;

[0061] The offset aperture ranges corresponding to all imaging points are stored, is the left boundary coordinate of the offset aperture, and is the right boundary coordinate of the offset aperture.

[0062] S8 selects the underground imaging point, extracts the data in the final offset aperture range from the screened seismic data trace set as the local reflection data body;

[0063] For example, the underground imaging point is selected , the offset aperture range of the underground imaging point is calculated , the local reflection data is extracted from the screened seismic data trace set according to the following conditions, and the local reflection data is represented as: :

[0064] ;

[0065] The obtained local reflection data body is a three-dimensional array with a small amount of data, and the dimensions are , wherein , is the local source, is the number of local receivers, is the source interval or receiver interval, and is the number of seismic data time samples. It should be noted that the source is the shot point, and the receiver is the receiving point.

[0066] S9 uses the local reflection data body to solve the local Marchenko equation, calculates the local focusing function and the local Green function corresponding to the imaging point;

[0067] S10 performs imaging operation on the underground imaging point according to the local focusing function and the local Green function using the imaging condition;

[0068] S11 completes the operation of all underground imaging points to obtain the final imaging result.

[0069] In an embodiment, the obtained local reflection data body is used to solve the Marchenko equation to calculate the local focusing function and the local Green function corresponding to the imaging point, including:

[0070] The travel time field information corresponding to the surface receiving point of the underground imaging point located in the final offset aperture range is convolved with the wavelet to obtain the local direct wave Green function excited at the underground imaging point and received at the surface receiving point coordinate;

[0071] The local direct wave Green function is inverted on the time axis as the initial local downlink focusing function, and the initial local uplink focusing function is directly set to 0;

[0072] A local time window function is created, and the corresponding surface receiving point of the local time window function is located in the final offset aperture range;

[0073] ​Using the local reflection data volume, the initial local down-going focusing function and the local up-going focusing function and the local time window function, the up-going local focusing function and the down-going local focusing function are obtained by iterative calculation according to the local Marchenko iterative formula;

[0074] The down-going local Green's function is calculated using the up-going local focusing function and the down-going local focusing function;

[0075] The imaging is performed using the down-going local Green's function and the down-going local focusing function.

[0076] Wherein, the travel time field information is obtained according to the travel time field information , the underground imaging point is obtained by folding the travel time field information and the wavelet received on the surface, and the local direct wave Green's function excited at the underground imaging point and received at the surface receiving point coordinate is obtained, and the local direct wave Green's function is expressed as , , wherein the underground imaging point is represented by , and the surface receiving point coordinate is represented by , which satisfies: , and the local direct wave Green's function is inverted on the time axis as the initial local down-going focusing function, that is . The superscript in the formula represents the local variable within the offset aperture range, the superscript in the formula represents the down-going propagation direction, in the formula represents the initial local down-going focusing function, and the subscript in the formula represents the iteration number, and the initial local down-going focusing function is determined according to the above method, and the initial local up-going focusing function is directly set to 0.

[0077] According to the travel time field information between the underground imaging point and the surface receiving point coordinate, the local time window function is created according to the following rules:

[0078] ,

[0079] , wherein the local time window function is represented by , the time is represented by , and the local time window function is a relatively small positive number, which is generally set to half of the wavelet time length, and the surface receiving point coordinate of the local time window function is also located within the offset aperture range, and the meaning of the local time window function is to assign zero to the values exceeding the time limit within the local offset aperture range.

[0080] Using the local reflection data volume, the initial local downlink focusing function, the initial local uplink focusing function, and the local time window function obtained in the previous step, the local uplink focusing function and the local downlink focusing function are iteratively calculated according to the following local Marchenko iterative formula:

[0081] ,

[0082] ,

[0083] Wherein, the spatial integral corresponds to It is the range of offset apertures corresponding to the underground imaging points. Represents a time variable. express Local up-focusing function of -1 iterations, express The local upward focusing function of the next iteration. express Local downlink focusing function in the next iteration and Both represent local time window functions. express Local downlink focusing function of -1 iterations, express The reversal, This represents the direct portion of the local downlink cluster function, which is equal to the initial local downlink cluster function.

[0084] Based on the iterated upper local focusing function and the lower local focusing function, the lower local Green's function is calculated using the following formula:

[0085] , This represents the downlink local Green's function.

[0086] Calculate the imaging values ​​at the underground imaging points;

[0087] , Indicates underground imaging points The image value at that location.

[0088] Repeat this process until all underground imaging points have been processed to obtain the final imaging result.

[0089] Taking model testing as an example, the following methods are used: Figure 2 The velocity model and density model shown are as follows, in which Figure 2 (a) in the figure represents the velocity model. Figure 2(b) in the diagram represents the density model, which applies the imaging method proposed in this application to demonstrate its feasibility and advantages. The model has a lateral distance of 32 km and a depth of 9 km. This model is large-scale, and using full-coverage sampled seismic data would result in significant memory and computational resource consumption.

[0090] Preprocessed seismic data after removing direct waves and performing deconvolution, such as... Figure 3 As shown in the left-hand figure, the wavelet morphology after deconvolution is as follows: Figure 3 As shown in the right-hand image;

[0091] Prepare a background velocity model, such as Figure 4 As shown, the background velocity model is derived from the conventional seismic imaging process. Based on this velocity model, the travel time field information is calculated by solving the equations using finite difference calculus. Figure 4 The image shows the travel time field information of an earthquake located 15 km from the surface, which served as the shot point.

[0092] Based on the travel time field information, the longitudinal travel time difference field was calculated. Taking underground spatial points at lateral distances of 15km and depths of 1km, 2km, 3km, and 4km as examples, the travel times of all surface shot points corresponding to the four underground imaging points were subtracted from the travel times corresponding to 10m below the four underground imaging points. The results are as follows: Figure 6 As shown, it can be observed that the larger the offset, the smaller the travel time difference. Based on the change in travel time difference, preset parameters... According to the formula If the coordinates of the shot points that meet this condition are within the offset aperture, the maximum and minimum coordinates of the corresponding offset aperture are obtained after judgment, thus determining the size of the offset aperture. The offset aperture sizes for all underground imaging points are then calculated.

[0093] The gradient of the background velocity model is calculated, and the gradient result is as follows: Figure 7 As shown, Figure 7 (a) and Figure 7 In (b), the horizontal and vertical components of the gradient are shown respectively. Figure 7 (c) shows the gradient vector field plotted directly using the two components. The intersection of the extension line along the gradient vector direction and the Earth's surface is selected as the offset aperture center position. After calculating the offset aperture center positions for all imaging points... Figure 8 This is a schematic diagram for determining the position of the offset aperture center;

[0094] Set the maximum offset distance to filter the reflection data volume, such as Figure 9 As shown, Figure 9 As shown in (a), each point represents one seismic data point in the full-coverage reflection data volume. The horizontal axis corresponds to the receiver coordinates, and the vertical axis corresponds to the shot location. Here, with the maximum offset set to 10, as shown... Figure 9As shown in (b) of FIG. 1, the reflection data volume becomes the migration aperture data volume, at this time, the reflection data volume saves a large amount of storage space compared to the previous one, and the reflection data volume is stored for subsequent imaging calculation. Then, according to the set maximum offset distance, the final migration aperture size and the migration aperture center position of the previously calculated underground imaging point, the final migration aperture range corresponding to all underground imaging points is calculated.

[0095] Select an underground imaging point, extract a local reflection data volume according to the final migration aperture range, for example, a point at a horizontal distance of 3 km and a depth of 1 km, according to the travel time field information (such as Figure 5 As shown in the travel time field information, the initial local downlink focusing function and the local time window function are calculated by using the travel time field information and the wavelet convolution, and the wavelet here refers to the Ricker wavelet.

[0096] The initial local downlink focusing function and the local uplink focusing function and the local time window function are used to solve the uplink local focusing function and the downlink local Green function by using the extracted local reflection data volume and applying the local Marchenko equation.

[0097] The downlink local Green function is calculated by using the uplink local focusing function and the downlink local focusing function.

[0098] The downlink local Green function and the downlink local focusing function are used for imaging, and the result is shown in Figure 10 Figure 10 The imaging result of the whole model is shown in (a) of FIG. 2, Figure 10 (b) of FIG. 2, Figure 10 (c) of FIG. 2, Figure 10 (d) of FIG. 2, Figure 10 (e) of FIG. 2 is a zoomed-in display diagram of a local part of the large diagram, which contains the local imaging results of each depth range. The imaging result is clear, and there is no low-frequency noise interference caused by far offset distance data.

[0099] The above only describes the preferred embodiments of the present application and is not used to limit the present application. Any modification, equivalent replacement and improvement made within the spirit and principle of the present application shall be included in the protection scope of the present application.​

Claims

1. An adaptive migration aperture Marchenko imaging method for large-scale subsurface regions, characterized in that, include: After removing direct waves from the original seismic data, deconvolution calculation is performed to obtain preprocessed seismic data including reflected waves and multiple waves. Based on the background velocity model, the travel time field information of the preprocessed seismic data is calculated; Based on the calculated travel time field information, longitudinal subtraction is performed to obtain the longitudinal travel time difference field; For underground imaging points, by comparing with the longitudinal travel time field, it is determined whether the shot point coordinates are within the offset aperture. The size of the offset aperture is obtained based on the difference between the farthest shot point and the nearest shot point within the offset aperture. By traversing all underground imaging points, the size of the offset aperture corresponding to all underground imaging points is obtained. Calculate the gradient of the background velocity model to obtain the normal direction of the underground structural interface. Extend the normal direction to the surface to obtain the intersection point, which is used as the offset aperture center position corresponding to the underground imaging point. Set the maximum offset distance and, for the entire seismic data volume, filter seismic data gathers where the distance between the shot point coordinates and the receiver point coordinates is less than or equal to the maximum offset distance. The smaller value between the maximum offset distance and the offset aperture size is taken as the final offset aperture size. The final offset aperture range is calculated from the final offset aperture size and the center position of the offset aperture. Subsurface imaging points are selected, and data within the final migration aperture range are extracted from the filtered seismic data trace set as local reflection data volumes; The local Marchenko equation is solved using local reflection data volumes to calculate the local focusing function and local Green's function corresponding to the underground imaging points. Based on the local focusing function and the local Green's function, imaging operations are performed on the underground imaging points using imaging conditions; Complete the calculation of all underground imaging points to obtain the final imaging result.

2. The adaptive migration aperture Marchenko imaging method for large-scale underground regions according to claim 1, characterized in that, Determining whether the shot point coordinates are within the offset aperture includes: matching the longitudinal travel time difference fields at different offset distances with the vertical cell mesh travel time at the subsurface imaging point. Make a judgment based on the multiple.

3. The adaptive migration aperture Marchenko imaging method for large-scale underground regions according to claim 2, characterized in that, The travel time of the vertical direction cell grid is represented as follows: , The vertical depth distance The horizontal distance. This is the background velocity model.

4. The adaptive migration aperture Marchenko imaging method for large-scale underground regions according to claim 1, characterized in that, Calculate the gradient of the background velocity model, including: calculating the partial derivatives of the background velocity model along the horizontal and vertical directions.

5. The adaptive migration aperture Marchenko imaging method for large-scale underground regions according to claim 1, characterized in that, The final offset aperture range is calculated from the final offset aperture size and the offset aperture center position, including: The minimum value of the final offset aperture range is obtained by taking the difference between the center position of the offset aperture and half the size of the final offset aperture. The maximum value of the final offset aperture range is obtained by summing the offset aperture center position with half of the final offset aperture size.

6. The adaptive migration aperture Marchenko imaging method for large-scale subsurface regions according to claim 1, characterized in that, Solving the Marchenko equation using local reflection data volumes, and calculating the local focusing function and local Green's function corresponding to the imaging point, including: The travel time field information received at the surface corresponding to the underground imaging point located within the final offset aperture range is convolved with the wavelet to obtain the local direct-reach Poinlin function excited at the underground imaging point and received at the coordinates of the receiving point on the surface. The local direct Bogle function is inverted on the time axis and used as the initial local downlink focusing function, while the initial local uplink focusing function is set directly to 0. Create a local time window function, wherein the surface receiving point corresponding to the local time window function is located within the final offset aperture range; Using the local reflection data volume, the initial local downlink focusing function and the local uplink focusing function, and the local time window function, the uplink local focusing function and the downlink local focusing function are obtained by iterative calculation based on the local Marchenko equation. The downlink local Green's function is calculated using the uplink local focusing function and the downlink local focusing function; Imaging is performed using the downlink local Green's function and the downlink local focusing function.

Images