Earth and rockfill dam hidden danger full-waveform inversion method based on simulation inter-well observation system

Abstract

The invention discloses an earth and rockfill dam hidden danger full-waveform inversion method based on a simulated inter-well observation system, and the method comprises the steps: laying excitation points and detection points at a dam crest and upstream and downstream dam slopes along the cross section of an earth and rockfill dam, and constructing the simulated inter-well observation system to collect seismic data; establishing a hidden-danger-free initial model, performing forward modeling as a current model, and constructing a target function based on observation and simulation records; when the target function value does not reach the preset threshold value, performing wave field back-stepping through the residual record, and calculating the physical property parameter gradient of the current model in combination with the forward modeling wave field; processing by adopting a conjugate gradient method to obtain a conjugate gradient, and then constructing a precondition operator for suppressing surface wave interference and accelerating deep updating in a dam body to correct the precondition operator; and obtaining an optimal inversion step length through step length search, updating the model, and carrying out iteration until an objective function is converged, so as to obtain a physical property parameter model reflecting the earth and rockfill dam hidden danger. According to the invention, the accuracy and inversion efficiency of earth and rockfill dam hidden danger detection can be effectively improved.

Description

Technical Field

[0001] This invention relates to the field of seismic wave full waveform inversion, and more specifically, to a method for full waveform inversion of earth-rock dam hazards based on a pseudo-well observation system. Background Technology

[0002] Earth-rock dams are the most widely used dam type in my country, accounting for over 90% of the country's total dams. Their safety is crucial for the normal operation of various water conservancy projects. Compared to concrete dams, earth-rock dam filling materials are inherently more susceptible to erosion and have weaker resistance to damage. Furthermore, early earth-rock dam design standards were often low, construction techniques and equipment were relatively outdated, and effective management and maintenance measures were lacking during long-term operation. This has resulted in many dams exhibiting hidden dangers such as seepage and slow-moving loose layers after decades of operation. If these issues are not detected and addressed promptly, they will not only significantly reduce the stability of the earth-rock dam structure and accelerate aging and damage, but also seriously affect flood control and irrigation functions, posing a potential threat to the safety of people and property downstream and the ecological environment.

[0003] Seismic wave methods, as an efficient, non-destructive, and low-cost geophysical exploration method, are widely used for dam hazard detection. The key to using seismic wave methods to detect dam hazards and assess dam safety lies in obtaining the P-wave or S-wave velocity profiles within the dam body through the analysis and processing of seismic data. Common methods for obtaining dam velocity profiles include reflected wave velocity analysis, surface wave inversion, and full waveform inversion. Full waveform inversion, with its comprehensive utilization of complete wavefield information and high spatial resolution, is more suitable for constructing velocity models of earth-rock dams compared to other methods.

[0004] For detecting potential hazards in earth-rock dams, the conventional method of seismic data acquisition involves arranging geophones along the longitudinal profile of the dam body, such as... Figure 1 As shown. Due to the small size of the potential hazard area inside the dam body and the small difference in physical properties between it and the dam body, the reflected wave field energy of the hazard area is weak. Therefore, the seismic data acquired using traditional observation systems contains relatively little wave field information related to the dam hazard. Using this seismic data for full waveform inversion results in low accuracy in identifying potential hazard areas inside the dam body. Well-based seismic methods are one of the highest resolution methods in seismic exploration. Compared to traditional surface-based seismic exploration, well-based observation systems have a wider coverage area and acquire richer seismic data. Using this data for full waveform inversion yields more accurate velocity profiles. However, conducting well-based seismic exploration requires drilling holes on both sides of the exploration area, which is not only costly but also damages the dam structure.

[0005] The nonlinear preconditioning conjugate gradient method has become one of the most widely used inversion algorithms due to its high inversion accuracy and good stability. One of the most important aspects of this method in the full waveform inversion of earth-rock dam hazards is calculating the model update gradient. Gradient calculation mainly involves two steps: first, calculating the gradient values ​​for the entire model region normally; second, constructing preconditioning operators to correct the gradients, highlighting the gradient values ​​in key areas to improve the efficiency and accuracy of the inversion. The preconditioning operators need to be set according to the actual inversion situation and the propagation characteristics of seismic waves inside the dam body. Currently, there is no preconditioning operator calculation scheme specifically for the full waveform inversion of earth-rock dam hazards.

[0006] In summary, the traditional data acquisition and inversion methods for detecting potential hazards in earth-rock dams using seismic wave methods have shortcomings, which hinder the development of high-precision and high-efficiency technologies for detecting potential hazards in earth-rock dams. Summary of the Invention

[0007] The purpose of this invention is to propose a full waveform inversion method for hidden dangers in earth-rock dams based on a pseudo-well observation system, so as to improve the efficiency and accuracy of hidden danger inversion in earth-rock dams.

[0008] To achieve the above objectives, this invention proposes a full waveform inversion method for hidden dangers in earth-rock dams based on a pseudo-well observation system, comprising: S1. Along the cross section of the earth-rock dam, excitation points and detector points are set up on the dam crest and the upstream and downstream slope surfaces to construct a pseudo-well observation system and collect data to obtain observed seismic records. S2. Establish an initial physical property parameter model of the dam body without hidden dangers, and use this initial model as the current model; S3. Perform forward modeling on the current model to obtain simulated seismic records and forward wavefields; S4. Construct an objective function based on the observed seismic records and the simulated seismic records, and calculate the current objective function value; S5. Determine whether the current objective function value is less than the preset threshold; if yes, output the current model as the inversion result; if no, proceed to step S6. S6. Based on the residuals between the observed seismic records and the simulated seismic records, the wavefield is inversely derived to obtain the inverse wavefield, and the gradient of the physical property parameters of the current model is calculated in combination with the forward wavefield. S7. Use the conjugate gradient method to process the gradient of the physical property parameters of the current model to obtain the conjugate gradient; S8. Construct a preconditional operator for suppressing surface wave interference and accelerating the updating of physical property parameters related to hidden dangers in the deep region of the dam body. Use the preconditional operator to correct the conjugate gradient to obtain the updated gradient. S9. Determine the optimal inversion step size for updating the current model, and update the current model according to the update gradient and the optimal inversion step size to obtain a new model; S10. The new model is used as the new current model, and the process returns to step S3 for iteration until the current objective function value in step S5 is less than the preset threshold, thereby obtaining a physical property parameter model that reflects the hidden dangers of earth-rock dams.

[0009] Optionally, in step S1, the simulated well observation system uses the upstream and downstream slopes of the dam as natural excitation and receiving wells, and sets up excitation points and receiver points along the surface of the dam at preset intervals. The source uses Ricker wavelet with a preset dominant frequency to collect direct longitudinal and transverse waves penetrating the dam, reflected and converted waves from the bedrock interface, reflected waves from the seepage area, and multiple reflected waves generated by the free boundary between the dam slope and the dam crest.

[0010] Optionally, in step S2, the geometry of the initial physical property parameter model is consistent with that of the earth-rock dam being probed.

[0011] Optionally, in step S3, the forward modeling simulation uses the staggered grid finite difference method to simulate the seismic wavefield; In the differential calculation process, the normal stress is set at integer grid points, and the shear stress and particle vibration velocity components are set at half grid points; The top surface of the dam body and the upstream and downstream slopes are set as free boundaries, which are set at half grid points, and the area above the free boundaries is set as vacuum; the bedrock around the dam body and the bottom interface are set as absorbing boundaries.

[0012] Optionally, in step S4, the objective function is the normalized L2 norm based on the waveform amplitude, which is defined as:

[0013] in, u i,j Indicates the first i Cannon j Simulated seismic records from receivers d i,j Indicates the first i Cannon j Seismic records observed at receiver points n.s. For the total number of cannons, nr The total number of detector points, E The objective function value represents the normalized error between observed and simulated seismic records.

[0014] Optionally, in step S6, calculating the gradient of the physical property parameters of the current model specifically includes: By using the residual between observed and simulated seismic records as the excitation source, the wave field is inversely derived to obtain the inverse normal stress wave field and shear stress wave field. The gradient of the Lamé parameters is obtained by cross-correlation calculation of the forward-modeled wavefield and the inverse-modeled wavefield:

[0015] in, l , m Let Lamé constant be . t xx , t zz , t xz These represent the horizontal normal stress wave field, the vertical normal stress wave field, and the shear stress wave field in the forward modeling simulation. or xx , or xx , or xz These are the horizontal normal stress wave field, the vertical normal stress wave field, and the shear stress wave field, respectively, derived in reverse. g λ Lamé constant l gradient, g μ Lamé constant m gradient, sources Indicates all firing points; The gradients of the Lamé parameters are converted into gradients of P-wave velocity, S-wave velocity, and density using the chain rule, thus obtaining the gradients of the physical property parameters of the current model:

[0016] in, V p For the longitudinal wave velocity, V s For transverse wave velocity, r For density, x , z These represent the horizontal and vertical displacement fields derived in reverse, respectively. v x , v z These are the horizontal and vertical velocity fields of the particle vibration, respectively. g Vp The gradient of the longitudinal wave velocity. g Vs The gradient of the transverse wave velocity. g ρ The gradient is the density.

[0017] Optionally, in step S7, the conjugate gradient method is used to process the gradients of the physical property parameters of the current model to obtain the conjugate gradients, specifically including: When the number of iterations n When = 1, the conjugate gradient δg 1 k = g 1 k ,in k Representing different physical property parameters, k When =1, g 1 1 g represents the first iteration Vp , k When =2, g 1 2 g represents the first iteration Vs , k When =3, g 1 3 g represents the first iteration ρ ; when n When ≥2, conjugate gradient δg n k The calculation formula is as follows:

[0018] in, g n k Let be the gradient of the physical property parameters in the nth iteration. The gradient of the physical property parameters in the (n-1)th iteration. δg n k The conjugate gradient of the nth iteration is... The conjugate gradient of the (n-1)th iteration. β n The conjugate correction coefficient. d For conjugate identifiers, This represents the square of the L2 norm.

[0019] Optionally, in step S8, the preconditioning operator is a piecewise function, and its expression is:

[0020] in, p For preconditioning operators, Z Z is the depth downwards from the dam crest as the origin. gradt1 Z is the first depth threshold. gradt2 The second depth threshold, H Z is the height of the dam. gradt1 <Z gradt2 < H, △ l Δ represents the depth range of the linear transition region. l =Z gradt2 -Z gradt1 , a、a 1 represents the depth adjustment factor.

[0021] Optionally, in step S9, determining the optimal inversion step size for updating the current model specifically includes: The optimal inversion step size is obtained by using the parabolic fitting method to search for the inversion step size.

[0022] Optionally, in step S9, the calculation formula for updating the current model based on the update gradient and the optimal inversion step size is as follows:

[0023] in, m n For the current model, m n+1 For the updated new model, α n To achieve the optimal inversion step size, p For preconditioning operators, d g n k is the conjugate gradient of the nth iteration.

[0024] The beneficial effects of this invention are as follows: 1. The full waveform inversion method for earth-rock dams based on a pseudo-well observation system proposed in this invention fully utilizes the geometric advantages of the dam body. By acquiring seismic data along the cross-section of the dam body, it can not only directly obtain direct shear and P waves penetrating the dam body, but also collect multiple reflected waves generated by the free boundaries of the dam slope and crest, thus including more information about potential hazard areas in the acquired seismic records. For the inversion of potential hazards in earth-rock dams, the more information about potential hazard areas contained in the seismic data, the more accurate the inversion results. Therefore, the method proposed in this invention can effectively improve the accuracy of potential hazard detection in earth-rock dams.

[0025] 2. The gradient of the inversion determines the magnitude and direction of the model update, significantly impacting the accuracy and efficiency of the inversion. The preconditioning operator constructed in this invention can effectively correct the inversion gradient, suppress the interference of surface waves on the inversion, and highlight the update speed of the model in the deep region of the dam body, thereby improving the efficiency and accuracy of the inversion. This preconditioning operator is suitable for inverting potential hazards in dams of different structural types. During the inversion process, only relevant parameters (depth adjustment factor) need to be adjusted according to the height of the dam body. a 1) values ​​are set to ensure that the gradient values ​​in the middle and deep regions are within a reasonable range.

[0026] The present invention has other features and advantages, which will be apparent from or will be set forth in detail in the accompanying drawings and the following detailed description, which together serve to explain the particular principles of the invention. Attached Figure Description

[0027] The above and other objects, features and advantages of the present invention will become more apparent from the accompanying drawings, in which like reference numerals generally denote like parts.

[0028] Figure 1 This is a schematic diagram of the survey line layout for seismic wave detection of potential hazards in earth-rock dams. The location of the excitation point is also arranged along the downstream dam slope-dam crest-upstream dam slope. Figure 2 This is a technical roadmap for the full waveform inversion method of earth-rock dam hidden dangers based on the pseudo-well seismic observation system; Figure 3a and Figure 3b These are the actual P-wave velocity model and the actual S-wave velocity model for an earth-rock dam with potential leakage. Figure 4 This is a schematic diagram showing the distribution of the positions of each parameter in the grid points during the differential calculation. The red dashed line represents the free boundary of the dam body. Figure 5 A schematic diagram of the boundary treatment scheme for an earth-rock dam; Figure 6a and Figure 6b These are the VX and VZ component seismic records of a single shot located at a horizontal distance of 34m from the shot point; Figure 7 The initial dam body physical property parameter model is free of hidden dangers; Figure 8a and Figure 8b The residual seismic records for the vx and vz components of a single shot located at a horizontal distance of 34m are shown respectively. Figure 9 This is a schematic diagram of the inversion preconditioner operator; Figure 10a and Figure 10b These are the inverted P-wave velocity model and S-wave velocity model of the dam body, respectively; Figure 11a and Figure 11b These are the horizontal P-wave velocity profile and the P-wave velocity profile at Z=21.6, respectively, of the inverted P-wave and P-wave velocity models. Detailed Implementation

[0029] This invention addresses the shortcomings of traditional seismic wave methods in detecting potential hazards in earth-rock dams. Based on the geometric characteristics of earth-rock dams, it proposes a full-waveform inversion method for detecting potential hazards in earth-rock dams using a pseudo-well-to-well observation system. This method uses the upstream and downstream slopes of the dam body as natural excitation and receiving wells, and deploys geophones along the cross-section of the dam body, such as... Figure 1 As shown (the location of the excitation point is the same as the location of the receiver point in the figure, also arranged along the downstream dam slope-dam crest-upstream dam slope), it can not only directly acquire direct shear and P waves penetrating the dam body, but also collect multiple reflected waves generated by the free boundaries of the dam slope and dam crest. This allows the acquired seismic records to contain more information about potential hazard areas, improving the accuracy of the inversion. During the full waveform inversion calculation, based on the propagation characteristics of seismic waves along the longitudinal profile of the dam body and the location distribution of potential hazards within the dam body, a precondition operator suitable for inversion of simulated well-to-well seismic data was constructed. This operator is used to suppress the influence of surface waves on gradient calculation, accelerate the update speed of physical parameters in the deep areas of the dam body, and improve the efficiency and accuracy of the inversion.

[0030] The invention will now be described in more detail with reference to the accompanying drawings. While preferred embodiments of the invention are shown in the drawings, it should be understood that the invention can be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that the invention will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.

[0031] According to the present invention, a full waveform inversion method for hidden dangers in earth-rock dams based on a pseudo-well observation system is proposed. First, according to Figure 1 The survey line method shown involves setting up excitation and receiving points along the cross-section of the dam to collect seismic records. Then, a dam without hidden dangers is used as the initial model. The seismic records of the initial model are simulated using forward modeling, and the residuals between the observed records and the initial records are calculated. The inverse wavefield is then obtained using the residual records. By cross-correlating the forward and inverse wavefields of the initial model, the conventional inversion gradient value is calculated. Subsequently, based on the propagation characteristics of seismic waves along the longitudinal section of the dam and the location distribution of hidden dangers within the dam, a precondition operator for inversion is constructed. This operator is used to correct the gradient value, obtaining the model update gradient. The initial model is then iteratively updated until the difference between the observed records and the forward records of the updated model is less than a given threshold. The updated model then becomes the physical property parameter model required for inversion. The technical flow of this invention is as follows: Figure 2 As shown. The specific steps for implementing the method of the present invention are as follows: S1. Along the cross section of the earth-rock dam, excitation points and detector points are set up on the dam crest and the upstream and downstream slope surfaces to construct a pseudo-well observation system and collect data to obtain observed seismic records. In this step, the simulated well observation system uses the upstream and downstream slopes of the dam as natural excitation and receiving wells. Excitation points and receiver points are arranged along the surface of the dam at preset intervals. The source uses Ricker wavelet with a preset dominant frequency to collect direct longitudinal and transverse waves penetrating the dam, reflected and converted waves from the bedrock interface, reflected waves from the seepage area, and multiple reflected waves generated by the free boundary between the dam slope and the dam crest.

[0032] In one example, based on the structural characteristics of an earth-rock dam, the following layout is provided along the cross-section of the dam body: nr At the geophone detection point, seismic waves are then generated along the cross section of the dam at a certain shot spacing, and a total of [number] seismic waves are obtained. n.s. Seismic records observed by artillery d .

[0033] S2. Establish an initial physical property parameter model of the dam body without hidden dangers, and use this initial model as the current model; In this step, the geometry of the initial physical property parameter model is consistent with that of the earth-rock dam being probed.

[0034] In one example, an initial physical property parameter model of the dam body without potential hazards is established. m 1. Initial model physical property values ​​(longitudinal wave velocity) V p transverse wave velocity V s ,density r The physical property parameters of the detected dam body are consistent with those of the areas without hidden dangers, and the geometric structure is consistent with that of the detected dam body.

[0035] S3. Perform forward modeling on the current model to obtain simulated seismic records and forward wavefields; In this step, the forward modeling uses the staggered grid finite difference method to simulate the seismic wave field. During the difference calculation, the normal stress is set at integer grid points, and the shear stress and particle vibration velocity components are set at half grid points. The top surface of the dam and the upstream and downstream slopes of the dam are set as free boundaries, and the free boundaries are set at half grid points. The area above the free boundaries is set as vacuum. The bedrock around the dam and the bottom interface are set as absorbing boundaries.

[0036] In one example, a spatial sixth-order, temporal second-order staggered-grid finite-difference method was used to simulate the seismic wavefield along the cross-section of the dam body for the initial model. The normal stress wavefields of the horizontal and vertical components were obtained. t xx , t zz With shear stress t xz Wavefield and initial seismic record u In the differential calculation process, the normal stress is... txx , t zz Set at integer grid points, shear stress t xz With particle vibration velocity components v x , v z The boundary is set at a half-grid point. Based on the actual situation of the dam's cross-sectional structure, the top surface and the upstream and downstream slopes are set as free boundaries, the front and rear surfaces of the bedrock and the bottom interface are set as absorbing boundaries, the free boundary between the dam crest and the dam slope is set at a half-grid point, and the area above the free boundary is set as a vacuum.

[0037] S4. Construct an objective function based on the observed seismic records and the simulated seismic records, and calculate the current objective function value; In this step, the objective function is the standardized L2 norm based on the waveform amplitude, which is defined as:

[0038] in, u i,j Indicates the first i Cannon j Simulated seismic records from receivers d i,j Indicates the first i Cannon j Seismic records observed at receiver points n.s. For the total number of cannons, nr The total number of detector points, E The objective function value represents the normalized error between observed and simulated seismic records.

[0039] S5. Determine whether the current objective function value is less than the preset threshold; if yes, output the current model as the inversion result; if no, proceed to step S6. In one example, given a preset inversion threshold e If the calculated objective function value E Less than e If the initial model is the dam body physical parameter model required for inversion, the inversion process terminates. If the objective function value E Greater than e If the result is positive, proceed to step S6 to continue the inversion iterative calculation.

[0040] S6. Based on the residuals between the observed seismic records and the simulated seismic records, the wavefield is inversely derived to obtain the inverse wavefield, and the gradient of the physical property parameters of the current model is calculated in combination with the forward wavefield. This step specifically includes: First, the residuals between observed and simulated seismic records are used as the excitation source to perform wavefield inverse calculation, thereby obtaining the inverse normal stress wavefield. or xx , or zz With shear stress wave field or xz ; Then, the forward wavefield and the inverse wavefield are cross-correlated to obtain the gradient of the Lamé parameters:

[0041] in, l , m Let Lamé constant be . t xx , t zz , t xz These represent the horizontal normal stress wave field, the vertical normal stress wave field, and the shear stress wave field in the forward modeling simulation. or xx , or xx , or xz These are the horizontal normal stress wave field, the vertical normal stress wave field, and the shear stress wave field, respectively, derived in reverse. g λ Lamé constant l gradient, g μ Lamé constant m gradient, sources Indicates all firing points; Then, the gradients of the Lamé parameters are converted into gradients of P-wave velocity, S-wave velocity, and density using the chain rule, thus obtaining the gradients of the physical property parameters of the current model:

[0042] in, V p For the longitudinal wave velocity, V s For transverse wave velocity, r For density, x , z These represent the horizontal and vertical displacement fields derived in reverse, respectively. v x , v z These are the horizontal and vertical velocity fields of the particle vibration, respectively. g Vp The gradient of the longitudinal wave velocity. gVs The gradient of the transverse wave velocity. g ρ The gradient is the density.

[0043] S7. Use the conjugate gradient method to process the gradient of the physical property parameters of the current model to obtain the conjugate gradient; To improve the inversion convergence speed, this step uses the conjugate gradient method to process the gradients of the physical parameters of the current model, obtaining the conjugate gradients, specifically including: When the number of iterations n When = 1, the conjugate gradient δg 1 k = g 1 k ,in k Representing different physical property parameters, k When =1, g 1 1 g represents the first iteration Vp , k When =2, g 1 2 g represents the first iteration Vs , k When =3, g 1 3 g represents the first iteration ρ .

[0044] when n When ≥2, conjugate gradient δg n k The calculation formula is as follows:

[0045] in, g n k Let be the gradient of the physical property parameters in the nth iteration. The gradient of the physical property parameters in the (n-1)th iteration. δg n k The conjugate gradient of the nth iteration is... The conjugate gradient of the (n-1)th iteration. β n The conjugate correction coefficient. d For conjugate identifiers, This represents the square of the L2 norm.

[0046] S8. Construct a preconditional operator for suppressing surface wave interference and accelerating the updating of physical property parameters related to hidden dangers in the deep region of the dam body. Use the preconditional operator to correct the conjugate gradient to obtain the updated gradient. Seismic waves propagate within the dam body in a spherical spread manner. The farther the propagation distance, the greater the energy attenuation of the wave field; therefore, spherical spread compensation is required when calculating the gradient. Furthermore, considering that potential hazards within the dam body mainly occur in the middle or bottom, and rarely appear in shallow depths at the dam crest (e.g., within 1m), this step uses a piecewise function as the preconditioner to accelerate the update of physical parameters in the deep regions of the dam body. Its expression is:

[0047] in, p For preconditioning operators, Z Z is the depth downwards from the dam crest as the origin. gradt1 Z is the first depth threshold. gradt2 The second depth threshold, H Z is the height of the dam. gradt1 <Z gradt2 < H , △ l Δ represents the depth range of the linear transition region. l =Z gradt2 -Z gradt1 , a、a 1 represents the depth adjustment factor.

[0048] In a preferred embodiment, the above formula... a =3.0, depth adjustment factor a 1 can be set according to the maximum height of the dam body, Z gradt1 =1m, Z gradt2 =2m. Rayleigh surface waves have the strongest wavefield energy in seismic records, but they are not the key wavefield for inverting potential hazards in earth-rock dams. Therefore, to avoid the influence of surface waves on gradient calculations, the preconditioner values ​​within 1m of the dam slope depth are set to 0. The preconditioner values ​​gradually increase from the dam crest to the dam base, with identical preconditioner values ​​at the same depth. In regions where the preconditioner value is 0, the calculated gradient value is 0, and the physical parameters in this region are not updated.

[0049] During model update, this preconditioning operator is used. P For conjugate gradients Make corrections to obtain the updated gradient. .

[0050] S9. Determine the optimal inversion step size for updating the current model, and update the current model according to the update gradient and the optimal inversion step size to obtain a new model; In this step, determining the optimal inversion step size for updating the current model specifically includes: using a parabolic fitting method to search for the inversion step size to obtain the optimal inversion step size.

[0051] Specifically, the parabolic fitting method is used to search for the inversion step size, and three test step sizes are calculated. α 1. α 2 and α 3 and the corresponding objective function values E 1, E 2 and E 3. The test step sizes and the objective function need to satisfy the following constraints:

[0052] Since the current n -th model parameters m n are known, the objective function value E 1 can be calculated, and the corresponding α 1 = 0. Then, the search for α 2 and α 3 starts. According to the formula α 2 n+1 = α 2 n / scale, the value of α 2 is continuously changed until E 2 < E 1, n = 0, 1, 2..., N - 1, where N represents the maximum number of search times, scale is the step size update scale, and the initial α 2 0 and scale are given according to experience. Then, according to the formula α 3 n+1 = α 3 n + α 3 n / scale, the value of α 3 is continuously changed until E3 < E2. The initial α 3 0 = α 2. <​​​​​​​​​​​​​​​​​​​​​​​​​​​3. This is a large-step trial model.

[0053] In search results α 1. α2 and α3 and their corresponding E 1. E 2 and E After step 3, by solving the coefficients of the two linear equations using parabolic fitting, the optimal step size can be obtained. α n .

[0054] Then, the calculation formula for updating the current model based on the update gradient and the optimal inversion step size is as follows:

[0055] in, m n For the current model, m n+1 For the updated new model, α n To achieve the optimal inversion step size, p For preconditioning operators, d g n k is the conjugate gradient of the nth iteration.

[0056] For example, using the above formula for the initial model m After updating, a new model is obtained. m 2.

[0057] S10. The new model is used as the new current model, and the process returns to step S3 for iteration until the current objective function value in step S5 is less than the preset threshold, thereby obtaining a physical property parameter model that reflects the hidden dangers of earth-rock dams.

[0058] In one example, based on the updated model m 2. Repeat steps S3-S5 to obtain the new objective function value. E If E is less than e ,but m 2 represents the inverted dam body physical property parameter model. Otherwise, continue updating the model through steps S6-S10. Repeat the above cycle until the value of the objective function E is less than... e The last updated model m n The model for the physical property parameters of the dam body to be inverted.

[0059] The method of the present invention will be further explained and illustrated below through a specific embodiment.

[0060] Example: This embodiment provides a full waveform inversion method for hidden dangers in earth-rock dams based on a pseudo-well observation system, which specifically includes the following steps: Step 1: Establish a physical property model of the earth-rock dam with potential leakage, including realistic P-wave velocity and S-wave velocity models, such as... Figure 3a and Figure 3b As shown in the model, the dam height is 25m, the bedrock depth is 5m, the dam crest width is 6m, the upstream and downstream slope ratio is 1:2, and the potential for leakage is located at the bottom of the dam. The longitudinal wave velocity of the dam body... V p1 =1160m / s, shear wave velocity V s1 =520m / s, density r 1 = 1790 kg / m 3 Longitudinal wave velocity in the leakage area V p2 =1400m / s, shear wave velocity V s2 =250m / s, density r 2 = 1400 kg / m 3 .

[0061] Seismic wavefield simulations were performed on the aforementioned dam model, and the simulations were used as inversion data from observed seismic records. d In the differential calculation process, the normal stress is... t xx , t zz Set at integer grid points, shear stress t xz With particle vibration velocity components v x , v z Set at half-grid points, such as Figure 4 As shown. Based on the actual situation of the dam's cross-sectional structure, the top surface and upstream and downstream slopes are set as free boundaries, the bedrock front and back and bottom interfaces are set as absorbing boundaries, the free boundaries between the dam crest and slopes are set at half-grid points, and the area above the free boundaries is set as a vacuum, such as... Figure 5 As shown.

[0062] The single-shot seismic record of the shot point in the forward modeling is located at a horizontal distance of 34m, as shown below. Figure 6a and Figure 6b As shown, where Figure 6a For VX component seismic records, Figure 6bThe figure shows the VZ component seismic record. R represents Rayleigh surface waves, DS represents excited direct shear waves, DP represents excited direct P-waves, PS represents converted shear waves at the bedrock interface, SP represents converted P-waves at the bedrock interface, PP represents reflected P-waves at the bedrock interface, SS represents reflected shear waves at the bedrock interface, USS represents multiple reflections of the reflected shear waves at the bedrock interface on the upstream dam slope, DSS represents multiple reflections of the reflected shear waves at the bedrock interface on the downstream dam slope, and MP represents multiple reflected P-waves generated by surface waves at the downstream dam slope angle.

[0063] Grid spacing Δ in forward modeling x =Δz=0.1m, the source is a Ricker wavelet, the dominant frequency is f m =100Hz, sampling time interval dt=2×10 -5 The bedrock absorbing boundary thickness is 40 grid points. Geophone points are arranged horizontally at 0.5m intervals along the dam surface, for a total of 212 geophones. Shot points are arranged horizontally at 5m intervals along the dam surface, firing sequentially from the upstream slope to the downstream slope, for a total of 22 shot points. Figure 6a and Figure 6b Seismic records of different components can identify not only direct waves, reflected and converted waves from the bedrock interface, and reflected shear waves from the seepage area, but also multiple reflected waves generated by the reflected shear waves from the bedrock interface on the upstream and downstream dam slopes, as well as multiple reflected P waves generated by surface waves at the downstream slope angle. Since the difference in shear wave velocity between the dam body and the bedrock is greater, the energy of the shear wave field is stronger than that of the P wave field.

[0064] Step 2: Establish an initial dam body physical property parameter model free of hidden dangers. m 1, such as Figure 7 As shown. Initial model physical property values ​​(longitudinal wave velocity) V p transverse wave velocity V s ,density r The geometry is the same as that in Figure 3.

[0065] Step 3: Using the spatial sixth-order and temporal second-order staggered-grid finite difference method, a seismic wavefield simulation along the cross section of the dam body is performed on the initial model. The normal stresses of the horizontal and vertical components are obtained. t xx , t zz With shear stress t xz Wavefield and initial seismic record u .

[0066] Step 5: Construct the objective function using the normalized L2 norm based on waveform amplitude, as defined below:

[0067] u i,j and d i,j They represent the first i Cannon j Simulated and observed seismic records at receiver points. Given an inversion threshold. e =0.05, if the calculated objective function value E Less than e If the initial model is the dam body physical parameter model required for inversion, the inversion process terminates. If the objective function value E Greater than e If so, then continue with the inversion iterative calculation.

[0068] Step 4: Calculate the residuals of the observed seismic records and the simulated seismic records. The residual record for a single shot located at a horizontal distance of 34m is as follows: Figure 8a and Figure 8b As shown, where Figure 8a For VX component seismic records, Figure 8b This is a seismic record of the vz component. The residual is used as the excitation point for wavefield inverse calculation, yielding the inverse normal stress wavefield. or xx , or zz With shear stress wave field or xz Then, the gradient of the Lamé parameters is calculated by cross-correlation between the forward and inverse wavefields:

[0069] In the formula, l , m Let be Lamé's constant.

[0070] Step 5: Use the chain rule to convert the gradients of the Lamé parameters into gradients of P-wave velocity, S-wave velocity, and density:

[0071] In the formula, v x , v z The velocity component of the particle vibration. x , z This represents the displacement calculated by inverse wave field calculation.

[0072] Step 6: To improve the inversion convergence speed, conjugate gradients are used to update the model. When the number of iterations... n When = 1, the conjugate gradient δg 1k = g 1 k , k These represent different physical property parameters. When n When ≥2, conjugate gradient δg n k The calculation formula is as follows:

[0073] Step 7: Seismic waves propagate within the dam body in a spherical spread manner. The farther the propagation distance, the greater the energy attenuation of the wave field. Therefore, spherical spread compensation is required when calculating the gradient. Furthermore, considering that potential hazards within the dam body mainly occur in the middle or bottom, and rarely appear within a depth of 1m at the dam crest, this paper adopts the following precondition operator to accelerate the update speed of physical property parameters in the deep region of the dam body. P :

[0074] In the formula, a =3.0, a 1 = 0.05, Z gradt1 =1m, Z gradt2 =2m, △ l =Z gradt2 -Z gradt1 Dam height H =25m. The Rayleigh surface wave has the strongest wavefield energy in the seismic record, but it is not the key wavefield for inverting potential hazards in earth-rock dams. To avoid the influence of surface waves on gradient calculations, the values ​​of the precondition operators within a dam slope depth of 1m are set to 0. The constructed inversion precondition operators are as follows: Figure 9 As shown, the values ​​of the preconditioning operators gradually increase from the dam crest to the dam base, and the preconditioning operators at the same depth are identical. In regions where the preconditioning operator is 0, the calculated gradient value is 0, and the physical parameters in this region are not updated.

[0075] Step 8: Use the parabolic fitting method to perform inversion step size search and calculate three test step sizes. α 1. α 2 and α 3 and the corresponding objective function value E 1. E 2 and E 3. The test step size and objective function need to meet the following constraints:

[0076] Due to the current number n Sub-model parameters m n Given the objective function value EIt can be calculated that the corresponding α 1 = 0. Then start searching α 2 and α 3. According to the formula α 2 n+1 = α 2 n / scale, continuously change α the value of 2 until E 2 < E 1, n = 0, 1, 2..., N - 1, N represents the maximum number of search times, scale is the step-size update scale, the initial α 2 0 and scale are given according to experience, α 2 0 = 0.02, β = 5, N = 5. Then according to the formula α 3 n+1 = α 3 n + α 3 n / scale, continuously change α the value of 3 until E3 < E2. The initial α 3 0 = α 2.

[0077] Search for α 1, α2 and α3 and the corresponding E 1, E 2 and E 3. After that, by fitting a parabola to solve the coefficients of the binary linear equation, the optimal step size α n can be obtained.

[0078] Step 9: Use the calculated gradient and inversion step size above to update the initial physical property parameter model to obtain a new model m 2. The update calculation is as follows:

[0079] Step 10: Based on the updated model m 2, repeat steps 3 - 4 to obtain a new objective function value E . If E is less than e , then m 2 is the inverted dam body physical property parameter model. Otherwise, continue to update the model through steps 5 - 11. Repeat the above loop until the value of the objective function E is less than e , then the last updated model m n is the inverted dam body physical property parameter model.

[0080] After 42 iterations, the value of the objective function E is less than... e The inversion model of the longitudinal wave velocity and transverse wave velocity of the dam body is as follows: Figure 10a and Figure 10b As shown, where Figure 10a For the P-wave velocity model, Figure 10b This is a transverse wave velocity model. (Comparison) Figure 3a and Figure 3b The actual shear and p-wave velocity models show that the inversion convergence is good, accurately reproducing the spatial range and geometry of the seepage channels inside the dam. The convergence of the shear wave inversion is better than that of the p-wave inversion. The horizontal profile curves at Z=21.6 of the inverted shear and p-wave velocity models are shown below. Figure 11a and Figure 11b As shown, where Figure 11a This is the longitudinal wave velocity profile. Figure 11b This is a shear wave velocity profile. As can be seen from the figure, the location and velocity of the leakage area in the inverted shear and p-wave models are basically consistent with the actual model. Because the shear wave field in the seismic record has stronger energy and contains more information about the hidden danger area, the convergence speed of the shear wave velocity inversion is faster than that of the p-wave, and the inverted velocity values ​​are also more accurate.

[0081] It should be noted that, to verify the effectiveness of the method of the present invention, step S1 of this embodiment simulates the data acquisition process of a real earth-rock dam through numerical simulation. Specifically, firstly, a physical property parameter model of an earth-rock dam with potential leakage is established, which represents the actual dam body to be investigated; then, a forward modeling simulation is performed on the model, and the obtained simulated seismic records are used as the observed seismic records required for inversion. d This simulates the real data acquisition process. In practical applications, step S1 can also obtain observed seismic records by acquiring seismic data from a real earth-rock dam, which is the same principle as using synthetic data in this embodiment. Based on the verification results of the embodiment, those skilled in the art can reasonably expect that the method of the present invention can also effectively detect hidden dangers in earth-rock dams in practical applications.

[0082] The above examples demonstrate that in the detection of potential hazards in earth-rock dams, the full-waveform inversion method for earth-rock dams based on a pseudo-well observation system proposed in this invention can fully utilize the geometric advantages of the dam body. By acquiring seismic data along the dam's cross-section, it can directly obtain not only direct shear and p-waves penetrating the dam body, reflected and converted waves from the bedrock interface, and reflected shear waves from seepage areas, but also multiple reflected shear waves generated by the reflected shear waves from the bedrock interface on the upstream and downstream dam slopes, and multiple reflected p-waves generated by surface waves at the downstream slope angle. The acquired seismic records contain more information about the hazard area. Therefore, the method proposed in this invention can effectively improve the accuracy of potential hazard detection in earth-rock dams. The gradient of the inversion determines the magnitude and direction of the model update, significantly impacting the accuracy and efficiency of the inversion. The preconditioning operator constructed in this invention can effectively suppress the interference of surface waves on the inversion, highlighting the update speed of the model in the deep regions of the dam body, and improving the efficiency and accuracy of potential hazard inversion in earth-rock dams.

[0083] The various embodiments of the present invention have been described above. These descriptions are exemplary and not exhaustive, nor are they limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments.

Claims

1. A method for full waveform inversion of hidden dangers in earth-rock dams based on a pseudo-well observation system, characterized in that, include: S1. Along the cross section of the earth-rock dam, excitation points and detector points are set up on the dam crest and the upstream and downstream slope surfaces to construct a pseudo-well observation system and collect data to obtain observed seismic records. S2. Establish an initial physical property parameter model of the dam body without hidden dangers, and use this initial model as the current model; S3. Perform forward modeling on the current model to obtain simulated seismic records and forward wavefields; S4. Construct an objective function based on the observed seismic records and the simulated seismic records, and calculate the current objective function value; S5. Determine whether the current objective function value is less than the preset threshold; if yes, output the current model as the inversion result; if no, proceed to step S6. S6. Based on the residuals between the observed seismic records and the simulated seismic records, the wavefield is inversely derived to obtain the inverse wavefield, and the gradient of the physical property parameters of the current model is calculated in combination with the forward wavefield. S7. Use the conjugate gradient method to process the gradient of the physical property parameters of the current model to obtain the conjugate gradient; S8. Construct a preconditional operator for suppressing surface wave interference and accelerating the updating of physical property parameters related to hidden dangers in the deep region of the dam body. Use the preconditional operator to correct the conjugate gradient to obtain the updated gradient. S9. Determine the optimal inversion step size for updating the current model, and update the current model according to the update gradient and the optimal inversion step size to obtain a new model; S10. The new model is used as the new current model, and the process returns to step S3 for iteration until the current objective function value in step S5 is less than the preset threshold, thereby obtaining a physical property parameter model that reflects the hidden dangers of earth-rock dams.

2. The method according to claim 1, characterized in that, In step S1, the simulated well observation system uses the upstream and downstream slopes of the dam as natural excitation and receiving wells. Excitation points and receiver points are arranged along the surface of the dam at preset intervals. The source uses Ricker wavelet with a preset dominant frequency to collect direct longitudinal and transverse waves penetrating the dam, reflected and converted waves from the bedrock interface, reflected waves from the seepage area, and multiple reflected waves generated by the free boundary between the dam slope and the dam crest.

3. The method according to claim 1, characterized in that, In step S2, the geometry of the initial physical property parameter model is consistent with that of the earth-rock dam being probed.

4. The method according to claim 1, characterized in that, In step S3, the forward modeling uses the staggered grid finite difference method to simulate the seismic wavefield; In the differential calculation process, the normal stress is set at integer grid points, and the shear stress and particle vibration velocity components are set at half grid points; The top surface of the dam body and the upstream and downstream slopes are set as free boundaries, which are set at half grid points, and the area above the free boundaries is set as vacuum; the bedrock around the dam body and the bottom interface are set as absorbing boundaries.

5. The method according to claim 1, characterized in that, In step S4, the objective function is the normalized L2 norm based on the waveform amplitude, which is defined as: in, u i,j Indicates the first i Cannon j Simulated seismic records from receivers d i,j Indicates the first i Cannon j Seismic records observed at receiver points ns For the total number of cannons, nr The total number of detector points, E The objective function value represents the normalized error between observed and simulated seismic records.

6. The method according to claim 1, characterized in that, In step S6, calculating the gradient of the physical property parameters of the current model specifically includes: By using the residual between observed and simulated seismic records as the excitation source, the wave field is inversely derived to obtain the inverse normal stress wave field and shear stress wave field. The gradient of the Lamé parameters is obtained by cross-correlation calculation of the forward-modeled wavefield and the inverse-modeled wavefield: in, λ , μ Let Lamé constant be . τ xx , τ zz , τ xz These represent the horizontal normal stress wave field, the vertical normal stress wave field, and the shear stress wave field in the forward modeling simulation. η xx , η xx , η xz These are the horizontal normal stress wave field, the vertical normal stress wave field, and the shear stress wave field, respectively, derived in reverse. g λ Lamé constant λ gradient, g μ Lamé constant μ gradient, sources Indicates all firing points; The gradients of the Lamé parameters are converted into gradients of P-wave velocity, S-wave velocity, and density using the chain rule, thus obtaining the gradients of the physical property parameters of the current model: in, V p For the longitudinal wave velocity, V s For transverse wave velocity, ρ For density, x , z These represent the horizontal and vertical displacement fields derived in reverse, respectively. v x , v z These are the horizontal and vertical velocity fields of the particle vibration, respectively. g Vp The gradient of the longitudinal wave velocity. g Vs The gradient of the transverse wave velocity. g ρ The gradient is the density.

7. The method according to claim 6, characterized in that, In step S7, the conjugate gradient method is used to process the gradients of the physical property parameters of the current model to obtain the conjugate gradients, specifically including: When the number of iterations n When = 1, the conjugate gradient δg 1 k = g 1 k ,in k Representing different physical property parameters, k When =1, g 1 1 g represents the first iteration Vp , k When =2, g 1 2 g represents the first iteration Vs , k When =3, g 1 3 g represents the first iteration ρ ; when n When ≥2, conjugate gradient δg n k The calculation formula is as follows: in, g n k Let be the gradient of the physical property parameters in the nth iteration. The gradient of the physical property parameters in the (n-1)th iteration. δg n k The conjugate gradient of the nth iteration is... The conjugate gradient of the (n-1)th iteration. β n The conjugate correction coefficient. δ For conjugate identifiers, This represents the square of the L2 norm.

8. The method according to claim 1, characterized in that, In step S8, the preconditioning operator is a piecewise function, and its expression is: in, p For preconditioning operators, Z Z is the depth downwards from the dam crest as the origin. gradt1 Z is the first depth threshold. gradt2 The second depth threshold, H Z is the height of the dam. gradt1 <Z gradt2 < H , △ l Δ represents the depth range of the linear transition region. l =Z gradt2 -Z gradt1 , a、a 1 represents the depth adjustment factor.

9. The method according to claim 8, characterized in that, In step S9, determining the optimal inversion step size for updating the current model specifically includes: The optimal inversion step size is obtained by using the parabolic fitting method to search for the inversion step size.

10. The method according to claim 9, characterized in that, In step S9, the calculation formula for updating the current model based on the update gradient and the optimal inversion step size is as follows: in, m n For the current model, m n+1 For the updated new model, α n To achieve the optimal inversion step size, p For preconditioning operators, δg n k is the conjugate gradient of the nth iteration.

Images