Full-wavefield forward method and system in fluid-solid layered medium based on spectral method

Abstract

The present application relates to a kind of full wave field forward method in fluid-solid layered medium based on spectral method, comprising the following steps: S1. establish the acoustic-elastic normal mode model in two-dimensional axisymmetric environment under column coordinate system;S2. using Chebyshev arrangement point method is discretized;S3. the discrete mode matrix of each medium layer is assembled to obtain global matrix, and using constraint condition is replaced by row;S4. the global matrix after completing constraint condition replacement is reconstructed;S5. after correctly normalizing full-depth mode function, it is superimposed to obtain displacement potential field, calculates fluid sound pressure field and solid displacement or stress field.The present application also relates to a kind of full wave field forward method in fluid-solid layered medium based on spectral method system.The present application establishes the joint propagation model of underwater sound and ground sound in fluid-solid layered medium, expands the applicable scope of spectral method, solves the problem that existing technology cannot calculate seabed elastic wave field, fills the existing technical blank.

Description

Technical Field

[0001] This invention relates to the field of marine acoustic modeling technology, and more specifically to a full-wavefield forward modeling method for fluid-solid layered media based on spectral methods. Background Technology

[0002] Oceans cover more than 70% of the Earth's surface, and both the water and seabed contain abundant resources. Marine acoustics, as an important component of marine science, mainly studies the structure and evolution of the ocean and seabed through mathematical and physical methods, numerical simulations, and instrumental observations. Sound waves radiated from underwater sound sources propagate through the water and simultaneously excite elastic wave fields in the layered elastic seabed, forming a full wavefield that couples the underwater acoustic field with the ground acoustic field. However, due to the complex and variable marine environment, and the complex elastic wave equations and interface conditions involved in coupling with the elastic seabed, it is difficult to obtain an analytical solution for the full wavefield in practice. Only approximate solutions can be found through numerical calculations. Therefore, developing theoretical models and numerical simulation algorithms for sound propagation in acoustic-structure coupled systems is of practical significance and provides important support for marine environmental monitoring and seabed resource exploration.

[0003] The shortcomings of existing technology are: 1. Traditional techniques can be divided into two categories based on their different solution strategies: 1.1. Direct solution, that is, using numerical discretization methods to directly discretize and solve the fluid wave equation, the solid elastic wave equation and the corresponding boundary and interface conditions, has the advantages of high computational accuracy, but the computational cost is high and it is difficult to meet the needs of efficient and stable calculation.

[0004] 1.2 Simplified solution: The original governing equations are simplified into another differential equation that is easier to solve using mathematical physics methods, and then numerical methods are used to solve them discretely. This method has advantages such as high computational efficiency and ease of implementation, but it will result in a loss of computational accuracy.

[0005] 2. Among the numerical discretization methods used for solving problems, the finite difference method is widely used due to its convenience and practicality. However, it also suffers from problems such as low computational accuracy, difficulty in constructing high-precision difference schemes, difficulty in mesh generation to account for different trends in physical quantities, and susceptibility to the influence of acoustic or elastic wavelengths. In contrast, the spectral method is a niche numerical discretization method, but its applications in atmospheric acoustics, fluid mechanics, and computational ocean acoustics have demonstrated its advantages such as high computational accuracy and fast error convergence. However, these applications are mainly concentrated in fluid media and cannot meet the computational needs of real marine environments with layered elastic seabeds. Summary of the Invention

[0006] To address the aforementioned problems, this invention provides a full-wavefield forward modeling method and system for fluid-solid layered media based on spectral methods. The purpose is to establish a joint propagation model of underwater acoustics and ground acoustics in fluid-solid layered media, expand the applicability of spectral methods, solve the problem that existing technologies cannot calculate the elastic wavefield of the seabed, and fill the existing technological gap.

[0007] To solve the above problems, the technical solution provided by the present invention is as follows: A full-wavefield forward modeling method for fluid-solid layered media based on spectral methods includes the following steps: S1. Establish an acoustoelastic normal wave model in a two-dimensional axisymmetric environment in cylindrical coordinates, and obtain the modal equations, constraints and normal wave form solutions; S2. Discretize the modal equations and constraints using the Chebyshev collocation method; S3. Assemble the discrete mode matrices of each medium layer to obtain a global matrix, and perform row replacement on the global matrix using the constraints; S4. Reconstruct the global matrix after the constraint replacement is completed, and establish an eigenvalue problem to solve for the horizontal wavenumber and full-depth mode function; S5. After correctly normalizing the full-depth modal functions, the displacement potential field is obtained by superimposing them. The fluid displacement or sound pressure field and the solid displacement or stress field are calculated based on the relationship between the displacement potential, displacement, and stress.

[0008] Preferably, the construction of the modal equations in step S1 includes the following steps: Sa1. The elastic wave equation in solids is separated into longitudinal wave equation and transverse wave equation by Helmholtz decomposition, and then the corresponding longitudinal wave Helmholtz equation and transverse wave Helmholtz equation are obtained by Fourier transform. Sa2. The fluid layer control equations are obtained through the longitudinal wave Helmholtz equations; Sa3. Perform Hankel transformation on the longitudinal wave Helmholtz equation, the transverse wave Helmholtz equation, and the fluid layer control equation to obtain the corresponding modal equations and Sommerfeld integrals; The Helmholtz equation for the longitudinal wave is expressed as follows: ; The Helmholtz equation for the transverse wave is expressed as follows: ; in: and Used to characterize the first i Displacement potential of longitudinal and transverse waves in the layer; and Used to characterize the first i The wave numbers of the longitudinal and transverse waves of the layer, and , , , ; and Used to characterize the first i The attenuation coefficients of longitudinal and transverse waves of the layer, wherein the unit of the attenuation coefficient is dB / wavelength; r Used to characterize horizontal distance; z Used to characterize vertical depth; Used to characterize the i The density of the layer medium; The governing equations for the fluid layer are expressed as follows: ; in: Used to characterize the Dirac function; Used to characterize the depth of a sound source located in a fluid layer; The modal equations are expressed as follows: ; The Sommerfeld integral is expressed as follows: ; in: and Used to characterize the first i The modal functions corresponding to the longitudinal wave displacement potential and the transverse wave displacement potential of the layer; Used to characterize the 0th order Bessel function; Used to characterize the horizontal wavenumber.

[0009] Preferably, the constraint conditions include upper boundary conditions, lower boundary pressure relief conditions, lower boundary rigidity conditions, semi-infinite radiation conditions, fluid-solid interface conditions, and solid-solid interface conditions, wherein: The upper boundary condition is expressed by the following formula: The lower boundary pressure release condition is expressed by the following formula: in: Used to characterize the l The depth of the lower interface of the layer; , The Lamé constant is used to characterize the following formula: ; in: and Used to characterize the first i The wave velocities of longitudinal and transverse waves in the layer; The lower boundary rigidity condition is expressed by the following formula: ; The semi-infinite radiation condition is simulated by setting an artificial absorption layer below the depth of interest to simulate a semi-infinite environment. The lower boundary condition of the artificial absorption layer is calculated by the lower boundary pressure release condition and the lower boundary rigidity condition. The fluid-solid interface condition is expressed by the following formula: in: Used to characterize the angular frequency of a time-harmonic point source, the point source is located at a depth on the vertical axis. Place; The solid-solid interface condition is expressed by the following formula: .

[0010] Preferably, the solution in the normal mode form is expressed by the following formula: in: M Used to characterize the maximum modal order of the cut-off mode; Used for characterization The conjugate; Used to characterize the first type of 0th order Hankel function.

[0011] Preferably, in step S2, the discrete mode equations obtained by discretizing the mode equations using the Chebyshev collocation method are expressed as follows: in: Used to characterize the i The truncation order of the layer; and Used to characterize discrete sequences and Formed 3D column vector, , Used to characterize the first i The mode functions of longitudinal and transverse waves in the first layer are... n Function values ​​at each collocation point; ; Used to characterize the product matrix, it is a N A +1 order diagonal matrix, the elements on its diagonal are the known functions. discrete vectors The value of the element at the corresponding position; Used to characterize the derivative matrix, it is a N A +1 order square matrix, its specific form is expressed by the following formula: ; The upper boundary condition is discretized as follows: The lower boundary pressure release condition is discretized as follows: The lower boundary rigidity condition is discretized as follows: The fluid-solid interface condition is discretized as follows: The solid-solid interface conditions are discretized as follows: in: , for 3D row vectors , , where vector Only the first element is 1, and the rest are 0. (Vector) Only the last element is 1, and the rest are 0.

[0012] Preferably, assembling the discrete mode matrices of each medium layer to obtain the global matrix in step S3 includes the following specific steps: Sb1. Merge the two modal function vectors of each layer into one. A column vector of dimension is expressed as follows: Sb2. The discrete mode equations of each layer are combined to obtain discrete equations, expressed as follows: in: The discrete mode matrix is ​​used to characterize the discrete mode matrix. Square array and All Square matrix; Sb3. The discrete equations across all layers in the entire physical space are simultaneously expressed in matrix form, as follows: in: The global matrix used to characterize the global matrix is ​​a 3D square matrix, column vectors Used to characterize the global mode vector; Used to characterize the horizontal wavenumber.

[0013] Preferably, the rule for performing row replacement on the global matrix using the constraints in step S3 is as follows: When the number of media layers is 1 and the number of rows in the matrix is ​​1, the constraint condition is: ; When the number of media layers is 1, and the number of rows in the matrix is... Then the constraint conditions are: ; When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: ; When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: ; When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: ; When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: or ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: or .

[0014] Preferably, the rule for reconstructing the global matrix after the constraint replacement in step S4 is as follows: all condition replacement rows are moved sequentially to the bottom of the global matrix, and the first and last columns of the longitudinal wave matrix and transverse wave matrix of each layer are moved sequentially to the right side of the matrix; simultaneously, the global mode vector and the elements corresponding to the replacement rows in the global matrix are moved sequentially to the bottom of the global mode vector; the reconstruction result is expressed by the following formula: in: It is Square array It is Square matrix; The eigenvalue problem is expressed as follows: ; Solve the eigenvalue problem and rearrange them according to the order before reconstruction to obtain the horizontal wavenumber and the full-depth mode function.

[0015] Preferably, in step S5, based on the horizontal wavenumber and the full-depth mode function obtained in step S4, the full-depth mode function is normalized by normalized integration to obtain the normalized mode function. Then, the longitudinal wave displacement potential and the transverse wave displacement potential are calculated using the normal form solutions of the longitudinal wave displacement potential and the transverse wave displacement potential. Then, the displacement field or the stress field of the solid layer is calculated by the relationship between the solid layer displacement potential function and the field quantity, and the fluid layer displacement field and the sound pressure field are calculated by the relationship between the fluid layer displacement potential function and the field quantity. The normalized integral is expressed as follows: .

[0016] The relationship between the displacement potential function of the solid layer and the field quantity is expressed by the following formula: in: Used to characterize horizontal displacement in solids; Used to characterize vertical displacement in solids; Used to characterize shear stress in solids; Used to characterize normal stress in solids; The relationship between the fluid layer displacement potential function and the field quantity is expressed by the following formula: in: Used to characterize horizontal displacement in a fluid layer; Used to characterize vertical displacement in a fluid layer; Used to characterize the sound pressure in a fluid layer.

[0017] A system based on the aforementioned spectral method-based full-wavefield forward modeling method in fluid-solid layered media includes a parameter reading module, a parameter interpolation module, an input initialization module, a discrete solution module, a mode selection module, a mode normalization module, a field quantity conversion module, a full-wavefield synthesis module, a plotting module, and a computer device, wherein: The parameter reading module is used to read environmental parameters, sound source parameters, receiver parameters, calculation control parameters, and plotting parameters provided by the user; The parameter interpolation module is used to calculate the ChebyshevGaussLobatto node in each depth interval and interpolate to obtain the value of each parameter at the node. The input initialization module is used to calculate the input variables that subsequent modules need to use based on the parameters. The discrete solution module is used to calculate the discrete values ​​of the modal equations and the constraints, assign values ​​according to the reconstruction form to construct the matrix eigenvalue problem, solve it using a standard solver, and reorder the calculation results of the horizontal wavenumber and modal function according to the form before reconstruction. The modal filtering module is used to calculate the phase velocity of each mode and filter the modes according to the input phase velocity range; The modal normalization module is used to calculate the normalized integral and normalize the filtered modes; The field quantity conversion module is used to calculate the corresponding sound pressure, displacement, and stress modes based on the displacement potential modes. The full-wave field synthesis module is used to calculate the full-wave field through modal superposition, and at the same time calculate the propagation loss field corresponding to each field quantity; The plotting module is used to plot modal function results, propagation loss results, and propagation loss curve results at a specific depth within a specific range, according to user requirements. The computer device includes a memory and a processor; the memory is used to store the programs of each module of the system; the processor is used to execute the programs of each module of the system.

[0018] Compared with the prior art, the present invention has the following advantages: 1. Compared with the direct discretization strategy, this invention simplifies the elastic wave equation through Helmholtz decomposition and further simplifies the Helmholtz equation into a modal equation based on normal mode theory, which is more efficient in calculation and can provide additional modal physical information.

[0019] 2. Compared with the finite difference method in the prior art, the present invention has faster computational convergence and higher computational accuracy.

[0020] 3. Compared with the traditional method of using equivalent reflection coefficient to process elastic media, this method strictly considers the control equation and constraint conditions of the elastic layer, and the eigenvalue calculation is more stable and accurate. It can simultaneously obtain the elastic mode and the elastic wave field.

[0021] 4. The full-wave field forward modeling system in fluid-solid layered media proposed in this invention has a modular functional design that facilitates daily management and continuous upgrades of the system. Attached Figure Description

[0022] Figure 1 This is a schematic diagram of a fluid-layered solid environment model according to a specific embodiment of the present invention; Figure 2This is a schematic diagram of marine environmental sound velocity information according to a specific embodiment of the present invention; Figure 3 This is a schematic diagram of the modal functions according to a specific embodiment of the present invention; Figure 4 This is a schematic diagram of the full-wave field propagation loss results in a specific embodiment of the present invention; Figure 5 This is a schematic diagram of the sound pressure propagation loss curve at a depth of 20m in a specific embodiment of the present invention; Figure 6 This is a schematic diagram of the normal stress propagation loss curve at a depth of 140m, according to a specific embodiment of the present invention. Detailed Implementation

[0023] The present invention will be further illustrated below with reference to specific embodiments. It should be understood that these embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. After reading the present invention, any modifications of the present invention in various equivalent forms by those skilled in the art will fall within the scope defined by the appended claims.

[0024] This invention application claims a full-wavefield forward modeling method for fluid-solid layered media based on spectral methods, comprising the following steps: S1. Establish an acoustoelastic normal wave model in a two-dimensional axisymmetric environment in cylindrical coordinates, and obtain the modal equations, constraints, and normal wave form solutions.

[0025] S2. Discretize the modal equations and constraints using the Chebyshev collocation method.

[0026] S3. Assemble the discrete mode matrices of each medium layer to obtain the global matrix, and perform row substitution on the global matrix using constraints.

[0027] S4. Reconstruct the global matrix after completing the constraint replacement, and establish an eigenvalue problem to solve for the horizontal wavenumber and full-depth mode function.

[0028] S5. After correctly normalizing the full-depth modal functions, the displacement potential field is obtained by superimposing them. The fluid displacement or sound pressure field and the solid displacement or stress field are calculated based on the relationship between the displacement potential, displacement, and stress.

[0029] It should be noted that, as Figure 1 The environmental model shown, the construction of modal equations in step S1 includes the following steps: Sa1. The elastic wave equation in a solid is separated into longitudinal wave equation and transverse wave equation by Helmholtz decomposition, and then the corresponding longitudinal wave Helmholtz equation and transverse wave Helmholtz equation are obtained by Fourier transform.

[0030] Sa2. The governing equations for the fluid layer are obtained through the longitudinal wave Helmholtz equation.

[0031] Sa3. The Helmholtz equations for longitudinal waves, transverse waves, and fluid layer control equations are subjected to Hankel transformation to obtain the corresponding modal equations and Sommerfeld integrals.

[0032] The Helmholtz equation for longitudinal waves is expressed as Equation 1: (1).

[0033] The Helmholtz equation for transverse waves is expressed as Equation 2: (2).

[0034] in: and Used to characterize the first i Displacement potential of longitudinal and transverse waves in the layer; and Used to characterize the first i The wave numbers of the longitudinal and transverse waves of the layer, and , , , ; and Used to characterize the first i The attenuation coefficients of longitudinal and transverse waves in the layer, with the unit of attenuation coefficient being dB / wavelength; r Used to characterize horizontal distance; z Used to characterize vertical depth; Used to characterize the i The density of the layer medium.

[0035] Sound propagation in the fluid layer is governed only by the longitudinal wave Helmholtz equation. Due to the presence of sound sources, the governing equation for the fluid layer is expressed as Equation 3: (3).

[0036] in: Used to characterize the Dirac function; Used to characterize the depth of a sound source located in a fluid layer.

[0037] The modal equations are expressed as equations 4a, 4b, and 4c: (4a) (4b) (4c).

[0038] Sommerfeld's integral expression is expressed as follows: 5a, 5b, 5c: (5a) (5b) (5c).

[0039] in: and Used to characterize the first i The mode functions corresponding to the longitudinal wave displacement potential and the transverse wave displacement potential of the layer; Used to characterize the 0th order Bessel function; Used to characterize the horizontal wavenumber.

[0040] It should be noted that the constraints include upper boundary conditions, lower boundary pressure relief conditions, lower boundary rigidity conditions, semi-infinite radiation conditions, fluid-solid interface conditions, and solid-solid interface conditions, among which: The upper boundary condition is expressed as Equation 6: (6) The conditions for releasing the lower boundary pressure are expressed by equations 7a and 7b: (7a) (7b) in: Used to characterize the l The depth of the lower interface of the layer; , Used to characterize Lamé constant , .

[0041] in: and Used to characterize the first i The wave velocities of longitudinal and transverse waves in the layer.

[0042] The lower boundary rigidity condition is expressed by equations 8a and 8b: (8a) (8b); The infinite radiation condition can be accurately calculated by discretizing the equation into a quadratic eigenvalue problem, but the solution faces problems such as high complexity and unstable solutions. To avoid this problem, this patent simulates a semi-infinite environment by setting an artificial absorption layer below the depth of interest. The lower boundary condition of the artificial absorption layer is calculated using the lower boundary pressure release condition and the lower boundary rigidity condition.

[0043] The fluid-solid interface conditions are expressed by equations 9a, 9b, and 9c: (9a) (9b) (9c) in: Used to characterize the angular frequency of a time-harmonic point source, the point source is located at a depth on the vertical axis. Place.

[0044] The solid-solid interface conditions are expressed by equations 10a, 10b, 10c, and 10d: (10a) (10b) (10c) (10d).

[0045] It should be noted that the normal mode wave formal solution is expressed according to equations 11 and 12: (11) (12) in: M Used to characterize the maximum modal order of the cut-off mode; Used for characterization The conjugate; Used to characterize the first type of 0th order Hankel function.

[0046] It should be noted that the Chebyshev collocation method uses Chebyshev polynomials as basis functions. The spectral method uses the basis functions and weight functions as weight functions. The specific forms of the basis functions and weight functions are expressed according to Equation 13: (13) Discrete points in the weighting function These are called collocations, taken as ChebyshevGaussLobatto (CGL) nodes, and expressed in Equation 14: (14) Assume the truncation order is N The expansion coefficient is Then, using basis functions can transform the function Approximately expressed as Equation 15: (15) As can be seen from the principle of the weighted residual method, the Chebyshev collocation method only requires that the equation be in this... The collocation holds strictly on a given set of points. That is, when discretizing a function using the Chebyshev collocation method, it is equivalent to representing the function using the values ​​of the function at a set of CGL nodes. For example, a function can be used... exist N +1 node Discrete values ​​on Composition N +1 dimension column vector To represent a function.

[0047] m The relationship between the second derivative and the antiderivative is: ,in The derivative matrix is ​​a N +1 order square matrix.

[0048] It should be further explained that the discrete mode equations obtained by discretizing the mode equations using the Chebyshev collocation method in step S2 are expressed as Equations 16 and 17: (16) (17) in: Used to characterize the i The truncation order of the layer; and Used to characterize discrete sequences and Formed 3D column vector, , Used to characterize the first i The mode functions of longitudinal and transverse waves in the first layer are... n Function values ​​at each collocation point; Suppose a known function and the function to be solved The product is Then the product term discrete vectors satisfy , Used to characterize the product matrix, it is a N A +1 order diagonal matrix, the elements on its diagonal are the known functions. discrete vectors The value of the element at the corresponding position; Used to characterize the derivative matrix, it is a N A +1 order square matrix, its specific form is expressed as Equation 18: (18).

[0049] It should be further noted that when using CGL nodes for discretization, the integral of the function to be solved can be approximately calculated using Equation 19: (19) The modal equations are discretized using the Chebyshev collocation method, and coordinate transformation is introduced. Transform the domain of the original depth to the domain allowed by the Chebyshev polynomial. Above. After coordinate transformation, the Helmholtz equations for longitudinal and transverse waves become equations 20 and 21: (20) (twenty one) use Indicates the first i The truncation order of the layer, at the th i Take from layer Use CGL nodes as collocation points to discretize the above equation.

[0050] It should be further explained that the upper boundary condition is discretized as Equation 22: (twenty two) The lower boundary pressure release condition is discretized into equations 23a and 23b: (23a) (23b) The lower boundary rigidity condition is discretized into equations 24a and 24b: (24a) (24b) The fluid-solid interface conditions are discretized as equations 25a, 25b, and 25c: (25a) (25b) (25c) The solid-solid interface conditions are discretized as equations 26a, 26b, 26c, and 26d: (26a) (26b) (26c) (26d) in: , for 3D row vectors , , where vector Only the first element is 1, and the rest are 0. (Vector) Only the last element is 1, and the rest are 0.

[0051] It should be noted that the assembly of the discrete mode matrices of each medium layer to obtain the global matrix in step S3 includes the following specific steps: Sb1. Merge the two modal function vectors of each layer into one. A column vector of dimension 27 is expressed as follows: (27) Sb2. Combine the discrete mode equations of each layer to obtain the discrete equations, expressed as Equation 28: (28) in: Used to characterize the discrete mode matrix, it is Square array and All A square matrix of order.

[0052] Sb3. Represent the discrete equations across all layers in the entire physical space in matrix form. Since there are no transverse waves in the fluid layer, therefore... Discarded, expressed according to equations 29 and 30: (29) (30) in: Used to represent the global matrix, it is 3D square matrix, column vectors Used to characterize the global mode vector; Used to characterize the horizontal wavenumber.

[0053] It should be noted that the rule for performing row replacement on the global matrix using constraints in step S3 is as follows: When the number of media layers is 1 and the number of rows in the matrix is ​​1, the constraint condition is: .

[0054] When the number of media layers is 1, and the number of rows in the matrix is... Then the constraint conditions are: .

[0055] When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: .

[0056] When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: .

[0057] When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: .

[0058] When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: .

[0059] When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: .

[0060] When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: .

[0061] When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: .

[0062] When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: .

[0063] When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: .

[0064] When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: or .

[0065] When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: .

[0066] When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: or .

[0067] It should be noted that the rules for reconstructing the global matrix after constraint replacement in step S4 are as follows: all condition replacement rows are moved sequentially to the bottom of the global matrix, and the first and last columns of the P-wave and S-wave matrices of each layer are moved sequentially to the right side of the matrix; simultaneously, the global mode vectors and the elements corresponding to the replacement rows in the global matrix are moved sequentially to the bottom of the global mode vectors; the reconstructed matrices and vectors are respectively... and The reconstruction result is expressed according to Equation 31: (31) in: It is Square array It is A square matrix of order.

[0068] The eigenvalue problem is expressed by Equation 32: (32).

[0069] Solve the eigenvalue problem and rearrange the values ​​according to their original order to obtain the horizontal wavenumber and full-depth mode functions. Set the phase velocity range. Screen the primary modes.

[0070] It should be noted that in step S5, based on the horizontal wavenumber and full-depth modal function obtained in step S4, the normalized integral given by equation 33 is solved using equation 19. The full-depth modal function is normalized by the normalized integral to obtain the normalized modal function. Then, the longitudinal wave displacement potential and transverse wave displacement potential are calculated using the normal form solutions of the longitudinal wave displacement potential and transverse wave displacement potential. Then, the displacement field or stress field of the solid layer is calculated by the relationship between the solid layer displacement potential function and the field quantity. The displacement field and sound pressure field of the fluid layer are calculated by the relationship between the fluid layer displacement potential function and the field quantity.

[0071] The normalized integral is expressed as follows: (33).

[0072] The relationship between the displacement potential function of the solid layer and the field quantity is expressed by equations 34, 35, and 36: (34) (35) (36) in: Used to characterize horizontal displacement in solids; Used to characterize vertical displacement in solids; Used to characterize shear stress in solids; Used to characterize normal stress in solids.

[0073] The relationship between the fluid layer displacement potential function and the field quantity is expressed by equations 37 and 38: (37) (38) in: Used to characterize horizontal displacement in a fluid layer; Used to characterize vertical displacement in a fluid layer; Used to characterize the sound pressure in a fluid layer.

[0074] Finally, the propagation loss of the field quantity is used to illustrate the wavefield calculation results. For example, the propagation loss TL of the sound pressure is expressed by Equation 39: (39) Where: TL is in decibels (dB); Used to characterize the sound pressure value at a distance of 1m from the sound source.

[0075] A system based on a full-wavefield forward modeling method for fluid-solid layered media using spectral methods includes a parameter reading module, a parameter interpolation module, an input initialization module, a discrete solution module, a mode selection module, a mode normalization module, a field quantity conversion module, a full-wavefield synthesis module, a plotting module, and a computer device, wherein: The parameter reading module is used to read user-provided environmental parameters, sound source parameters, receiver parameters, calculation control parameters, and plotting parameters.

[0076] The parameter interpolation module is used to calculate the ChebyshevGaussLobatto nodes in each depth range and interpolate the values ​​of each parameter at the node.

[0077] The input initialization module is used to calculate the input variables that subsequent modules need to use based on the parameters.

[0078] The discrete solution module is used to calculate the discrete values ​​of the modal equations and constraints, assign values ​​according to the reconstruction form to construct the matrix eigenvalue problem, solve it using the standard solver, and reorder the calculation results of the horizontal wavenumber and modal function according to the form before reconstruction.

[0079] The modal filtering module is used to calculate the phase velocity of each mode and filter the modes according to the input phase velocity range.

[0080] The modal normalization module is used to calculate the normalized integral and normalize the filtered modes.

[0081] The field quantity conversion module is used to calculate the corresponding sound pressure, displacement, and stress modes based on the displacement potential modes.

[0082] The full-wave field synthesis module is used to calculate the full-wave field through modal superposition, and at the same time calculate the propagation loss field corresponding to each field quantity.

[0083] The plotting module is used to plot modal function results, propagation loss results, and propagation loss curve results at a specific depth within a specific range, according to user requirements.

[0084] Computer equipment includes memory and processor; memory is used to store the programs of the various modules of the system; processor is used to execute the programs of the various modules of the system.

[0085] It should be noted that, in order to provide a more detailed explanation of the implementation method, a typical full-wavefield forward modeling example in a shallow sea environment is given.

[0086] Example: Environmental parameters: sound speed information, such as Figure 2 As shown, the seawater depth is 127 m and the density is 1.0 g / cm³. 3 There was no absorption attenuation; the sedimentary layer was 23 m thick and the density ranged from 1.7 g / cm³. 3 It increases linearly with depth to 1.8 g / cm 3 The attenuation coefficients for longitudinal and transverse waves are 0.2 and 0.2, respectively. and 0.4 Below the sedimentary layer is semi-infinite bedrock with a density of 2.2 g / cm³. 3 The attenuation coefficients for longitudinal and transverse waves are 0.1 and 0.1, respectively. and 0.2 .

[0087] Sound source parameters: The sound source is located at a depth of 20 m in the seawater, and the sound source frequency is 100 Hz.

[0088] Receiver parameters: Receivers are deployed at depths of 20 m and 140 m, respectively, at horizontal intervals of 10 m within a horizontal distance range of 0 to 10 km.

[0089] Calculation control parameters: the cutoff order of the seawater layer is 100, and the cutoff order of the sedimentary layer is 20; the semi-infinite bedrock is simulated by setting a 50 m thick bedrock layer and a 200 m thick absorbing layer, with the P-wave and S-wave attenuation coefficients of the absorbing layer both being 5. The remaining parameters are the same as those of the bedrock. The lower boundary of the absorption layer is selected as the pressure release boundary. The cutoff order of the bedrock layer is 30, and the cutoff order of the absorption layer is 20. The phase velocity screening range is 0~4500 m / s. The maximum calculation distance is 10 km, the horizontal interval is 10 m, and the vertical interval is 1 m.

[0090] Plotting parameters: Select to plot the fluid acoustic pressure mode and solid normal stress mode within the depth range of 0~180 m; select to plot the full wave field propagation loss within the depth range of 0~180 m and the distance range of 0~10 km, including the propagation loss of fluid acoustic pressure and solid normal stress; select to plot the propagation loss curves of fluid acoustic pressure and solid normal stress received by the receiver.

[0091] Step 1 In this example, the first layer of seawater is a fluid layer, and the following three layers of elastic seabed are solid layers. According to Helmholtz decomposition, the elastic wave equation in the elastic seabed can be separated into two wave equations, and then the corresponding Helmholtz equations can be obtained through Fourier transform: in, , The environmental parameters have already been provided. and The sound velocity and attenuation coefficient are given by environmental parameters, and the frequency is given by sound source parameters.

[0092] The governing equations for the ocean layer are: in, The environmental parameters have already been provided. The parameters have already been given by the sound source.

[0093] Performing the Hankel transformation on equations (40) to (42) yields the corresponding modal equations. and Sommerfeld integral The Lamé constant is calculated based on the sound velocity and density given by environmental parameters. and In the environment considered in the example, the upper boundary is the seawater surface, which can be regarded as a pressure relief boundary. After the Hankel transformation, it becomes... The lower boundary of the absorption layer is also set as a pressure relief boundary, with the condition being: In the formula, , The control parameters have already been calculated and given.

[0094] The seawater-sediment interface is a fluid-solid interface, and its conditional form is as follows: The sedimentary layer-bedrock interface is a solid-solid interface, and its conditional form is as follows: The normal wave solutions of the displacement potentials of longitudinal and transverse waves are: in, M Determined during the modality screening process, , and This is given from the solution process of the matrix eigenvalue problem.

[0095] The normalized integral is Step 2 Discretize the modal equation in equation (43) into the following matrix form. Discretize the upper boundary condition in equation (45) as follows: Discretize the lower boundary pressure release condition in equation (46) as follows: Discretize the fluid-solid interface condition in equation (47) as follows: Discretize the solid-solid interface condition in equation (48) as follows: Step 3 The discrete equations across all layers in the entire physical space are simultaneously expressed in matrix form. Since there are no transverse waves in the seawater layer, therefore... Discarding, we get the following form Calculate the CGL nodes corresponding to each layer, interpolate to obtain the values ​​of the input parameters at the CGL nodes, form a parameter vector, and then calculate the discrete matrix and discrete constraint vector of each layer according to the given common discrete terms.

[0096] Step 4 Based on the reconstructed matrix form, the global matrix is ​​assigned values ​​using the calculated discrete matrices and discrete constraint vectors of each layer. The horizontal wavenumber and mode functions are calculated using a standard eigenvalue solver and rearranged according to their original order. The main modes are selected based on the phase velocity range, and each mode is normalized. The horizontal wavenumber results are shown in Table 1, and the real part results of the corresponding mode functions are shown in Table 2. Figure 3 As shown.

[0097] Table 1. Partial Results of Horizontal Wavenumber Step 5 The longitudinal and transverse wave displacement potentials are calculated using equations (11) and (12). Based on the relationships given in equations (34) to (36), the displacement or stress field of the solid layer can be calculated from the displacement potential. The displacement and sound pressure fields of the fluid layer can be calculated using equations (37) and (38). Finally, the corresponding propagation loss is calculated. The full-wavelength propagation loss results are as follows: Figure 4 As shown, the seawater layer is described by sound pressure, while the sedimentary layer and bedrock are described by normal stress. The sound pressure propagation loss curve at a depth of 20 m is shown below. Figure 5 The normal stress propagation loss curve at a depth of 140 m is as follows: Figure 6 .

[0098] In the above detailed description, various features are combined together in a single embodiment to simplify this disclosure. This approach to disclosure should not be construed as reflecting an intention that embodiments of the claimed subject matter require more features than are explicitly stated in each claim. Rather, as reflected in the appended claims, the invention is presented with fewer features than all of the features of the single disclosed embodiment. Therefore, the appended claims are hereby explicitly incorporated into the detailed description, wherein each claim stands alone as a preferred embodiment of the invention.

[0099] The disclosed embodiments have been described above to enable any person skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to those skilled in the art, and the general principles defined herein can be applied to other embodiments without departing from the spirit and scope of this disclosure. Therefore, this disclosure is not limited to the embodiments given herein, but is consistent with the broadest scope of the principles and novel features disclosed in this application.

[0100] The foregoing description includes examples of one or more embodiments. It is certainly impossible to describe all possible combinations of components or methods in order to describe the above embodiments, but those skilled in the art will recognize that further combinations and arrangements of the various embodiments are possible. Therefore, the embodiments described herein are intended to cover all such changes, modifications, and variations that fall within the scope of the appended claims. Furthermore, the term "comprising" as used in the specification or claims is interpreted in a manner similar to the term "including," as interpreted when used as a conjunction in the claims. Additionally, the use of any term "or" in the specification of the claims is intended to mean "non-exclusive or."

[0101] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A full-wavefield forward modeling method for fluid-solid layered media based on spectral methods, characterized in that: Includes the following steps: S1. Establish an acoustoelastic normal wave model in a two-dimensional axisymmetric environment in cylindrical coordinates, and obtain the modal equations, constraints and normal wave form solutions; S2. Discretize the modal equations and constraints using the Chebyshev collocation method; S3. Assemble the discrete mode matrices of each medium layer to obtain a global matrix, and perform row replacement on the global matrix using the constraints; S4. Reconstruct the global matrix after the constraint replacement is completed, and establish an eigenvalue problem to solve for the horizontal wavenumber and full-depth mode function; S5. After correctly normalizing the full-depth modal functions, the displacement potential field is obtained by superimposing them. The fluid displacement or sound pressure field and the solid displacement or stress field are calculated based on the relationship between the displacement potential, displacement, and stress.

2. The full-wavefield forward modeling method for fluid-solid layered media based on spectral methods according to claim 1, characterized in that: The construction of the modal equations in step S1 includes the following steps: Sa1. The elastic wave equation in solids is separated into longitudinal wave equation and transverse wave equation by Helmholtz decomposition, and then the corresponding longitudinal wave Helmholtz equation and transverse wave Helmholtz equation are obtained by Fourier transform. Sa2. The fluid layer control equations are obtained through the longitudinal wave Helmholtz equations; Sa3. Perform Hankel transformation on the longitudinal wave Helmholtz equation, the transverse wave Helmholtz equation, and the fluid layer control equation to obtain the corresponding modal equations and Sommerfeld integrals; The Helmholtz equation for the longitudinal wave is expressed as follows: ； The Helmholtz equation for transverse waves is expressed as follows: ； in: and Used to characterize the first i Displacement potential of longitudinal and transverse waves in the layer; and Used to characterize the first i The wave numbers of the longitudinal and transverse waves of the layer, and , , , ; and Used to characterize the first i The attenuation coefficients of longitudinal and transverse waves of the layer, wherein the unit of the attenuation coefficient is dB / wavelength; r Used to characterize horizontal distance; z Used to characterize vertical depth; Used to characterize the i The density of the layer medium; The governing equations for the fluid layer are expressed as follows: ； in: Used to characterize the Dirac function; Used to characterize the depth of a sound source located in a fluid layer; The modal equations are expressed as follows: ； The Sommerfeld integral is expressed as follows: ； in: and Used to characterize the first i The modal functions corresponding to the longitudinal wave displacement potential and the transverse wave displacement potential of the layer; Used to characterize the 0th order Bessel function; Used to characterize the horizontal wavenumber.

3. The full-wavefield forward modeling method for fluid-solid layered media based on spectral methods according to claim 2, characterized in that: The constraints include upper boundary conditions, lower boundary pressure relief conditions, lower boundary rigidity conditions, semi-infinite radiation conditions, fluid-solid interface conditions, and solid-solid interface conditions, wherein: The upper boundary condition is expressed by the following formula: The lower boundary pressure release condition is expressed by the following formula: in: Used to characterize the l The depth of the lower interface of the layer; , The Lamé constant is used to characterize the following formula: ； in: and Used to characterize the first i The wave velocities of longitudinal and transverse waves in the layer; The lower boundary rigidity condition is expressed by the following formula: ； The semi-infinite radiation condition is simulated by setting an artificial absorption layer below the depth of interest to simulate a semi-infinite environment. The lower boundary condition of the artificial absorption layer is calculated by the lower boundary pressure release condition and the lower boundary rigidity condition. The fluid-solid interface condition is expressed by the following formula: in: Used to characterize the angular frequency of a time-harmonic point source, the point source is located at a depth on the vertical axis. Place; The solid-solid interface condition is expressed by the following formula: 。 4. The full-wavefield forward modeling method for fluid-solid layered media based on spectral methods according to claim 3, characterized in that: The solution in the normal mode form is expressed by the following formula: in: M Used to characterize the maximum modal order of the cut-off mode; Used for characterization Conjugate; Used to characterize the first type of 0th order Hankel function.

5. The full-wavefield forward modeling method for fluid-solid layered media based on spectral methods according to claim 4, characterized in that: In step S2, the discretization of the modal equations using the Chebyshev collocation method yields the discrete modal equations, expressed as follows: in: Used to characterize the i The truncation order of the layer; and Used to characterize discrete sequences and Formed 3D column vector, , Used to characterize the first i The mode functions of longitudinal and transverse waves in the first layer are... n Function values ​​at each collocation point; ; Used to characterize the product matrix, it is a N A +1 order diagonal matrix, the elements on its diagonal are the known functions. discrete vectors The value of the element at the corresponding position; Used to characterize the derivative matrix, it is a N A +1 order square matrix, its specific form is expressed by the following formula: ； The upper boundary condition is discretized as follows: The lower boundary pressure release condition is discretized as follows: The lower boundary rigidity condition is discretized as follows: The fluid-solid interface condition is discretized as follows: The solid-solid interface conditions are discretized as follows: in: , for 3D row vectors , , where vector Only the first element is 1, and the rest are 0. (Vector) Only the last element is 1, and the rest are 0.

6. The full-wavefield forward modeling method for fluid-solid layered media based on spectral methods according to claim 5, characterized in that: Step S3, which involves assembling the discrete mode matrices of each dielectric layer to obtain the global matrix, includes the following specific steps: Sb1. Merge the two modal function vectors of each layer into one. A column vector of dimension is expressed as follows: Sb2. The discrete mode equations of each layer are combined to obtain discrete equations, expressed as follows: in: The discrete mode matrix is ​​used to characterize the discrete mode matrix. Square array and All Square matrix; Sb3. The discrete equations across all layers in the entire physical space are simultaneously expressed in matrix form, as follows: in: The global matrix used to characterize the global matrix is ​​a 3D square matrix, column vectors Used to characterize the global mode vector; Used to characterize the horizontal wavenumber.

7. The full-wavefield forward modeling method for fluid-solid layered media based on spectral methods according to claim 6, characterized in that: The rule for performing row replacement on the global matrix using the aforementioned constraints in step S3 is as follows: When the number of media layers is 1 and the number of rows in the matrix is ​​1, the constraint condition is: ; When the number of media layers is 1, and the number of rows in the matrix is... Then the constraint conditions are: ; When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: ; When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: ; When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: ; When the number of media layers is 2, and the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: or ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: ; When the number of dielectric layers is And the number of rows in the matrix is Then the constraint conditions are: or .

8. The full-wavefield forward modeling method for fluid-solid layered media based on spectral methods according to claim 7, characterized in that: The rules for reconstructing the global matrix after the constraint replacement in step S4 are as follows: all condition replacement rows are moved sequentially to the bottom of the global matrix, and the first and last columns of the P-wave and S-wave matrices of each layer are moved sequentially to the right side of the matrix; simultaneously, the global mode vectors and the elements corresponding to the replacement rows in the global matrix are moved sequentially to the bottom of the global mode vectors; the reconstruction result is expressed by the following formula: in: It is Square array It is Square matrix; The eigenvalue problem is expressed as follows: ； Solve the eigenvalue problem and rearrange them according to the order before reconstruction to obtain the horizontal wavenumber and the full-depth mode function.

9. The full-wavefield forward modeling method for fluid-solid layered media based on spectral methods according to claim 8, characterized in that: In step S5, based on the horizontal wavenumber and the full-depth mode function obtained in step S4, the full-depth mode function is normalized by normalized integration to obtain the normalized mode function. Then, the longitudinal wave displacement potential and the transverse wave displacement potential are calculated using the normal form solutions of the longitudinal wave displacement potential and the transverse wave displacement potential. Then, the displacement field or the stress field of the solid layer is calculated by the relationship between the solid layer displacement potential function and the field quantity, and the fluid layer displacement field and the sound pressure field are calculated by the relationship between the fluid layer displacement potential function and the field quantity. The normalized integral is expressed as follows: ； The relationship between the displacement potential function of the solid layer and the field quantity is expressed by the following formula: in: Used to characterize horizontal displacement in solids; Used to characterize vertical displacement in solids; Used to characterize shear stress in solids; Used to characterize normal stress in solids; The relationship between the fluid layer displacement potential function and the field quantity is expressed by the following formula: in: Used to characterize horizontal displacement in a fluid layer; Used to characterize vertical displacement in a fluid layer; Used to characterize the sound pressure in a fluid layer.

10. A system based on the full-wavefield forward modeling method for fluid-solid layered media according to any one of claims 1 to 9, characterized in that: It includes a parameter reading module, a parameter interpolation module, an input initialization module, a discrete solution module, a mode selection module, a mode normalization module, a field quantity conversion module, a full-wave field synthesis module, a plotting module, and computer equipment, among which: The parameter reading module is used to read environmental parameters, sound source parameters, receiver parameters, calculation control parameters, and plotting parameters provided by the user; The parameter interpolation module is used to calculate the ChebyshevGaussLobatto node in each depth interval and interpolate to obtain the value of each parameter at the node. The input initialization module is used to calculate the input variables that subsequent modules need to use based on the parameters. The discrete solution module is used to calculate the discrete values ​​of the modal equations and the constraints, assign values ​​according to the reconstruction form to construct the matrix eigenvalue problem, solve it using a standard solver, and reorder the calculation results of the horizontal wavenumber and modal function according to the form before reconstruction. The modal filtering module is used to calculate the phase velocity of each mode and filter the modes according to the input phase velocity range; The modal normalization module is used to calculate the normalized integral and normalize the filtered modes; The field quantity conversion module is used to calculate the corresponding sound pressure, displacement, and stress modes based on the displacement potential modes. The full-wave field synthesis module is used to calculate the full-wave field through modal superposition, and at the same time calculate the propagation loss field corresponding to each field quantity; The plotting module is used to plot modal function results, propagation loss results, and propagation loss curve results at a specific depth within a specific range, according to user requirements. The computer device includes a memory and a processor; the memory is used to store the programs of each module of the system; the processor is used to execute the programs of each module of the system.

Images