Three-dimensional tunnel seismic wave full waveform inversion method and system

Abstract

According to the three-dimensional tunnel seismic wave full-waveform inversion method and system, a finite difference parallel computing method based on a 2.5 d space blocking algorithm replaces conventional cyclic finite difference and convolution finite difference, a three-dimensional data volume is divided into two parts, the two parts are arranged in a shared memory, and the third part is arranged in a register, so that the three-dimensional tunnel seismic wave full-waveform inversion is realized. And the thread blocks and the calculation units in cuda are utilized to the maximum extent. According to the method, the problem of overlarge occupation of a memory / video memory based on convolution tunnel full-waveform inversion can be solved, the operation resources of inversion are saved, and the calculation efficiency of inversion is improved.

Description

Technical Field

[0001] The present invention belongs to the field of tunnel seismic wave advance detection data processing and inversion, and specifically relates to a three-dimensional tunnel seismic wave full waveform inversion method and system, especially a three-dimensional tunnel seismic wave full waveform inversion method and system based on accompanying wave field checkpoints. Background Art

[0002] The statements in this section merely provide background information related to the present invention and do not necessarily constitute prior art.

[0003] In tunnel earthquake prediction, obtaining a relatively accurate velocity distribution ahead of the tunnel is a crucial prerequisite for accurately imaging the adverse geological structures ahead. Full waveform inversion (FWI) is currently one of the most accurate velocity inversion methods and has achieved initial application in tunnel earthquake prediction. Unlike conventional oil and gas seismic exploration methods, tunnel earthquake prediction, particularly those constructed by TBMs, can achieve speeds of up to 30 meters per day. This requires a higher level of efficiency in processing and interpreting tunnel earthquake prediction data. Conventional FWI methods, however, are computationally time-consuming and struggle to meet the demands of rapid tunnel detection.

[0004] To improve the efficiency of full-waveform inversion for tunnel seismic analysis, researchers have proposed a full-waveform inversion method based on convolution and GPU parallelism, effectively reducing the computational cost of full-waveform inversion. On the other hand, for three-dimensional tunnel seismic wave multi-shot active source detection, this full-waveform inversion process based on convolution and GPU parallelism also has drawbacks, namely, excessive memory and video memory usage. When the model is large, for example, reaching 200*200*200, with 2000 time steps, it will occupy a very large amount of video memory. Therefore, convolution-based full-waveform inversion is currently mostly used for full-waveform inversion of two-dimensional models. For large three-dimensional models, it will exceed the physical limits of the hardware's memory / video memory, seriously affecting the promotion and application of this method. Summary of the Invention

[0005] In order to solve the above problems, the present invention proposes a three-dimensional tunnel seismic wave full waveform inversion method and system. The present invention can solve the problem of excessive memory / video memory usage based on convolution tunnel full waveform inversion, save inversion computing resources, and improve the computational efficiency of inversion.

[0006] According to some embodiments, the present invention adopts the following technical solutions:

[0007] A three-dimensional tunnel seismic wave full waveform inversion method includes the following steps:

[0008] Set the parameters for tunnel seismic wave full waveform inversion based on the on-site data acquisition situation;

[0009] According to the wave velocity of the surrounding rock around the tunnel, a uniform initial wave velocity model is constructed to calculate the forward propagation process of the active source of multiple shots in the tunnel;

[0010] Set checkpoints according to the earthquake forward modeling time step;

[0011] The initial model seismic records are used to calculate the difference between the shot points and receiver positions of the actual seismic records. The residual reflects the difference between the model prediction and the actual observation. New source data is generated based on the residual. The difference is used as the source, and the forward propagation wave field is calculated again to obtain the accompanying wave field.

[0012] Using the wave field data at the checkpoint, the entire forward time step is divided into multiple segments, and the wave field data at each node is forward calculated sequentially to obtain the reverse propagation wave field data;

[0013] While calculating the companion wavefield, the companion wavefield is integrated, and the acquired back-propagated wavefield data is cross-correlated with the companion wavefield to obtain the wavefield gradient of the current shot.

[0014] Accumulate the wave field gradients of multiple shots to obtain the final gradient;

[0015] The gradient is optimized using the stochastic optimization method with adaptive momentum, and the gradient of each point in the model is continuously iterated into the velocity model to optimize the velocity model.

[0016] As an optional implementation, the parameters include: grid size, sampling rate, sampling time, observation system and surrounding rock wave velocity.

[0017] As an optional implementation method, a finite difference CUDA parallel method is used to calculate the forward propagation process of multiple active sources in the tunnel. The finite difference CUDA parallel method is an algorithm that is a fusion of the 2.5D blocking algorithm and the finite difference algorithm.

[0018] As an optional implementation, the process of fusing the 2.5D blocking algorithm with the finite difference method includes:

[0019] Before performing finite difference calculations, the data required for the calculation process is split into two parts: the xy dimension and the z dimension. The xy dimension is placed in shared memory, and the z dimension is placed in a register. In each time step, the wave field data and the difference coefficients are split and calculated according to the above segmentation method.

[0020] As an optional implementation, the 2.5D blocking algorithm performs blocking on the xy-dimensional plane and places data into registers through the z-dimension for data flow transmission.

[0021] As an optional implementation, the process of setting checkpoints according to the seismic forward modeling time step includes: truncating the entire time step into multiple time nodes, storing a corresponding number of forward propagating wavefield data at the time nodes according to the differential order, and calculating the reverse propagating wavefield data within the time segment between the two nodes based on the forward propagating wavefield data of the two time nodes.

[0022] As a further step, let (s, t) be the maximum length of any computation chain, i.e., the maximum number of time steps. Reversing the computation chain, there are at most s checkpoints, and there are at most t forward steps from any state. Construct an equation to determine the logarithmic dependence of memory requirements and the number of operations with respect to the running time of function evaluation, and determine the checkpoints according to the equation.

[0023] As a further example, the equation is:

[0024]

[0025] The equation shows a logarithmic dependence of the memory requirement and the number of operations relative to the running time of the function evaluation. If we have the values of two of the three quantities, namely the maximum number of steps, the maximum number of checkpoints, and the maximum number of repeated advances, we can calculate the value of the third quantity.

[0026] As an optional implementation method, using the wavefield data at the checkpoint, a whole forward modeling time step is divided into multiple segments, and the wavefield data at each node is sequentially forward modeled to obtain the back-propagated wavefield data. The specific process includes:

[0027] During the inversion operation, reserve Checkpoint memory space, set the first checkpoint as the initial state, where n is the step size of the forward modeling;

[0028] Starting from the last specified checkpoint, advance to the second-to-last time step, execute forward time steps, and do not record intermediate steps; if one or more checkpoints are idle, set as many checkpoints as possible to intermediate states;

[0029] Step forward in time, record the middle part to the step end, scan backward to the second-to-last time step, calculate the adjoint, and if it is a checkpoint, release it for subsequent use;

[0030] The forward propagation wave field is recalculated using the finite difference CUDA parallel method;

[0031] Split a whole forward time step into The wavefield data at each node in the forward sequence are forward calculated to obtain the reverse propagation wavefield data.

[0032] A three-dimensional tunnel seismic wave full waveform inversion system, comprising:

[0033] The initial model building module is configured to construct a uniform initial wave velocity model based on the wave velocity of the surrounding rock around the tunnel and calculate the forward propagation process of the active source of multiple shots in the tunnel;

[0034] A checkpoint setting module is configured to set a checkpoint according to a time step of an earthquake forward modeling;

[0035] The companion wavefield calculation module is configured to use the obtained initial model seismic records and the actual seismic records to calculate the difference between the corresponding shot point and receiver point positions to obtain the residual. The residual reflects the difference between the model prediction and the actual observation. New source data is generated based on the residual, and the difference is used as the source to recalculate the forward propagation wavefield to obtain the companion wavefield.

[0036] The reverse propagation wavefield calculation module is configured to use the wavefield data at the checkpoint to divide a whole forward modeling time step into multiple segments, and sequentially perform forward modeling calculations on the wavefield data at each node in each segment to obtain reverse propagation wavefield data;

[0037] The gradient calculation module is configured to integrate the companion wavefield while calculating it, cross-correlate the acquired back-propagated wavefield data with the companion wavefield to obtain the wavefield gradient of the current shot, and accumulate the wavefield gradients of multiple shots to obtain the final gradient.

[0038] The iteration module is configured to optimize the gradient by adopting a stochastic optimization method of adaptive momentum, and continuously iterate the gradient of each point in the model into the velocity model to optimize the velocity model.

[0039] A computer-readable storage medium is used to store computer instructions, and when the computer instructions are executed by a processor, the steps in the above method are completed.

[0040] An electronic device includes a memory and a processor, and computer instructions stored in the memory and executed on the processor. When the computer instructions are executed by the processor, the steps in the above method are completed.

[0041] Compared with the prior art, the present invention has the following beneficial effects:

[0042] The present invention proposes a three-dimensional tunnel seismic wave full waveform inversion method based on the accompanying wave field checkpoint. In the process of calculating the forward propagation wave field, the full waveform of the seismic wave is intercepted and saved. The wavefield data of each time node is partitioned and the back-propagated wavefield data is calculated in blocks, which effectively reduces the space occupied by the back-propagated wavefield data, solves the problem of excessive space occupation of full waveform inversion for large models with long time steps, and expands the scope of application of full waveform inversion.

[0043] The present invention proposes a finite difference parallel computing method based on a 2.5D spatial blocking algorithm, which replaces conventional cyclic finite differences and convolutional finite differences. By dividing the three-dimensional data volume into two parts, with the second dimension placed in shared memory and the third dimension placed in registers, the thread blocks and computing units in CUDA are utilized to the maximum extent. While achieving the maximum throughput of the CUDA core, the computational efficiency of finite differences is improved as much as possible, and the time required for forward modeling is greatly reduced.

[0044] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, preferred embodiments are given below and described in detail with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] The accompanying drawings, which constitute a part of the present invention, are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute improper limitations on the present invention.

[0046] Figure 1 A schematic diagram of a process flow provided for an embodiment;

[0047] Figure 2 A transformation curve of the objective function versus time offset provided by an embodiment;

[0048] Figure 3 An embodiment provides The optimal value of Schematic diagram of the effective domain;

[0049] Figure 4 An evaluation process provided for one embodiment uses an inversion function schematic;

[0050] Figure 5 A schematic diagram of the runtime for different numbers of checkpoints provided for one embodiment. DETAILED DESCRIPTION

[0051] The present invention will be further described below with reference to the accompanying drawings and embodiments.

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

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

[0054] In the absence of conflict, the embodiments and features in the embodiments of this application can be combined with each other.

[0055] Example 1

[0056] The present invention first introduces a 2.5D blocking algorithm into the finite difference method for solving the three-dimensional acoustic wave equation, and performs segmentation processing and calculation on the data in the finite difference calculation process.

[0057] In order to solve the problem of excessive space occupation by full 3D waveform inversion in tunnel seismic detection mode, the present invention optimizes the traditional full waveform inversion adjoint state method and no longer saves the entire reverse propagation wavefield data. Instead, the reverse propagation wavefield data is calculated and saved in blocks by means of checkpoints.

[0058] The following is a detailed introduction to a three-dimensional tunnel seismic wave full waveform inversion method. Figure 1 As shown, the following steps are included:

[0059] (1) A uniform initial wave velocity model is constructed based on the wave velocity of the surrounding rock around the tunnel, and the forward propagation process of multiple active sources in the tunnel is calculated using the finite difference CUDA parallel method;

[0060] (2) Set checkpoints according to the earthquake forward modeling time step, that is, use The method is used to evenly intercept time nodes, and according to the order of finite difference, a corresponding number of wave field data records and seismic records of all steps are saved at each time node;

[0061] (3) Using the initial model earthquake record obtained in step (2) and the actual earthquake record, the corresponding shot point and receiver point positions are subtracted to calculate the residual. The residual reflects the difference between the model prediction and the actual observation. New source data is generated based on the residual. The difference is used as the source. The forward propagation wave field is calculated again using the finite difference CUDA parallel method to obtain the accompanying wave field.

[0062] (4) While calculating the accompanying wave field in step (3), use the wave field data at the checkpoint stored in step (2) to split the entire forward time step into Segment, forward calculation is performed on the wave field data at each node in the order from back to front to obtain the reverse propagation wave field data;

[0063] (5) While calculating the companion wavefield, integrate the companion wavefield, then cross-correlate the back-propagated wavefield data obtained in step (4) with the companion wavefield to obtain the wavefield gradient of the current shot;

[0064] (6) Accumulate the wave field gradients of multiple shots to obtain the final gradient;

[0065] (7) The gradient is optimized using the stochastic optimization method of adaptive momentum, and the gradient of each point in the model is continuously iterated into the velocity model to achieve the purpose of optimizing the velocity model.

[0066] Among them, in step (1), the forward propagation process of multiple active sources in the tunnel is calculated using the finite difference CUDA parallel method. The finite difference CUDA parallel method is to integrate the 2.5D blocking algorithm with the finite difference. The specific content is introduced below.

[0067] In traditional finite difference and convolution-based finite difference calculations, the throughput of CUDA kernels is limited by the available bandwidth, so CUDA kernels cannot fully utilize the available computing power. Improving kernel throughput is the key to improving the efficiency of calculating three-dimensional data.

[0068] First, let R denote the radius (in units) of the template range, typically defined as the Manhattan distance (e.g., k-point template) or the Lo norm (e.g., LBM). We assume that the 3D data is laid out with the z-axis as the most frequently changing dimension, followed by the x and y directions. Without loss of generality, let (0,0,0) (the origin) be the grid point with the smallest coordinates, and let Nx, Ny, and Nz be the input grid sizes in the x, y, and z directions, respectively. Let P be a grid point, and [P] be the distance from that grid point to the origin. The size of each grid cell is represented by E.

[0069] Assume that all mesh elements in the XY plane are calculated for a specific Z value (initialized to 0). By scaling the Z value from 0 to Zmax, all node data within the 3D mesh is included in the calculation. For any specific Z value (Zs), only the mesh data within the (Zs+ / -R) Z planes is required to reside in the cache, for a total of 2R+1 planes. In theory, this will increase kernel throughput by requiring blocks with larger dimx and dimy to reside in the cache, thereby reducing the additional bandwidth consumed, as described below.

[0070] We use the assumption that only the (2R+1)XY planes need to be cached. Therefore, we perform blocking on the two-dimensional (XY) plane and loop through the third dimension (Z) to put the data into registers for data flow. This is called 2.5D blocking, borrowing the term 2.5D ( Figure 2).make and Indicates the size of the blockage in the X and Y directions respectively.

[0071] The area within each plane is represented as an XY subplane. Since the buffer needs to be completely stored in the cache, the following formula must be satisfied:

[0072]

[0073] In formula (1), R represents the radius of the template range (unit), and denote the block size in the X and Y directions respectively, C denotes the size of the available GPU shared memory and registers, and ε denotes the size of each grid cell.

[0074] 1. First, load the grid elements in the XY sub-plane respectively. Put in middle.

[0075] 2. Cache-friendly template calculation

[0076] (1) For each Next, The XY subplane is loaded into the buffer middle;

[0077] (2) Storage The template calculation is performed on the XY sub-plane in the image and the result is stored in the external memory.

[0078] Note that unlike 3D blocking, there is no additional bandwidth requirement in the Z direction. In addition, the access pattern from external memory is very regular and can be captured by various hardware prefetchers. For the additional bandwidth of the ghost layer (k 2.5D ), the required memory is:

[0079]

[0080] Minimizing the overestimation yields:

[0081]

[0082] In summary, 2.5D blocking helps reduce the extra bandwidth (compared to loading each element once in a small fixed-size on-chip memory) and the corresponding extra computation required to perform stencil calculations, which may significantly speed up computational efficiency and reduce runtime.

[0083] Before performing the finite difference calculation, this embodiment divides the data required for the calculation into two parts, namely the xy dimension and the z dimension. The xy dimension is placed in the shared memory, and the z dimension is placed in the register. In each time step, the wavefield data and the difference coefficient are divided and calculated according to the above division method.

[0084] It should be noted that the setting values of the various parameters in the above embodiments can be adjusted according to specific circumstances.

[0085] In step (2), checkpoints are set according to the earthquake forward modeling time step.

[0086] First, it should be noted that the reverse mode of computing differentials is the discrete analog of the adjoint method known from the calculus of variations. The gradient of a scalar-valued function is generated using the reverse mode (in its basic form) with no more than five times the number of operations required to compute the function itself. Generally speaking, this mode allows the computation of the Jacobian matrix to be at most five times the number of correlations multiplied by the computation of the underlying vector function. However, the spatial complexity of the wavefield backpropagation computation, i.e., its memory requirement, is proportional to the time complexity of the function evaluation itself, as all intermediate results need to be recorded. Therefore, wavefield backpropagation is severely limited by the amount of memory available.

[0087] This embodiment uses the checkpoint method to first intercept the entire time step time nodes, storing a corresponding number of forward propagating wavefield data at the time nodes according to the differential order, and calculating the reverse propagating wavefield data within the time segment between the two nodes based on the forward propagating wavefield data of the two time nodes.

[0088] By choosing the right checkpoints, the spatial complexity of wavefield backpropagation can be reduced from proportional to logarithmic growth. Furthermore, intermediate values at checkpoints are mostly recalculated rather than recorded. To apply this algorithm to wavefield backpropagation calculations, it is necessary to decompose the forward finite difference process into a series of time steps with roughly equal computational effort.

[0089] The main purpose of the "Resolve" checkpoint technology is to better analyze the underground geological structure and characteristics by collecting and processing data at specific checkpoint locations during seismic exploration. The checkpoint scheduling implemented by revolve executes the following loop to implement the inverse function evaluation, such as Figure 4 As shown, the specific process includes:

[0090] (1) Reserve s checkpoint spaces and set the first checkpoint as the initial state.

[0091] (2) Start from the last time step and then decrease by 1 each time until it is less than 2;

[0092] ①Multiple forward: Starting from the last specified checkpoint, forward to the second-to-last step, execute forward time steps, do not record intermediate steps, if one or more checkpoints are idle, then set as many checkpoints as possible to intermediate states;

[0093] ② Combined reverse: forward time step, record the middle part to the step end, reverse scan to the penultimate step, calculate the adjoint, and if the penultimate step is a checkpoint, release it for subsequent use.

[0094] ③End the loop.

[0095] Here, setting as many checkpoints as possible as intermediate states refers to the minimum value of the number of idle checkpoints and the number of states between the current state and the penultimate state.

[0096] Typically, the second restriction only applies to the final stage of the computation, when most checkpoints have been released. In the early stages, the number of available checkpoints may not be large, and the key question is where to set the checkpoints.

[0097] This embodiment adopts an approach of maintaining equal intervals between intermediate states between the current state and the penultimate state.

[0098] Specifically, let (s, t) be the maximum length of any computation chain, i.e., the maximum number of time steps, which can be reversed by the above process and has at most s checkpoints, such as Figure 5 As shown, there are at most t forward steps from any state:

[0099]

[0100] This equation shows a logarithmic dependence of the memory requirements and the number of operations on the runtime of the function evaluation. In addition, if we have values for two of the three quantities: the maximum number of steps, the maximum number of checkpoints, and the maximum number of repeated advances (without recording intermediate steps), we can calculate the value of the third quantity.

[0101] For example, if there are 10 time steps and assuming 3 checkpoints are stored, then plugging into the equation yields t = 2. The figure illustrates this process for the values of s and t and a given number of time steps. The value of β(s, t-1) (here β(3,1) = 4) determines the next checkpoint after the initial checkpoint.

[0102] Therefore, state 4 is also marked as a checkpoint. Finally, the execution of action 1 marks state 7 as a checkpoint because β(s-1, t-1) defines the number of steps from the second checkpoint to the third checkpoint, and β(2, 1) is equal to 3.

[0103] After the checkpoint location and number are set, further considerations include the time and number of steps required to restore the data and how to minimize the number of times the current state is saved to the checkpoint stack.

[0104] To address these issues, this embodiment adopts the following solutions:

[0105] 1. Checkpoint recovery time

[0106] To determine the time required to recover data from a checkpoint, this example considers the following factors:

[0107] (1) Checkpoint frequency: The frequency of checkpoints directly affects the recovery time. If a checkpoint saves the state every T time steps, then the time required to recover from the current state to the second-to-last state will be T.

[0108] (2) Recovery strategy: The recovery strategy will also affect the time. If the stored checkpoints can be accessed quickly, the recovery time will be shorter.

[0109] 2. Reverse the number of steps of the entire chain

[0110] If you want to reverse the state of the entire chain, you need to consider the following factors:

[0111] (1) Time step: Assuming that each time step is Δt, the total number of time steps can be defined as N.

[0112] (2) Forward steps and reverse steps: Forward steps and reverse steps are recorded only once, so they can be expressed by the following formula:

[0113] Number of backward steps = checkpoint interval N = kN, where k is the number of steps between checkpoints.

[0114] 3. Minimum number of checkpoints to save state

[0115] To determine the minimum number of times the current state should be saved to the checkpoint stack, consider the following aspects:

[0116] (1) State change rate: If the state changes quickly, checkpoints need to be saved more frequently; if the change is slow, the save frequency can be reduced.

[0117] (2) Error tolerance: Based on the system's tolerance for errors, the maximum number of steps after each checkpoint is saved can be determined.

[0118] If the tolerance is high, the number of checkpoint saves can be reduced.

[0119] Suppose m represents the total number of time steps to be reversed, then m reverse steps are required. Before each combined reverse step (excluding the first step), the forward scan starts from the current last checkpoint and ends at the penultimate state. Therefore, the data of the last checkpoint needs to be restored.

[0120] Therefore, the number of times to read data from the checkpoint is m - 1 times.

[0121] Let p(m, s) denote the minimum number of additional forward steps required to reverse a sequence of m time steps with s checkpoints stored at any time. Then p(m, s) must satisfy:

[0122]

[0123] Take the explicit form:

[0124] p(m, s) = tm - β(s + 1, t - 1) (6)

[0125] where t is the unique integer such that β(s, t - 1) < m ≤ β(s, t).

[0126] For reversing m > 1 time steps, the first checkpoint is set to the initial state, and the second checkpoint is set to the state of m. Then we must reverse time steps to the right of the second checkpoint, and set at most s - 1 checkpoints, and reverse time steps to the left of the second checkpoint, and set at most s checkpoints. Similarly, we need to take steps forward to reach the state Therefore, p(m, s) must satisfy (5). Equation (6) is represented by induction:

[0127]

[0128] The checkpoint must be set to the initial state. For the first reverse step, at least t forward steps are required. To perform the next reverse step, t - 1 steps forward are needed, and so on. Therefore, the minimum number of forward steps is equal to:

[0129] [[ID=•40]]

[0130] When m = 1, no forward steps are necessary without recording intermediate products because one forward step and one reverse step are sufficient to reverse one time step. Similarly, when t = 0:

[0131] β(s, -1) = 0 < 1 = m = β(s, 0),

[0132]

[0133] This gives s>1 and m>1. They determine the uniqueness of t∈N:

[0134] β(s,t-1) <m≤β(s,t)(10)

[0135] Assume that or and All The assertion is true. Now, if the second checkpoint is set to a state that satisfies both conditions (7) and (8), then equation (6) is valid. If the second checkpoint is set to a state that does not satisfy both conditions, more than tm-β(s+1,t-1) steps forward are required.

[0136] Set the second checkpoint to:

[0137]

[0138] From formulas (4) and (5), we can conclude that:

[0139]

[0140] The valid domain is Figure 2 , which is shaded gray in the figure. From this figure, we can see that there are usually m ranges of choices. We will now show that to reverse m time steps, we need tm-β(s+1,t-1) forward steps.

[0141] From formula (7) and the inductive hypothesis, we can know that:

[0142]

[0143] is equal to the minimum number of forward steps of m time steps. Similarly, from Equation (8) and the inductive hypothesis, we can obtain:

[0144]

[0145] To reverse Time step, need to step forward. Finally, need to take m steps forward to reach

[0146] Therefore, if the second checkpoint is set to a state that satisfies (7) and (8) Then the minimum number of forward steps to reverse m time steps is:

[0147]

[0148] If the second checkpoint is set to a state m that does not satisfy (7) and (8), then the number of forward steps for reversing m time steps will be greater than tm - β(s + 1, t - 1). Note that and is Therefore:

[0149]

[0150] is a convex function, where A convex function can have at most one interval of function minimum. Therefore, there is sufficient evidence to show that when:

[0151] [[ID=…]] [[ID=…]] [[ID=…]]

[0152] The minimum number of forward steps is greater than tm - β(s + 1, t - 1).

[0153] First, we will examine the [[ID=…]] case.

[0154] One possibility is β(s, t - 1) < m - β(s - 1, t - 1). Then we get: …

[0155] … … …

[0156] … and … Therefore, we need the minimum … number of forward steps to reverse … time steps. Similarly, we can conclude that … forward steps are necessary to reverse m - m time steps. Therefore, if the second checkpoint is set to … then the number of forward steps for reversing m time steps is …

[0157] … … …

[0158] … From this, we get: …<…]]

[0159] … … …

[0160] … Therefore, if the second checkpoint is set to the state … then more than p(m, s) forward steps are required. …

[0161] … Now consider the second possibility, when β(s, t - 1) = m - β(s - 1, t - 1), …

[0162] … … …

[0163] … Similarly, when s > 1 and t > 0, … (注:原文中部分连续的省略号内容在翻译时保留原样,因为不清楚完整内容无法准确翻译。)is valid. From the induction hypothesis, we have the need for a total of forward steps to reverse the state to the left of the second checkpoint, forward steps to reverse the position to the right of the second checkpoint time steps. This rate of return is

[0164]

[0165] Therefore, if the second checkpoint is set to state β(s,t-1)+1, we need more than p(m,s) forward steps to reverse m time steps. Finally, let β(s,t-1) be greater than m-β(s-1,t-1), then It can be obtained that when s>1, t>0, it satisfies:

[0166]

[0167] We can now conclude from the inductive hypothesis that the minimum number of forward steps to reverse m time steps is given by Given. Similarly, Forward steps require reversal time steps. Finally, we get:

[0168]

[0169] This shows that if the second checkpoint is set to:

[0170]

[0171] It takes m time steps to go backwards in time to go more than p(m,s) steps forward.

[0172] Secondly, we must also consider In this case, let β(s,t-2) be less than m-β(s-1,t):

[0173]

[0174] Using the inductive hypothesis, we get is the necessary forward step to reverse the left side of the second checkpoint, is the necessary forward step to reverse the right side of the second checkpoint. This leads to:

[0175]

[0176] Therefore, if you select the state As a second checkpoint, more than p(m,s) forward steps must be performed to reverse m time steps. Now assume

[0177] β(s,t-2)=m-β(s-1,t)(31)

[0178] Depend on We can conclude that t>2:

[0179]

[0180] Therefore, at least forward steps to reverse time steps, at least forward steps to reverse time steps. The rate of return is:

[0181]

[0182] Therefore, if the second checkpoint is set to state β(s,t-2)-1, more than p(m,s) forward steps are performed to reverse m time steps. Finally, the case where β(s,t-2)>m-β(s-1,t) will be investigated. Under the conditions. It can be concluded that when t>2

[0183]

[0184] From the inductive hypothesis we get that at least forward steps to reverse time steps, at least forward steps to reverse time steps. Therefore:

[0185]

[0186] Therefore, identity (3) is proven, because if the second checkpoint is set to state Then the minimum number of forward steps to reverse m time steps is also greater than p(m,s).

[0187] In step (7), the gradient can be optimized by using the stochastic optimization method with adaptive momentum, and an existing algorithm such as the Adam algorithm can be selected.

[0188] The checkpointing algorithm described in this embodiment allows for various trade-offs between time and storage requirements. Specifically, by choosing the appropriate inversion parameter, a logarithmic growth in space complexity can be achieved. In any case, a slight increase in the number of operations can significantly reduce storage requirements. Therefore, many applications that require large amounts of storage without requiring checkpoints can be addressed using the inversion function.

[0189] This embodiment optimizes the adjoint state checkpoint method using the 2.5D blocking algorithm, saving memory through checkpoint, that is, the three-dimensional matrix memory space of the n-step operation is changed to The three-dimensional matrix memory operation space is expanded; the 2.5D time blocking algorithm is used to divide the entire forward modeling operation time into several time blocks, which fully utilizes the CUDA parallel architecture to realize parallel computing of forward modeling and effectively improves computing efficiency.

[0190] Full-waveform inversion using the aforementioned adjoint wavefield checkpoint algorithm and 2.5D blocking algorithm was compared with convolution-based full-waveform inversion for a 100*100*100 model, 500 time steps, two shots with 12 channels, and 40 iterations. The proposed algorithm achieved a 9x improvement in time efficiency and a 90% reduction in video memory usage, as shown in Table 1. The test environment was Windows 11, processor: i5-13500, 16GB of RAM, 1TB of hard drive, graphics card: RTX3060, 12GB of video memory. The convolution-based full-waveform inversion could only run for an 80*80*80 model, 500 time steps, two shots with 12 channels, and 40 iterations. However, the improved full-waveform inversion algorithm in this paper could run for a 300*200*200 model, 500 time steps, two shots with 12 channels, and 40 iterations. This demonstrates the significant improvement in both space usage and time efficiency of the proposed algorithm.

[0191] Table 1 Comparison of efficiency and memory usage of different full waveform inversion algorithms

[0192]

[0193] like Figure 2-Figure 3 As shown in FIG, it is a schematic diagram of the execution result of the method provided in this embodiment. It can be seen that the method provided in this embodiment has a small offset of the objective function over time, and for The optimal value of The effective domain is larger.

[0194] Example 2

[0195] A three-dimensional tunnel seismic wave full waveform inversion system, comprising:

[0196] The initial model building module is configured to construct a uniform initial wave velocity model based on the wave velocity of the surrounding rock around the tunnel and calculate the forward propagation process of the active source of multiple shots in the tunnel;

[0197] A checkpoint setting module is configured to set a checkpoint according to a time step of an earthquake forward modeling;

[0198] The companion wavefield calculation module is configured to calculate the target function by subtracting the shot point and receiver point positions corresponding to the initial model seismic record from the actual seismic record, and use the difference as the earthquake source to calculate the forward propagation wavefield again to obtain the companion wavefield;

[0199] The reverse propagation wavefield calculation module is configured to use the wavefield data at the checkpoint to divide a whole forward modeling time step into multiple segments, and sequentially perform forward modeling calculations on the wavefield data at each node in each segment to obtain reverse propagation wavefield data;

[0200] The gradient calculation module is configured to integrate the companion wavefield while calculating it, cross-correlate the acquired back-propagated wavefield data with the companion wavefield to obtain the wavefield gradient of the current shot, and accumulate the wavefield gradients of multiple shots to obtain the final gradient.

[0201] The iteration module is configured to optimize the gradient by adopting a stochastic optimization method of adaptive momentum, and continuously iterate the gradient of each point in the model into the velocity model to optimize the velocity model.

[0202] Example 3

[0203] An electronic device includes a memory and a processor, and computer instructions stored in the memory and executed on the processor. When the computer instructions are executed by the processor, the steps of the method described in embodiment 1 are completed.

[0204] It will be understood by those skilled in the art that embodiments of the present invention may be provided as methods, systems, or computer program products. Thus, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware. Furthermore, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0205] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowcharts and / or block diagrams, as well as combinations of processes and / or blocks in the flowcharts 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, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowcharts and / or block diagrams. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0206] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0207] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0208] The foregoing description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Those skilled in the art will readily appreciate that various modifications and variations of the present invention are possible. Any modifications, equivalent substitutions, or improvements made by those skilled in the art that fall within the spirit and principles of the present invention and do not require creative effort are intended to be within the scope of protection of the present invention.

Claims

1. A three-dimensional tunnel seismic wave full waveform inversion method, characterized by: The following steps are involved: Set the parameters for tunnel seismic wave full waveform inversion based on the on-site data acquisition situation; According to the wave velocity of the surrounding rock around the tunnel, a uniform initial wave velocity model is constructed to calculate the forward propagation process of the active source of multiple shots in the tunnel; Set checkpoints according to the earthquake forward modeling time step; The initial model seismic records are used to calculate the difference between the shot points and receiver positions of the actual seismic records. The residual reflects the difference between the model prediction and the actual observation. New source data is generated based on the residual. The difference is used as the source, and the forward propagation wave field is calculated again to obtain the accompanying wave field. Using the wave field data at the checkpoint, the entire forward time step is divided into multiple segments, and the wave field data at each node is forward calculated sequentially to obtain the reverse propagation wave field data; While calculating the companion wavefield, the companion wavefield is integrated, and the acquired back-propagated wavefield data is cross-correlated with the companion wavefield to obtain the wavefield gradient of the current shot. Accumulate the wave field gradients of multiple shots to obtain the final gradient; The gradient is optimized using the stochastic optimization method with adaptive momentum, and the gradient of each point in the model is continuously iterated into the velocity model to optimize the velocity model.

2. A three-dimensional tunnel seismic wave full waveform inversion method according to claim 1, characterized in that: The parameters include: grid size, sampling rate, sampling time, observation system and surrounding rock wave velocity.

3. The method for full waveform inversion of three-dimensional tunnel seismic waves according to claim 1, wherein: The forward propagation process of multiple active sources of tunnel guns is calculated by using a finite difference CUDA parallel method, wherein the finite difference CUDA parallel method is an algorithm obtained by integrating the 2.5D blocking algorithm with the finite difference algorithm.

4. A three-dimensional tunnel seismic wave full waveform inversion method according to claim 1, characterized in that: The process of fusing the 2.5D blocking algorithm with the finite difference method includes: Before performing finite difference calculations, the data required for the calculation process is split into two parts: the xy dimension and the z dimension. The xy dimension is placed in shared memory, and the z dimension is placed in a register. In each time step, the wave field data and the difference coefficients are split and calculated according to the above segmentation method.

5. The three-dimensional tunnel seismic wave full waveform inversion method according to claim 1, characterized in that: The 2.5D blocking algorithm performs blocking on the xy-dimensional plane and places data into registers for data flow transmission through the z-dimension.

6. A three-dimensional tunnel seismic wave full waveform inversion method according to claim 1, characterized in that: The process of setting checkpoints according to the seismic forward modeling time step includes: cutting the entire time step into multiple time nodes, storing a corresponding number of forward propagating wavefield data at the time nodes according to the differential order, and calculating the reverse propagating wavefield data within the time segment between the two nodes based on the forward propagating wavefield data of two time nodes.

7. A three-dimensional tunnel seismic wave full waveform inversion method as claimed in claim 6, characterized in that: (s, t) is the maximum length of any computation chain, i.e., the maximum number of time steps. Reversing a computation chain results in at most s checkpoints, and there are at most t forward steps from any state. Checkpoints are determined based on equations that determine the logarithmic dependence of memory requirements and the number of operations on the running time of function evaluation.

8. A three-dimensional tunnel seismic wave full waveform inversion method according to claim 7, characterized in that: The equation is: The equation shows a logarithmic dependence of the memory requirement and the number of operations relative to the running time of the function evaluation. If we have the values of two of the three quantities, namely the maximum number of steps, the maximum number of checkpoints, and the maximum number of repeated advances, we can calculate the value of the third quantity.

9. The method for full waveform inversion of three-dimensional tunnel seismic waves according to claim 1, wherein: Using the wavefield data at the checkpoint, the entire forward time step is divided into multiple segments. The wavefield data at each node is forward calculated sequentially to obtain the back-propagated wavefield data. The specific process includes: During the inversion operation, reserve Checkpoint memory space, set the first checkpoint as the initial state, where n is the step size of the forward modeling; Starting from the last specified checkpoint, advance to the second-to-last time step, execute forward time steps, and do not record intermediate steps; if one or more checkpoints are idle, set as many checkpoints as possible to intermediate states; Step forward in time, record the middle part to the step end, scan backward to the second-to-last time step, calculate the adjoint, and if it is a checkpoint, release it for subsequent use; The forward propagation wave field is recalculated using the finite difference CUDA parallel method; Split a whole forward time step into The wavefield data at each node in the forward sequence are forward calculated to obtain the reverse propagation wavefield data.

10. A three-dimensional tunnel seismic wave full waveform inversion system, characterized by: include: The initial model building module is configured to construct a uniform initial wave velocity model based on the wave velocity of the surrounding rock around the tunnel and calculate the forward propagation process of the active source of multiple shots in the tunnel; A checkpoint setting module is configured to set a checkpoint according to a time step of an earthquake forward modeling; The companion wavefield calculation module is configured to use the obtained initial model seismic records and the actual seismic records to calculate the difference between the corresponding shot point and receiver point positions to obtain the residual. The residual reflects the difference between the model prediction and the actual observation. New source data is generated based on the residual, and the difference is used as the source to recalculate the forward propagation wavefield to obtain the companion wavefield. The reverse propagation wavefield calculation module is configured to use the wavefield data at the checkpoint to divide a whole forward modeling time step into multiple segments, and sequentially perform forward modeling calculations on the wavefield data at each node in each segment to obtain reverse propagation wavefield data; The gradient calculation module is configured to integrate the companion wavefield while calculating it, cross-correlate the acquired back-propagated wavefield data with the companion wavefield to obtain the wavefield gradient of the current shot, and accumulate the wavefield gradients of multiple shots to obtain the final gradient. The iteration module is configured to optimize the gradient by adopting a stochastic optimization method of adaptive momentum, and continuously iterate the gradient of each point in the model into the velocity model to optimize the velocity model.

Images