Spherical wave reflection sound field calculation and inversion method, system and apparatus based on global interpolation decoupling

Abstract

The present application relates to a kind of based on global interpolation decoupling spherical wave reflection sound field calculation and inversion method, system and equipment, belong to geophysical and acoustic detection field, the method includes the following steps: input observation data and initialize seabed parameter;Response kernel interpolation node is set in integral interval and interpolation mode is selected;Propagation matrix is calculated by interpolation weight and wave field propagation kernel;Response vector is calculated by response kernel;Spherical wave reflection sound field is calculated using propagation matrix and response vector;According to objective function, update seabed parameter value and condition judgment.The method uses global interpolation to decouple spherical wave reflection sound field, avoids the solution of each order derivative of existing local Taylor expansion method response kernel and corresponding propagation matrix;Avoid redundant calculation, improve the calculation efficiency of spherical wave inversion.

Description

Technical Field

[0001] This invention belongs to the field of geophysics and acoustic detection, and specifically relates to a method, system and equipment for calculating and inverting spherical wave reflection sound field based on global interpolation decoupling. Background Technology

[0002] The seabed reflected sound field carries rich information about the physical properties and stratification of seabed sediments. Currently, the calculation of the seabed reflection coefficient is mainly based on plane wave theory. However, in actual detection, the sound source (seismic source) emits spherical waves, and plane waves are only an approximation of spherical waves in the high-frequency, far-field situation. This makes it difficult for plane wave theory to accurately describe the seabed reflection coefficient in low-frequency, near-field, and large incident angle situations. Spherical wave theory can more accurately characterize the relationship between the seabed reflection coefficient and the incident angle, which can be represented by the Sommerfeld integral (wavenumber integral). However, its integrand contains a kernel function that oscillates with height depending on the transmission and reception distance. Numerical integration of a rapidly oscillating function makes the computational efficiency of spherical wave reflected sound field significantly lower than that of the plane wave method. In the inversion process, each evaluation of the objective function requires expensive integration calculations of the oscillating function. Therefore, inefficient spherical wave forward modeling has become a major bottleneck restricting the inversion process.

[0003] To address the aforementioned issues, one approach to accelerate the process is to separate the response kernel, which varies with medium parameters, from the propagation kernel, which oscillates violently with the transmitting and receiving geometric information, in the integrand of the spherical wave reflection acoustic field integral. This allows the latter to be calculated and stored offline, while online iterations during the inversion process only require inexpensive updates to the former. This decoupling framework effectively improves inversion efficiency. However, current work uses piecewise local Taylor expansions for decoupling, where the response end consists of the derivatives of the plane wave reflection coefficients at various sub-interval locations. Consequently, the offline phase requires pre-calculating and storing multiple sets of propagation kernel matrix information corresponding to the derivatives as the expansion order increases, and each online iteration relies on the repeated calculation of higher-order derivatives of the reflection coefficients. Summary of the Invention

[0004] To address the issues of low computational efficiency of seabed spherical wave reflection coefficients and the reliance of existing local Taylor expansion decoupling methods on high-order derivative information, this invention provides a method for calculating and inverting spherical wave reflection sound fields based on global interpolation decoupling. This method reconstructs plane wave reflection coefficients through global interpolation, decoupling the oscillatory part dependent on the observation geometry from the reflection coefficient part dependent on seabed parameters in the integral expression of the spherical wave reflection sound field. This allows the original integral to be written as the inner product of the propagation matrix and the seabed parameter response vector, avoiding the calculation of high-order derivatives. Furthermore, an efficient seabed parameter inversion method is established within this global interpolation decoupling framework.

[0005] This invention is achieved through the following technical solution: A method for calculating and inverting the acoustic field of spherical wave reflection based on global interpolation decoupling, the method comprising the following steps: Step 1: Input observation data Initialize seafloor strata properties and layering information ; Step 2: Set fixed interpolation nodes within the integration interval and select the interpolation method; the nodes can be equally spaced or non-equally spaced. Furthermore, in the second step, based on the observation data input in the first step... Set fixed response kernel interpolation nodes within the integration interval. Among them, the node spacing Upper limit of integration, whether the spacing is equal or non-equal. Determine using the following formula: ; In the formula, For wave number, for The maximum value of the response kernel within the interval. The total propagation distance of the wave field from the excitation point, through reflection from the seabed, to the receiving point. Representing observation data The sum of the distances from the corresponding excitation point and receiver point to the seabed reflector surface, Indicates the allowable truncation error.

[0006] Step 3: Calculate the propagation matrix using interpolation weights and wavefield propagation kernel; the propagation matrix depends only on the observation geometry, including frequency, transmission distance, and transmission / reception point location, and is independent of seabed parameters. It is pre-calculated and stored once in the offline stage. Furthermore, in the third step, the interpolation weights are selected according to the interpolation method set in the second step. Indicates the domain is The Interpolation weights; the propagation matrix is ​​constructed as follows: , Each element Calculated using the following formula: ; in The wave field propagation kernel determines the properties of wave field propagation; The imaginary unit, The index of the interpolation node. , The number of interpolation nodes set for the second step. The component index in the observed data. (The upper limit of the value is determined by the number of observed data components.) The first in the observation data The horizontal distance from the point source to the receiving point corresponding to each component. Indicates the first in the observation data The sum of the distances from the point source and the receiving point to the seabed corresponding to each component. The first in the observation data The angular frequency corresponding to each component The first in the observation data The vertical wavenumber corresponding to each component The wave number in water. Indicates the horizontal wave number. It is an exponential function. It is a zero-order Bessel function of the first kind.

[0007] Step 4: Calculate the response vector from the response kernel; the response vector contains the function value information of the response kernel at the above fixed interpolation nodes, depending on the seabed parameters; Furthermore, in the fourth step, the response kernel function value is calculated at the fixed interpolation nodes set in the second step. , Obtain the corresponding vector , This is a vector representing the physical properties and layering information of the seafloor strata.

[0008] Step 5: Calculate the decoupled spherical wave reflection sound field by taking the inner product of the propagation matrix obtained in Step 3 and the response vector obtained in Step 4. The calculation expression is as follows: ; matrix It includes information on the spatial propagation of the wave field and the distribution of discrete nodes within the integration interval, but does not include information on seafloor strata properties and layering that determine the plane wave reflection coefficient. ;vector It contains information on seafloor strata properties and layering that determine the plane wave reflection coefficient. ; Step 6: Construct the objective function using the spherical wave reflection sound field and update the seabed parameters based on the objective function; Step 7: Condition judgment: Return to step (4) or output the current seabed parameters.

[0009] The present invention also provides a spherical wave reflection sound field calculation and inversion system based on global interpolation decoupling, wherein the system runs the spherical wave reflection sound field calculation and inversion method based on global interpolation decoupling.

[0010] The present invention also provides a device for calculating and inverting the sound field of spherical wave reflection based on global interpolation decoupling, wherein the device is equipped with the inversion system and executes the inversion method.

[0011] The beneficial effects of this invention compared to the prior art are as follows: This invention addresses the problems of low computational efficiency in existing spherical wave reflection inversion methods and the reliance on high-order derivative information in existing local Taylor expansion decoupling methods. It proposes a global interpolation-based method for calculating and inverting the spherical wave reflection acoustic field: using plane wave reflection coefficients at a finite number of fixed nodes, the plane wave reflection coefficients over the entire integration interval are reconstructed through global interpolation. This decouples the spherical wave reflection acoustic field integral into an inner product between the propagation kernel interpolation weights, which depend only on the observation geometry, and the nodal values ​​of the reflection coefficients, which depend only on seabed parameters. This method requires only the nodal function values ​​of the reflection coefficients and eliminates the need for derivatives of any order, thus avoiding the calculation of high-order derivatives and their corresponding propagation and sum matrices, reducing computational complexity, and further improving computational efficiency. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of the calculation method of the present invention; Figure 2 The graph shows the observation data, where black dots represent the observation data, green lines represent the relationship between the final accurately inverted reflection coefficient and the incident angle, and red lines represent the relationship between the final decoupled inversion reflection coefficient and the incident angle. Figure 3 The graph shows the node values, where the solid black line represents the relationship between the node value and the node number, and the dashed black line represents the node value when the node number is 450. Figure 4 Here is a diagram of the propagation matrix, where A) represents the magnitude of the propagation matrix, B) represents the real part of the propagation matrix, and C) represents the imaginary part of the propagation matrix; Figure 5 Response vector diagram, where, The magnitude of the response vector Represents the real part of the response vector, Represents the imaginary part of the response vector; Figure 6 The diagram shows the decoupled spherical wave reflection coefficient and calculation error, where A represents the decoupled reflection coefficient and the accurate reflection coefficient, and B represents the calculation error of the decoupled reflection coefficient. Figure 7 The graph shows the relationship between the P-wave velocity of the first layer and the number of iterations, with the green line representing the exact inversion and the red line representing the decoupled inversion. Figure 8 The graph shows the relationship between the first-layer shear wave velocity and the number of iterations, with the green line representing the exact inversion and the red line representing the decoupled inversion. Figure 9 The graph shows the relationship between the density of the first layer and the number of iterations. The green line represents the exact inversion, and the red line represents the decoupled inversion. Figure 10The graph shows the relationship between the thickness of the first layer and the number of iterations. The green line represents the exact inversion, and the red line represents the decoupled inversion. Figure 11 The graph shows the relationship between the longitudinal wave velocity in half-space and the number of iterations. The green line represents the exact inversion, and the red line represents the decoupled inversion. Figure 12 The graph shows the relationship between the transverse wave velocity in half-space and the number of iterations. The green line represents the exact inversion, and the red line represents the decoupled inversion. Figure 13 The graph shows the relationship between the density of the half-space and the number of iterations. The green line represents the exact inversion, and the red line represents the decoupled inversion. Figure 14 To show the relationship between the normalized objective function and the number of iterations, the green line represents the exact inversion and the red line represents the decoupled inversion. Detailed Implementation

[0013] The technical solution of the present invention will now be described in detail with reference to the accompanying drawings. It is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the protection scope of the present invention.

[0014] This invention relates to a method for calculating and inverting the sound field of spherical wave reflection based on global interpolation decoupling, such as... Figure 1 As shown, its main principle is as follows.

[0015] According to the Sommerfeld integral expression, the reflection coefficient of spherical waves on the seabed can be written as: (1) In equation (1), It is the response kernel of formation parameters. This is the wave field propagation kernel that determines the properties of wave field propagation. Among them, The imaginary unit, Provide information on the physical properties and layering of the seafloor strata. Indicates the horizontal wavenumber. For vertical wavenumber, The horizontal distance from the point source to the receiving point. This represents the sum of the distances from the point source and the receiving point to the seabed. ω is the angular frequency. Indicates the propagation distance of the wave field. The wave number in water. Indicates the horizontal wavenumber. It is an exponential function. It is a zero-order Bessel function of the first kind.

[0016] The accurate spherical wave reflection coefficient can be obtained from equation (1). However, the calculation process requires integration of the rapidly oscillating function, making the computational efficiency of the spherical wave reflected sound field significantly lower than that of the plane wave method. In the inversion process, each evaluation of the objective function requires expensive integration of the oscillating function, making the inefficient spherical wave forward modeling a major bottleneck restricting the inversion. To address the above problems, one acceleration approach is to decouple the response kernel of the formation parameters from the wave field propagation kernel, allowing the latter to be calculated and stored offline, while the online iteration in the inversion process only needs to update the former cheaply. This decoupling framework effectively improves the inversion efficiency. However, at present, this work achieves decoupling through piecewise local Taylor expansion, and its response end consists of the derivatives of the plane wave reflection coefficient at each sub-interval configuration point. Therefore, the offline stage must pre-calculate and store multiple sets of propagation kernel matrix information corresponding to the derivatives as the expansion order increases, and each online iteration depends on the repeated calculation of the higher-order derivatives of the reflection coefficient.

[0017] This invention addresses the problems of low computational efficiency in existing spherical wave reflection inversion methods and reliance on higher-order derivative information in existing local Taylor expansion decoupling methods. It proposes a method for calculating and inverting the acoustic field of spherical wave reflection based on global interpolation decoupling. The algorithm flow is as follows: Figure 1 As shown, the specific steps include: Step 1: Input observation data Initialize seabed elastic parameters .

[0018] This embodiment uses single-reflection point data under layered seabed conditions as an example to describe the technical solution. The observation data is as follows: Figure 2 As shown, the data is obtained by forward modeling of a solid seabed medium model consisting of a layer of seabed sediment over an infinite half-space using equation (1) and by adding 5% random noise. Figure 2 (Black dots) The minimum incident angle of the data is 1°, and the maximum incident angle is 79°, distributed at 2° intervals. The data frequency shown is 100Hz, and the distance from the source and receiver to the reflecting interface is 50m (i.e.,...). ), and It can be calculated based on geometric relationships. The ocean current P and S wave velocities and seafloor P and S wave attenuation information are used as known values ​​in both forward and inverse models. The seafloor stratigraphic properties and layering information are randomly initialized as follows: The components represent, in order, the first layer longitudinal wave velocity (m / s), the first layer transverse wave velocity (m / s), and the first layer density (kg / m³). 3 First layer thickness (m), half-space longitudinal wave velocity (m / s), half-space transverse wave velocity (m / s), half-space density (kg / m³) 3 ).

[0019] Step 2: Set fixed response kernel interpolation nodes within the integration interval and select the interpolation method; the nodes can be equally spaced or non-equally spaced.

[0020] As a specific implementation method, based on the input observation data information in the first step, the integration upper limit is... With the allowable truncation error The relationship between them can be determined by the following formula: ; In the formula, For wave number, for The maximum value of the response kernel within the interval. The total propagation distance of the wave field from the excitation point, through reflection from the seabed, to the receiving point. Representing observation data The sum of the distances from the corresponding excitation point and receiver point to the seabed reflector surface, This indicates the allowable truncation error. This embodiment sets an upper limit for integration. .

[0021] Then, fixed response kernel interpolation nodes are set within the integration interval. Among them, the node spacing The spacing can be either equal or non-equal. This embodiment... and Within the interval, nodes of 450 and 50 are set at equal intervals according to the response kernel angle, respectively. Nodes are numbered from 1 to 500. The relationship between node values ​​and node numbers is as follows: Figure 3 As shown, this relationship is affected by seabed parameters, which are initialized in the first step as an example. The interpolation method chosen in this embodiment is cubic spline interpolation.

[0022] Step 3: Calculate the propagation matrix using interpolation weights and wavefield propagation kernel; the propagation matrix depends only on the observation geometry (frequency, transmission distance, transmission point location, etc.) and is independent of seabed parameters, and is pre-calculated and stored once in the offline stage.

[0023] In one specific implementation, in this step, the interpolation weights are selected according to the interpolation method set in the second step. Indicates the domain is The The interpolation weights. The weights of the cubic spline interpolation used in this embodiment can be written as: In the formula, , , , For the Dirac function, The horizontal wave number; ; and Represent matrices respectively The , elements, matrix ,in , .

[0024] The propagation matrix can be constructed as follows ; Among the elements Calculated using the following formula: ; in The wave field propagation kernel determines the properties of wave field propagation; The imaginary unit, The index of the interpolation node. , The number of interpolation nodes set for the second step. The component index in the observed data. (The upper limit of the value is determined by the number of observed data components.) The first in the observation data The horizontal distance from the point source to the receiving point corresponding to each component. Indicates the first in the observation data The sum of the distances from the point source and the receiving point to the seabed corresponding to each component. The first in the observation data The angular frequency corresponding to each component The first in the observation data The vertical wavenumber corresponding to each component The wave number in water. Indicates the horizontal wave number. It is an exponential function. It is a zero-order Bessel function of the first kind.

[0025] Figure 4 The propagation matrix obtained by the above calculation in this embodiment is given, and the parameters used are determined by the relevant parameters in the first step.

[0026] Step 4: Calculate the response vector from the response kernel; the response vector contains the function value information of the response kernel at the above fixed interpolation nodes, which depends on the seabed parameters.

[0027] In one specific implementation, this step involves calculating the response kernel function value at the fixed interpolation node set in the second step. The response vector is obtained. The corresponding vector calculated from the parameters given in the first step of this embodiment is as follows: Figure 5 As shown.

[0028] Step 5: Calculate the decoupled spherical wave reflection sound field by taking the inner product of the propagation matrix obtained in Step 3 and the response vector obtained in Step 4. The calculation expression is as follows: ;in, This matrix contains information about the spatial propagation of the wave field and the distribution of discrete nodes within the integration interval, but does not contain information about the seafloor strata properties and layering that determine the plane wave reflection coefficient. In this embodiment, the observation parameters from the first step and the initial seabed parameters are used. The obtained decoupled spherical wave reflection coefficient and its calculation error are as follows: Figure 6 As shown.

[0029] Step 6: Construct the objective function using the spherical wave reflection sound field and update the seabed parameters based on the objective function.

[0030] This embodiment uses a least squares objective function. in . , , The value corresponds to the observed data. The corresponding values ​​of each component.

[0031] The method for updating seabed parameters is as follows: ; in, Indicates the first iteration real number ; For the first The step size of the next iteration can be obtained by the following formula: ,in It is the damping coefficient. For Jacobian matrices, Depend on Constructing diagonal elements. With greater damping, the update more closely approximates a stable gradient descent; with decreasing damping, the update gradually approaches the standard Gaussian-Newton direction. In this embodiment, the parameters to be inverted... With the number of iterations The relationship of change is as follows Figures 7-13 As shown in the figure. The relationship between the normalized objective function and the number of iterations during the iteration process is as follows. Figure 14 As shown.

[0032] from Figures 7-13 and Figure 14The given inversion process is basically consistent overall, except for the errors introduced by the decoupled forward modeling calculation (such as...). Figure 6 The number of nodes will have some impact on the convergence process, but it will not significantly affect the final inversion result. Of course, further increasing the number of nodes can further reduce the error introduced by the forward modeling.

[0033] Step 7: Condition judgment: Return to step (4) or output the current seabed parameters.

[0034] In this embodiment, the number of iterations is set to 10, that is, when Stop iterative calculation and output when the time comes The inversion results are shown in Table 1. Table 1 Inversion Results .

[0035] As shown in Table 1, the inversion results are essentially the same as the exact inversion under the parameters given in the examples, and the minor computational errors are negligible. This indicates that decoupling inversion can replace exact inversion in terms of computational accuracy. The above example was calculated on an Intel(R) Core(TM) i9-9980HK CPU @ 2.4GH platform. The exact inversion took 28.717s, while the decoupling inversion took 5.8652s in the offline propagation matrix calculation stage and 2.5953s in the online inversion stage, resulting in a total speedup of 3.3942 times and an online inversion speedup of 11.065 times. As the number of computation iterations increases, the total speedup will approach the online inversion speed. It should be noted that the above speedup will vary depending on the specific computing environment, but this does not affect the overall conclusion that decoupling inversion improves computational efficiency.

Claims

1. A method for calculating and inverting the sound field of spherical wave reflection based on global interpolation decoupling, characterized in that, The method Includes the following steps: Step 1: Input observation data Initialize seafloor strata properties and layering information ; Step 2: Set fixed interpolation nodes within the integration interval and select the interpolation method; the nodes can be equally spaced or non-equally spaced. Step 3: Calculate the propagation matrix using interpolation weights and wavefield propagation kernel; the propagation matrix depends only on the observation geometry, including frequency, transmission distance, and transmission / reception point location, and is independent of seabed parameters. It is pre-calculated and stored once in the offline stage. The interpolation weights are selected according to the interpolation method set in the second step. Indicates the domain is The Interpolation weights; the propagation matrix is ​​constructed as follows: , Each element Calculated using the following formula: ； in The wave field propagation kernel determines the properties of wave field propagation; The imaginary unit, The index of the interpolation node. , The number of interpolation nodes set for the second step; The component index in the observed data. The upper limit of the value is determined by the number of observed data components. The first in the observation data The horizontal distance from the point source to the receiving point corresponding to each component. Indicates the first in the observation data The sum of the distances from the point source and the receiving point to the seabed corresponding to each component. The first in the observation data The angular frequency corresponding to each component The first in the observation data The vertical wavenumber corresponding to each component The wave number in water. Indicates the horizontal wavenumber. It is an exponential function. It is a zero-order Bessel function of the first kind. For the first The total propagation distance from the excitation point corresponding to each component to the receiving point via reflection from the seabed; Step 4: Calculate the response vector from the response kernel; the response vector contains the function value information of the response kernel at the above fixed interpolation nodes, depending on the seabed parameters; Calculate the response kernel function value at the fixed interpolation nodes set in the second step. ,in In response to the kernel interpolation node, , The product is the upper limit of the integral, which yields the corresponding vector. , This is a vector containing information on the physical properties and stratification of the seafloor strata. Step 5: Calculate the decoupled spherical wave reflection sound field by taking the inner product of the propagation matrix obtained in Step 3 and the response vector obtained in Step 4. The calculation expression is as follows: ； matrix It includes information on the spatial propagation of the wave field and the distribution of discrete nodes within the integration interval, but does not include information on seafloor strata properties and layering that determine the plane wave reflection coefficient. , For the transpose operator, Indicates the matrix Transpose; vector It contains information on seafloor strata properties and layering that determine the plane wave reflection coefficient. ; Step 6: Construct the objective function using the spherical wave reflection sound field and update the seabed parameters based on the objective function; Step 7: Conditional judgment: Return to step 4 or output the current seabed parameters.

2. The method for calculating and inverting the sound field of spherical wave reflection based on global interpolation decoupling according to claim 1, characterized in that, In the second step, based on the observation data input in the first step... Set fixed response kernel interpolation nodes within the integration interval. Among them, the node spacing Upper limit of integration, whether the spacing is equal or non-equal. Determine using the following formula: ； In the formula, For wave number, for The maximum value of the response kernel within the interval. The total propagation distance of the wave field from the excitation point, through reflection from the seabed, to the receiving point. Representing observation data The sum of the distances from the corresponding excitation point and receiver point to the seabed reflector surface, Indicates the allowable truncation error.

3. A system for calculating and inverting the sound field of spherical wave reflection based on global interpolation decoupling, characterized in that, The system operates the spherical wave reflection sound field calculation and inversion method based on global interpolation decoupling as described in any one of claims 1-2.

4. A device for calculating and inverting the sound field of spherical wave reflection based on global interpolation decoupling, characterized in that, The device is equipped with the spherical wave reflection sound field calculation and inversion system based on global interpolation decoupling as described in claim 3.

Images