A fast calculation method for broadband sound field at complex undulating boundaries in shallow seas

Abstract

This application belongs to the field of marine acoustics technology and provides a rapid calculation method for broadband sound fields of complex undulating boundaries in shallow seas. The method includes: establishing an equivalent source model; extracting several discrete frequency points from the target frequency band to generate a discrete frequency point combination; sequentially extracting discrete frequency points from the discrete frequency point combination as first discrete frequency points and determining the source intensity of each virtual sound source in the equivalent source model at each first discrete frequency point by solving the dimensional matrix equation, until the trend of the source intensity of each virtual sound source changing with the first discrete frequency point satisfies the frequency domain extension condition of the source intensity; determining the source intensity of each virtual sound source at each second discrete frequency point based on the frequency domain extension of the source intensity; and calculating the broadband sound field of the target shallow sea area based on the source intensity of each virtual sound source at each discrete frequency point in the discrete frequency point combination. The rapid calculation method provided by this application can significantly reduce the computational workload of calculating broadband sound fields in shallow seas using equivalent source models under complex undulating boundary conditions.

Description

Technical Field

[0001] This application belongs to the field of marine acoustics technology, specifically, it provides a fast calculation method for broadband sound fields of complex undulating boundaries in shallow seas. Background Technology

[0002] Sound field computation, as a core foundation of underwater acoustic propagation research, plays a crucial role in sonar system design, underwater target detection, and marine communication assessment. With the increasing demand for broadband sound field computation in underwater acoustic applications, modern underwater acoustic technologies, represented by broadband matched field processing, increasingly rely on high-density sound field computation across multiple frequencies and spatial locations.

[0003] Among existing numerical methods for calculating marine acoustic fields, the Equivalent Source Method (ESM) has been widely used due to its unique physical mechanism. This method configures a virtual equivalent source array at undulating boundaries (sea surface, seabed), uses the Green's function superposition principle to accurately solve the Helmholtz equation, and characterizes the boundary reflection and transmission effects based on the medium continuity condition. It effectively avoids the angle-dependent error of traditional reflection coefficient models and exhibits high accuracy in complex marine environments such as range-dependent environments.

[0004] Currently, mainstream numerical methods for calculating ocean acoustic fields based on equivalent source models all employ a frequency sweeping calculation strategy. This involves calculating the acoustic field of a virtual equivalent source array at each discrete frequency point within a set bandwidth, and then synthesizing it to obtain a broadband response. As the complexity of boundary conditions increases and the calculation frequency rises, the spatial density of the virtual equivalent source array generally needs to be increased to ensure calculation accuracy. In this case, the acoustic field calculation time will increase linearly with the bandwidth of the calculated frequency band. Since broadband acoustic field calculations typically account for 80% to 90% of the total computation in applications such as underwater target detection, the computational efficiency defects that existing numerical methods cannot overcome when facing complex undulating boundaries and broadband scenarios have severely restricted the real-time processing capability and practical application scope of broadband matched field processing methods. Summary of the Invention

[0005] This application provides a method for rapidly calculating the broadband sound field of a complex, undulating boundary in shallow seas, through embodiments, comprising the following steps:

[0006] S1. Establish an equivalent source model for the actual sound source. The equivalent source model includes three equivalent source groups for the sea surface reflection field, seabed transmission field, and seabed reflection field, which are used to represent the actual sound source, respectively. The number of virtual sound sources included in the equivalent source group is determined based on the upper frequency limit of the target frequency band and the severity of boundary undulations in the target shallow sea area.

[0007] S2, extract several discrete frequency points from the target frequency band to generate a combination of discrete frequency points;

[0008] S3, sequentially extract discrete frequency points from the discrete frequency point combination as the first discrete frequency point, and for each extracted first discrete frequency point, solve... The dimensional matrix equation determines the source intensity of each virtual sound source in the equivalent source model at the first discrete frequency point, until the trend of the source intensity of each virtual sound source changing with the first discrete frequency point satisfies the frequency domain extension condition of the source intensity, where, The total number of virtual sound sources;

[0009] S4. For each virtual sound source, perform frequency domain extension of the source strength based on the source strength at each first discrete frequency point to determine the source strength of the virtual sound source at each second discrete frequency point, wherein the second discrete frequency point is the discrete frequency point other than the first discrete frequency point in the combination of discrete frequency points.

[0010] S5 calculates the broadband sound field of the target shallow sea area based on the source intensity of each virtual sound source at each discrete frequency point in the discrete frequency point combination.

[0011] The fast calculation method for broadband sound fields in shallow seas with complex undulating boundaries provided in this application is based on an in-depth analysis of the frequency domain structural characteristics of each virtual sound source intensity in the equivalent source model for calculating broadband sound fields in shallow seas. It utilizes the modal information contained in the frequency domain structure of the virtual sound source intensity caused by the waveguide environment under complex boundary conditions in shallow seas, and then adaptively improves the modal decomposition technique for time-domain evolved signals. By judging the frequency domain extension conditions, it adaptively switches from the computationally intensive frequency-point matrix equation solving operation to the rapid extension of the source intensity frequency domain. While ensuring the accuracy of the broadband sound field calculation results, it greatly reduces the amount of computation and improves the speed of sound field calculation. Attached Figure Description

[0012] Figure 1 A basic principle diagram of a shallow sea undulating boundary environment and an equivalent source model;

[0013] Figure 2 This is a flowchart of a conventional method for calculating broadband acoustic fields in shallow seas using an equivalent source model.

[0014] Figure 3 For use Figure 2 A schematic diagram of the shallow sea broadband sound field obtained by the process calculation;

[0015] Figure 4 A flowchart illustrating a rapid calculation method for a broadband sound field at a complex undulating boundary in shallow water, according to an embodiment of this application;

[0016] Figure 5 for Figure 3 Chinese virtual sound source , and The diagram shows the curve of the real part of the source strength as a function of frequency, where (a) is... The schematic diagram of the curve of the real part of the source strength as a function of frequency is shown in (b). The schematic diagram of the curve of the real part of the source strength as a function of frequency is shown in (c). A schematic diagram of the curve showing the real part of the source strength as a function of frequency;

[0017] Figure 6 for Figure 3 Chinese virtual sound source , and The diagram shows the curve of the real part of the source strength as a function of frequency, where (a) is... The schematic diagram of the curve of the real part of the source strength as a function of frequency is shown in (b). The schematic diagram of the curve of the real part of the source strength as a function of frequency is shown in (c). A schematic diagram of the curve showing the real part of the source strength as a function of frequency;

[0018] Figure 7 for Figure 3 Chinese virtual sound source , and The diagram shows the curve of the real part of the source strength as a function of frequency, where (a) is... The schematic diagram of the curve of the real part of the source strength as a function of frequency is shown in (b). The schematic diagram of the curve of the real part of the source strength as a function of frequency is shown in (c). A schematic diagram of the curve showing the real part of the source strength as a function of frequency;

[0019] Figure 8 This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 100Hz, calculated using the method provided in this application in Specific Embodiment 1.

[0020] Figure 9 This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 100Hz, calculated using conventional methods.

[0021] Figure 10 This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 500Hz, calculated using the method provided in this application in Specific Embodiment 1.

[0022] Figure 11 This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 500Hz, calculated using conventional methods.

[0023] Figure 12This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 600Hz, calculated using the method provided in this application in Specific Embodiment 1.

[0024] Figure 13 This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 600Hz, calculated using conventional methods.

[0025] Figure 14 This is a schematic diagram showing the sound pressure level calculation results of a fixed receiving point in the frequency range of 50 Hz to 600 Hz in specific embodiment 1. Detailed Implementation

[0026] The present application will now be further described based on preferred embodiments and with reference to the accompanying drawings.

[0027] In the description of the embodiments of this application, it should be noted that the terms "upper," "lower," "inner," and "outer," etc., indicating orientation or positional relationships based on the orientation or positional relationships shown in the accompanying drawings, or the orientation or positional relationships commonly used when the product of this application is in use, are only for the convenience of describing this application and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation on this application. Furthermore, in the description of this application, the terms "first," "second," etc., are used to distinguish different units, but these are not limited by the manufacturing order, nor should they be construed as indicating or implying relative importance. Their names may differ in the detailed description and claims of this application. In addition, for ease of understanding, various components in the drawings have been enlarged or reduced, but this is not intended to limit the scope of protection of this application.

[0028] The vocabulary used in this specification is for illustrative purposes and is not intended to limit the scope of this application. It should also be noted that, unless otherwise expressly specified and limited, the terms "set," "connected," and "linked" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection, a direct connection, or an indirect connection via an intermediate medium; or they can refer to the internal communication between two components. Those skilled in the art will understand the specific meaning of these terms in this application.

[0029] To clearly illustrate the key improvements of the technical solution in this application, we first provide a detailed explanation of the current conventional method for calculating shallow sea broadband sound fields based on equivalent source models and its computational bottlenecks.

[0030] Figure 1This diagram illustrates the basic principle of a shallow sea undulating boundary environment and an equivalent source model. This equivalent source model is used to calculate the broadband sound field of the broadband underwater acoustic signal emitted by the sound source in seawater and the seabed medium. The calculation idea is to simulate the reflection and transmission fields on the undulating boundary by using the sound fields generated by several sets of equivalent sources, and to solve the source strength of each equivalent source based on the boundary continuity condition and the Green's function superposition principle. Finally, by superimposing the sound fields generated by these equivalent sources, the accurate broadband sound field distribution generated by the real sound source in seawater and the seabed medium is obtained.

[0031] by Figure 1 Taking the Pekeris waveguide in shallow water as an example, a broadband signal emitted from a real sound source (shown as a red dot in the figure) will be reflected by the sea surface and then reflected and transmitted by the undulating seabed during its propagation in the shallow water waveguide environment, thus forming a sound image in the seawater (where the speed of sound and density of sound are respectively...). and ) and the seabed medium (the sound velocity, density, and attenuation parameters of the seabed medium are respectively ) and the seabed medium. , and When calculating the broadband sound field at any point in the seawater or seabed using the equivalent source model, the real sound source can be replaced by three sets of discrete equivalent source combinations. The sound fields generated by these three sets of equivalent source combinations can then be used to replace the sea surface reflection field, seabed reflection field, and seabed transmission field generated by the real sound source, respectively, thereby finally determining the broadband sound field at any point in the seawater or seabed.

[0032] Specifically, the interface between seawater and air is considered to have a depth of... On a flat sea surface, several virtual sound sources are discretely generated. These virtual sound sources are then translated upwards along the normal direction from the sea surface to obtain the equivalent source set of the sea surface used to represent sea surface reflection. This refers to the combination of several equally spaced black dots located above the sea surface in the diagram; the interface between seawater and the seabed is considered to be dynamic. On the undulating seabed, two sets of discrete virtual sound sources are generated. One set of virtual sound sources is translated upwards along the normal direction from the seabed to obtain the first equivalent seabed sound source set used to represent seabed transmission. (That is, the combination of several equally spaced black dots above the undulating seabed in the figure), another set of virtual sound sources is translated downwards along the normal direction from the seabed to obtain the second equivalent seabed sound source group used to represent seabed reflection. (That is, the combination of several equally spaced black dots located below the undulating seabed in the figure).

[0033] Equivalent sound source group , and The spacing between the discrete virtual sound sources in the image. The translational distance of each virtual sound source relative to the sea surface, and and The distance that each virtual sound source moves relative to the seabed is generally set to... to , The wavelength is the frequency to be calculated; , and The number of virtual sound sources is determined based on the size of the target shallow sea area for sound field calculation, to ensure that the sound field generated by each equivalent sound source group can accurately reflect the propagation characteristics of the underwater acoustic signal emitted by the real sound source in the target shallow sea area.

[0034] Finish Figure 1 Once the equivalent source model shown is constructed, the sound field at any receiving point in the seawater and seabed medium can be determined through the following derivation process.

[0035] First, in an ocean waveguide, a point source should satisfy the Helmholtz equation shown in equation (1):

[0036] (1),

[0037] (1) In the formula, , These represent the position vectors of the receiving point and the point sound source, respectively. For wave number, The angular frequency of the underwater acoustic signal emitted by a point sound source. The sound pressure Green's function, This is the Dirac function. Figure 1 In the two-dimensional scene shown, the sound pressure Green's function for:

[0038] (2),

[0039] in, It is a zeroth-order Hankel function of the first kind. It is the imaginary unit.

[0040] and Corresponding displacement Green's function for:

[0041] (3),

[0042] in, For the density of the medium, It is the normal unit vector.

[0043] Sound field in seawater and seabed media in shallow sea areas We need to solve the Helmholtz equations satisfied by the sound field shown in equation (4) to obtain:

[0044] (4).

[0045] During the solution process, the sea surface and seabed satisfy the pressure release boundary condition and the continuity boundary condition, respectively, namely:

[0046] (5),

[0047] in, , These are receiving points located on the sea surface and the seabed, respectively. , These are the sound pressure function and the normal displacement function, respectively. , They represent Sound pressure on the seawater side and sound pressure on the seabed side. , They represent The normal displacement on the seawater side and the normal displacement on the seabed medium side.

[0048] When using Figure 1 shown , and When performing sound field calculations using the equivalent source model, the solution to equation (1) can be represented by superimposing basis functions (i.e., Green's functions). Specifically, let... , , Each of them contains , , A virtual sound source, then , , It can be represented by a vector formed by the source intensities of each virtual sound source contained therein. , , The locations of each virtual sound source can also be represented as position vectors. , , ,Right now:

[0049] (6),

[0050] (6) In the formula and They represent The Middle The source strength and location of a virtual sound source. and They represent The Middle The source strength and location of a virtual sound source. and They represent The Middle The source strength and location of each virtual sound source, indicated by the superscript. This is a transpose operation. It should be noted that, due to... and The virtual sound sources in the image are obtained by symmetrically shifting the image from the seabed to both sides, therefore they always exist. .

[0051] Furthermore, the sea surface and seabed can be discretized into... , Given a set of nodes, the following vector representation of the sea surface nodes is obtained. Vector representation of seabed nodes :

[0052] (7),

[0053] (7) In the formula For the first The location of each sea surface node For the first The location of each seabed node is determined by the number of equations in the system to ensure the uniqueness of the solution when solving the matrix equations. Generally, let... , That is, each sea surface node corresponds to A virtual sound source in one, each seabed node corresponds to and A pair of mirrored virtual sound sources.

[0054] for Figure 1 The shallow sea waveguide environment shown can be considered as consisting of a direct wave sound field and a superposition of sound fields reflected from the sea surface and seabed boundary; the seabed medium is a homogeneous single-layer medium, with only a seabed transmitted sound field existing, and the frequency is... Sound field in seawater produced by real sound sources and sound field in seabed medium It can be represented as:

[0055] (8).

[0056] in, , These are the wave numbers of seawater and the seabed medium, respectively. For the real sound source at the receiving point The direct-wave sound field at a given location can be determined by the source intensity of the real sound source and... , The distance between them is directly determined.

[0057] Obviously, the values ​​of each term in equation (8) are all frequency-dependent. Only by first calculating the frequency points of each virtual sound source can we obtain the values ​​of each virtual sound source. Only with the corresponding source strength can each frequency point be calculated. The sound field of the seawater and the sound field of the seabed.

[0058] To solve for the source intensity of each virtual sound source, the boundary conditions shown in equation (5) need to be applied at each sea surface node and seabed node, thus obtaining equation (9):

[0059] (9),

[0060] in, For the real sound source at the sea surface node The direct wave sound field at that location, , These represent the direct wave sound pressure and direct wave displacement at the actual sound source at the seabed node, respectively.

[0061] For each sea surface node and seabed nodes All of these can yield a set of equations as shown in equation (9), which can be used to... A system of equations in the form of The matrix equation form, as mentioned above, is usually represented as follows: Set to with Consistency ensures that the matrix equation can be solved uniquely; therefore, the boundary conditions can be used to obtain solutions for each frequency point. The matrix equation is expressed as shown in equation (10):

[0062] (10)

[0063] in:

[0064] , Equivalent sound source groups , The Green's function matrix of the sound pressure of each virtual sound source at each sea surface node is expressed as follows:

[0065] (11);

[0066] , , Equivalent sound source groups , , The Green's function matrix of the sound pressure of each virtual sound source at each seabed node is expressed as follows:

[0067] (12);

[0068] , , Equivalent sound source groups , , The displacement Green's function matrix of each virtual sound source at each seabed node is expressed as follows:

[0069] (13).

[0070] From equations (11) to (13), it can be seen that the number of rows in the coefficient matrix on the left side of the matrix equation (10) is... The number of columns is ,Depend on , It can be seen that the dimension of this matrix equation is dimension.

[0071] Therefore, determining the source strength of each virtual sound source is transformed into solving the problem shown in (10). The problem of 3D matrix equations. Equation (10) can be solved using standard numerical methods, such as Gaussian elimination or LU decomposition, thus obtaining... , , The source strength of each virtual sound source is obtained by substituting the source strength of the virtual sound source into equation (8). Corresponding sound field in seawater and sound field in seabed medium .

[0072] Figure 2 This demonstrates how, when the actual sound source is a broadband sound source, the equivalent source model described above can be directly used to calculate the signal at any receiving point in the target shallow sea area. A flowchart of a broadband sound field. For example... Figure 2 As shown, when the underwater acoustic signal emitted by the real sound source is a broadband signal, the target frequency band for which sound field calculation is required is obtained after discretization. discrete frequency points At each frequency point All of these require consideration of equation (10). Solve the dimensional matrix equation, and then substitute the virtual source strength into equation (8) to perform the summation operation of the sound field.

[0073] Figure 3 In one specific embodiment, the use of Figure 2 The diagram shows a broadband sound field of the target shallow sea area calculated using the shown process. Figure 3 The sound field shown comprises three dimensions. Two dimensions are depth and distance, used to characterize the positions of various receiving points in the target shallow sea area where sound field calculations are to be performed. The third dimension is frequency, i.e., for each receiving point, through... Figure 2 The process calculates each discrete frequency point. sound field or .

[0074] According to the knowledge of matrix operations, when using numerical methods such as Gaussian elimination and LU decomposition to solve equation (10), matrix decomposition is required first, and its computational complexity is O(n). Then, the trigonometric equations are solved, with a computational complexity of O(n log n). ,exist When the number of elements is large, the overall complexity of the solution can be expressed as: Use equation (8) to perform The calculation of the sea surface reflection field and seabed reflection field of the virtual sound source involves Multiplication and The second addition operation is performed. The calculation of the seabed transmission field of the virtual sound source involves Multiplication and The time complexity of the above calculations can be expressed as: (The original text contains several typographical errors and inconsistencies, which have been omitted from the translation.) .

[0075] The above analysis reveals that the speed bottleneck in calculating broadband sound fields in shallow sea environments with complex undulating boundary conditions using the equivalent source model mainly stems from the limitations of... , , Calculation of the intensity of each virtual sound source: Due to the significant undulations at the seabed boundary in shallow waters, it is generally necessary to reduce the spacing between the virtual source intensities to ensure accurate calculation results. , , The total number of virtual source strengths included Significantly increased, even reaching The order of magnitude, that is, in the matrix equation (10) established for each frequency point, the total number of unknowns reaches 3000, and the dimension of the coefficient matrix reaches Dimension, its computational scale reaches Therefore, when using the above method to calculate the broadband sound field, as the number of virtual source strengths increases, the computational workload of solving the matrix equations for all frequency points in the target broadband will expand rapidly.

[0076] To address the issue of drastically increased computational load when using virtual source models for broadband sound field calculations under complex boundary conditions in shallow seas, this application improves upon existing calculation processes that involve solving matrix equations for each frequency point, proposing a rapid calculation method for broadband sound fields in complex, undulating shallow sea boundaries. Figure 4 As shown, the method includes the following steps:

[0077] Step S1: Establish an equivalent source model of the actual sound source. The equivalent source model includes three equivalent source groups: the sea surface reflection field, the seabed transmission field, and the seabed reflection field, which are used to represent the actual sound source. The number of virtual sound sources in the equivalent source group is determined based on the upper limit of the target frequency band and the severity of the boundary undulation of the target shallow sea area.

[0078] Step S2: Extract several discrete frequency points from the target frequency band to generate a discrete frequency point combination;

[0079] Step S3: Extract discrete frequency points sequentially from the discrete frequency point combination as the first discrete frequency point. For each extracted first discrete frequency point, solve... The dimensional matrix equation determines the source intensity of each virtual sound source in the equivalent source model at the first discrete frequency point, until the trend of the source intensity of each virtual sound source changing with the first discrete frequency point satisfies the frequency domain extension condition of the source intensity, where, The total number of virtual sound sources;

[0080] Step S4: For each virtual sound source, perform frequency domain extension of the source strength based on the source strength at each first discrete frequency point to determine the source strength of the virtual sound source at each second discrete frequency point, wherein the second discrete frequency point is the discrete frequency point other than the first discrete frequency point in the combination of discrete frequency points.

[0081] Step S5: Calculate the broadband sound field of the target shallow sea area based on the source intensity of each virtual sound source at each discrete frequency point in the discrete frequency point combination.

[0082] The steps S1 to S5 are described in detail below with reference to the accompanying drawings.

[0083] <Establish an equivalent source model for a real sound source>

[0084] Step S1 is used to establish three equivalent source groups of the actual sound source, namely the sea surface equivalent source group mentioned above. First seabed equivalent source group and the second seabed equivalent source group These three equivalent source groups each contain several virtual sound sources, which are used to represent the sea surface reflection field, seabed transmission field, and seabed reflection field generated by the actual sound source in the target shallow sea area.

[0085] , , The specific construction process of each virtual sound source has been described in the previous text. Figure 1 A detailed introduction has been provided, and correspondingly, the expressions for the source strength and location of each virtual sound source are shown in equation (6). For example, it should be known that in equation (6) , , The number of virtual sound sources included is the same; however, in other embodiments, the number of virtual sound sources may also be... The number of virtual sound sources included is set to the same as , The number of virtual sound source pairs generated in pairs varies, for example , , The number of virtual sound sources included are respectively , , The total number of virtual sound sources included in the equivalent source model of the target shallow sea area is... It is equal to Furthermore, it is understandable that, due to , Virtual sound sources in the game are always generated in pairs, therefore it is necessary to ensure that .

[0086] In some embodiments, The value can be determined based on the upper frequency limit of the target frequency band (in this application, the target frequency band refers to the lower frequency limit of the calculated broadband sound field). Up to frequency limit The frequency range constituted by this calculation is used because, in the process of calculating a broadband sound field, each discrete frequency point in the target frequency band needs to pass through the equivalent reflection and transmission fields of each virtual sound source. Therefore, to ensure the accuracy of the sound field calculation results, the upper frequency limit can be used. Determine the lower limit of wavelength in seawater. And set the spacing between adjacent virtual sound sources to of to This is to avoid the problem of excessively large spacing and quantity of virtual sound sources, which could lead to insufficient resolution of the calculated broadband sound field.

[0087] Furthermore, in some preferred embodiments, the intensity of boundary undulations in the target shallow sea area can be further adjusted. The range of values ​​can be optimized. For example, when the seabed boundary of the target shallow sea area is relatively undulating, the spacing between virtual sound sources can be further reduced to depict the undulating boundary with higher resolution, thereby improving the accuracy of the sound field calculation results.

[0088] Constructing Discrete Frequency Point Combinations

[0089] Step S2 is used to construct the discrete frequency point combination corresponding to the target frequency band. In some specific embodiments, the lower limit of the frequency for which sound field calculation needs to be performed can be determined based on the effective frequency range of the underwater acoustic signal actively emitted by the real sound source, or based on the application scenario of the sound field calculation results (such as the analysis of the sound field structure of shallow sea waveguides in characteristic frequency bands). and frequency limit Thus, the target frequency band is obtained. ~ Then, based on the required frequency resolution of the broadband sound field in the target frequency band, several discrete frequency points are extracted to generate a discrete frequency point combination. ,in The total number of discrete frequency points. and .

[0090] <Fast Calculation of Broadband Sound Field Based on Frequency Domain Extension of Virtual Sound Source Intensity>

[0091] Compare Figure 4 and Figure 2 It can be seen that the core difference between the rapid calculation method for shallow sea broadband sound field provided in this application and the conventional method of solving matrix equations at each frequency point to obtain the source intensity and then superimposing them to calculate the broadband sound field lies in the fact that, through steps S3 and S4, the method combines discrete frequency points... The algorithm determines whether the frequency domain structure of the sound source intensity obtained by solving the matrix equation has shown obvious modal evolution characteristics. When this condition is met, the computationally complex matrix solution process is terminated, and the algorithm switches to modal decomposition and principal component analysis of the source intensity evolution characteristics with frequency in the solved frequency band. Using the modal correlation parameters obtained from the principal component analysis, the frequency domain structure of the source intensity in the unsolved frequency band is extended, thereby obtaining the source intensity of all discrete frequency points in the target frequency band.

[0092] The following is a mechanistic analysis of the feasibility of steps S3 and S4 in the method provided in this application. Dynamic mode decomposition technology is a data-driven unsupervised dimensionality reduction algorithm. Its core objective is to extract the dominant dynamic modes of system evolution over a past period by analyzing time series data, and to obtain prediction results of system evolution behavior in the future by weighted combination of these modes for time domain extension.

[0093] The applicant's use Figure 2 Analysis of the shallow sea broadband sound field calculated using the illustrated process reveals that the source intensity of each virtual sound source exhibits obvious modal characteristics that dynamically evolve with frequency in the frequency domain. Figure 5 schematically shown The three virtual sound sources in the middle (i.e. Figure 3 In , and The real part of the source strength of ) The curve showing how the frequency changes. Figure 6 schematically shown The three virtual sound sources in the middle (i.e. Figure 3 In , and The real part of the source strength of ) The curve showing how the frequency changes. Figure 7 schematically shown The three virtual sound sources in the middle (i.e. Figure 3 In , and The real part of the source strength of ) Curves showing frequency variation. Furthermore, curves showing the variation of the imaginary part of each virtual source strength with frequency can be plotted to obtain similar frequency domain variation curves.

[0094] pass Figures 5 to 7 As can be seen, the sound field at the receiving point in a shallow sea waveguide environment is the result of the superposition of direct waves, reflected waves, and transmitted waves propagating through multiple paths, exhibiting obvious structural interference characteristics. These structural interference characteristics are introduced into the source strength information of each virtual sound source during the solution of the matrix equations through the constraint conditions of undulating boundaries. Furthermore, since the boundary constraints are frequency-dependent, the frequency domain structure of the solved virtual sound source strength exhibits frequency-dependent modal characteristics. Therefore, time-domain dynamic mode analysis techniques can be introduced into the extension of the frequency domain structure of the virtual sound source strength, that is, extracting several discrete frequency points within a portion of the target frequency band (in this application, this frequency band is referred to as the first frequency band, and the discrete frequency points extracted within the first frequency band are referred to as the first discrete frequency points), and then... Figure 2 Using the same matrix equation solving method as shown, the frequency domain structure of the source strength in the frequency band is obtained; then, dynamic mode decomposition of the source strength frequency domain structure in the first frequency band is performed; then, several discrete frequency points are extracted in the remaining frequency band of the target frequency band (in this application, this frequency band is called the second frequency band, and the discrete frequency points extracted in the second frequency band are called the second discrete frequency points). For each second discrete frequency point, mode decomposition of the source strength frequency domain structure of the previous several frequency points is performed, and the source strength at each second discrete frequency point is obtained by weighted combination of the decomposed modes. That is, the source strength frequency domain structure of the second frequency band is obtained by using the source strength frequency domain structure obtained by solving the matrix equation in the first frequency band and extending it in the frequency domain.

[0095] Obviously, it should be noted that the prerequisite for performing frequency domain structural extension processing based on modal decomposition in the frequency domain is that the source strength of the virtual sound source has modal characteristics in the frequency domain. As can be seen from the previous analysis, the modal characteristics of the source strength of the virtual sound source are related to the environmental structure of the shallow sea target area (such as water depth, seabed undulation, etc.) and the target frequency band. Therefore, in the embodiments of this application, the division of the first frequency band and the second frequency band cannot be determined in advance, but is dynamically determined according to the matrix equation solution results of each discrete frequency point.

[0096] In some optional embodiments, step S3 further includes the following steps:

[0097] Step S31: Obtain the discrete frequency point combination ,in and , , These are the lower and upper limits of the target frequency band, respectively.

[0098] Step S32, extract all unsolved values ​​from the discrete frequency point combination. In the discrete frequency points of the dimensional matrix equation, the first discrete frequency point is extracted in sequence and marked as the first discrete frequency point;

[0099] Step S33: Solve for the first discrete frequency point marked in step S32. The 3D matrix equation is used to obtain the source intensity of each virtual sound source at the first discrete frequency point and extract its real and imaginary parts.

[0100] Step S34: Update the curves of the real and imaginary source strengths of each virtual sound source as a function of the first discrete frequency point.

[0101] Step S35: Determine whether the trend of the real part and imaginary part of the source strength of each virtual sound source changing with the first discrete frequency point satisfies the frequency domain extension condition of the source strength. If not, return to step S32. If so, mark all discrete frequency points in the discrete frequency point combination that are not marked as the first discrete frequency point as the second discrete frequency point, and proceed to step S4.

[0102] Specifically, in step S31, the discrete frequency point combination after discretization of the target frequency band is first obtained. Then, the discrete frequency point that has not yet undergone matrix equation solving is sorted in ascending or descending order from the discrete frequency point combination and marked as the first discrete frequency point. The matrix equation is then solved according to equation (10) (steps S32 and S33). It should be noted that since sound pressure is generally expressed in complex form, the source strength obtained by solving includes both real and imaginary parts.

[0103] As the number of frequency points marked as the first discrete frequency point increases, the curves of the source strength of each virtual sound source changing with frequency will gradually reflect the modal characteristics of the frequency domain. Therefore, after each matrix equation solution operation for a first discrete frequency point, the curves of the real part and imaginary part of the source strength of each virtual sound source changing with frequency are updated in step S34. Then, in step S35, it is determined whether the matrix equation solution operation can be terminated and the frequency domain extension operation of the source strength is performed on the remaining discrete frequency points that have not undergone matrix solution operation.

[0104] After completing step S3, The discrete frequency points in the array will be divided into a first discrete frequency point combination and a second discrete frequency point combination. Without loss of generality, we can assume that the number of the first discrete frequency point combinations is... (obviously) Then the first discrete frequency point combination is: The second discrete frequency point combination is For any virtual sound source, its frequencies at each of the first discrete frequencies are... The source strength at that location has been solved. The dimensional matrix equation yields a sequence that can be represented as a sequence. At each of the second discrete frequency points This needs to be obtained through frequency domain extension in step S4.

[0105] Frequency domain extension conditions can be rationally designed based on the boundary undulations of the virtual sound source intensity frequency domain structure and the shallow sea waveguide environment, as well as parameters such as the target frequency band size and upper and lower frequency limits. Figures 5 to 7 As shown, the real and imaginary parts of the source strength of the virtual sound source exhibit obvious periodic repetition structural characteristics as frequency changes. Therefore, based on the above frequency domain structural characteristics, the curves of the real and imaginary parts of the source strength of each virtual sound source exhibiting periodic repetition, and the number of repetitions reaching a preset value (e.g., more than 10 times), can be used as a frequency domain extension condition. After satisfying this condition, it indicates that the frequency domain structural characteristics of the source strength have shown stable modal evolution characteristics. At this point, the computationally intensive operation of solving the source strength through matrix equations can be stopped, and the source strength of the remaining frequency points can be obtained by extension using the source strength of each frequency point that has been calculated.

[0106] In some preferred embodiments, the frequency domain extension of the source strength for each virtual sound source in step S4 includes the following steps:

[0107] Step S41: Extract the source strength of the virtual sound source at each of the first discrete frequency points. And based on its actual part and the virtual part Construct source strong real part vector The source strong imaginary part vector ,in, The total number of the first discrete frequency points.

[0108] Specifically, as mentioned earlier, the source strength obtained by solving the matrix equation is in complex form, containing both real and imaginary parts. Therefore, the source strength at each of the first discrete frequency points... It can be represented in complex form: ,in, The imaginary unit is used. By extracting the real and imaginary parts of the source strength at each of the first discrete frequency points, the real part vector of the source strength can be obtained. The source strong imaginary part vector .

[0109] Step S42, respectively for and Perform dynamic mode decomposition operation to obtain The corresponding eigenvalue matrix eigenmode matrix sum coefficient vector ,as well as The corresponding eigenvalue matrix eigenmode matrix sum coefficient vector .

[0110] because and All have been The structure exhibits periodic oscillations within the bandwidth of each discrete frequency point. Therefore, through dynamic mode decomposition, several mode functions that determine its shape characteristics, as well as the corresponding mode eigenvalues ​​and weighting coefficients, can be extracted. In subsequent frequency domain extension, combinations of these mode functions can be used to approximate the source strength at each second discrete frequency point.

[0111] In some specific embodiments, for The dynamic mode decomposition operation includes the following steps:

[0112] Step A1, using a length of The windows are for the sequence respectively and Perform a matrix 2D transformation operation to obtain the first real part state matrix. Second real part state matrix .

[0113] Specifically, each with a length of The window in and Slide point by point in the middle, and capture each time. Each element is used as a column vector, and then the column vectors obtained from each sliding cut are concatenated to obtain... and :

[0114] (14)

[0115] (15).

[0116] Preferably, ,in, To perform a floor operation, thus... and The number of rows and columns should be kept as consistent as possible. From the first real part state matrix... To the second real part of the state matrix This reflects the characteristic of the real part of the source strength evolving with the frequency domain within the bandwidth of the first discrete frequency point. Therefore, it can be based on... and The modal decomposition results are further extended to the real part of the source strength at each discrete frequency point outside the bandwidth.

[0117] Step A2, for Perform truncation rank Truncated SVD decomposition:

[0118] (16)

[0119] in, , , They are respectively of A left-singular vector matrix, Victoiresian matrix and A right-singular vector matrix, for The conjugate transpose of the matrix has dimensions of . .

[0120] Step A3: Construct the low-dimensional mapping matrix of the real part based on the following formula. :

[0121] (17)

[0122] in, for The conjugate transpose of the matrix. for The inverse matrix.

[0123] Step A4: Solve the eigenvalue problem shown in the following equation to obtain the eigenvalue matrix. and the low-dimensional eigenvector matrix of the real part :

[0124] (18).

[0125] Specifically, in steps A2 to A4, the state representing the previous state is first used. Perform truncated SVD decomposition, and then combine the decomposition results with the subsequent state. Construct a low-dimensional mapping matrix Finally, through the By performing eigenvalue decomposition, we can obtain a value that reflects... Low-dimensional operators for frequency domain evolution characteristics.

[0126] truncated rank The size of this value determines the number of modes used for subsequent frequency domain extension. Preferably, The value should be greater than or equal to 20 and less than or equal to 30, in order to reduce the amount of computation while ensuring the accuracy of the continuation.

[0127] The results obtained through the above steps , The expressions are as follows:

[0128] (19)

[0129] (20).

[0130] Step A5, determine based on the following formula eigenmode matrix :

[0131] (twenty one).

[0132] Specifically, upon obtaining , Then, it is necessary to use equation (21) to... Increase the dimensionality to obtain an approximation. The frequency domain structure characteristics Each mode, combined with the dimensions of each matrix, indicates that... for A 3D matrix, its specific expression is:

[0133] (twenty one).

[0134] From equation (21), we can see that, Depend on It consists of feature mode vectors, each feature mode including Each element represents the mode in this The values ​​at discrete frequency points.

[0135] Step A6: Determine the coefficient vector by solving the following equation. :

[0136] (twenty two).

[0137] Specifically, in order to perform frequency domain extension of the source strong real part, besides , In addition, it is necessary to determine the contribution of each mode at the frequency points of the frequency domain extension. For this purpose, the calculated values ​​can be used. The real part information of the source strength at the first discrete frequency point is obtained by solving the matrix equation shown in equation (22). dimensional coefficient vector Its specific expression is:

[0138] (twenty three),

[0139] Each element in the vector represents the weight of each characteristic mode at the second discrete frequency point of the frequency domain extension.

[0140] right The steps for performing dynamic mode decomposition can be found in steps A1 to A6. In some specific embodiments, for The dynamic mode decomposition operation includes the following steps:

[0141] Step B1, using a length of The windows are for the sequence respectively and Perform a matrix 2D transformation operation to obtain the first imaginary state matrix. Second imaginary state matrix .

[0142] Step B2, for Perform truncation rank Truncated SVD decomposition:

[0143] (twenty four),

[0144] in, , , They are respectively of A left-singular vector matrix, Victoiresian matrix and A right-singular vector matrix, for The conjugate transpose of the matrix has dimensions of . .

[0145] Step B3: Construct the low-dimensional mapping matrix of the imaginary part based on the following formula. :

[0146] (25),

[0147] in, for The conjugate transpose of the matrix. for The inverse matrix.

[0148] Step B4: Solve the eigenvalue problem shown in the following equation to obtain the eigenvalue matrix. and the low-dimensional eigenvector matrix of the imaginary part :

[0149] (26).

[0150] Step B5, based on the following formula, determines eigenmode matrix :

[0151] (27).

[0152] Step B6: Determine the coefficient vector by solving the following equation. :

[0153] (28).

[0154] Step S43: Determine the virtual sound source at each of the second discrete frequency points based on the following formula. Source strength :

[0155] (29)

[0156] in, for of Power of 1 for of Power of 1 This indicates retrieving the real part of the element within the square brackets. This indicates the expression to retrieve the first element within the curly braces. There are elements, among which... The length of each modal eigenvector in the eigenmode matrix is ​​given. It is the imaginary unit.

[0157] Specifically, step S42 yields the frequency domain extension for performing the source strong real part. , , And frequency domain extension for the imaginary part of the source strength , , Then, it is necessary to determine the frequency for any second discrete frequency point. How to determine the real and imaginary parts? Taking the real part as an example, its... 1 characteristic mode vector It reflects the morphological characteristics of each mode. eigenvalues ​​in complex form This reflects the frequency domain oscillation characteristics. Furthermore, since the length of each characteristic mode is... Taking all the above into consideration, for any second discrete frequency point It can be done (Right now of The power of, ..., of The power () represents the frequency domain oscillation of each mode, and only the last element of each characteristic mode is taken, i.e., the power () Each element, in order to strengthen the real part The extension at that point yields the following equation: Chu Yuanqiang Department The expression:

[0158] (30)

[0159] in, This indicates that the real part of the element within the square brackets is retrieved.

[0160] For the same reasoning, we can obtain the equation shown below. Source of strength and virtuality The expression:

[0161] (31).

[0162] It should be noted that equation (31) is used to obtain... The value of the strong imaginary part of the source is still needed, but it still needs to be obtained through... Characteristic modes in complex form , Take the real part of the result of multiplication.

[0163] Integrating equations (30) and (31), we can obtain the virtual sound source after frequency domain extension. The complete expression for source strength is:

[0164] (32).

[0165] Equation (32) is the complete summation expression that is equivalent to the matrix multiplication form of equation (29).

[0166] <Calculation of Broadband Sound Field>

[0167] After calculating the source strength of the first discrete frequency point and extending the source strength of the second discrete frequency point in steps S3 and S4 respectively, the source strength of the virtual sound source at all discrete frequency points within the target frequency band can be obtained. Then, in step S5, the receiving points in the target sea area where the sound field needs to be calculated are... By using equation (8) to calculate the sound field value at each frequency point, the receiving point can be obtained. Broadband sound field.

[0168] <Specific Implementation Example 1>

[0169] This embodiment uses Figure 1 The shallow sea environment shown is used as the target shallow sea area. The broadband sound field is rapidly calculated using the method provided in this application to verify the effectiveness of the method. The average water depth is 200 m, the sound source depth is 50 m, and the water density is... Speed ​​of sound in water Seabed density speed of sound at the bottom of the sea Seabed attenuation .

[0170] The frequency band of the first discrete frequency point is located in the range of 200 Hz-400 Hz, and the frequency interval between adjacent first discrete frequency points is 1 Hz. The frequency domain extension of the source strength is performed in the frequency band range of 50 Hz to 600 Hz. In addition, in order to compare the method of this application with the conventional process, the sound field calculation based on matrix solution is also performed point by point in the frequency bands of 50 Hz-200 Hz and 400 Hz-600 Hz with a frequency interval of 1 Hz.

[0171] The method provided in this application was written as an executable m-file in Matlab and run on a desktop computer with 256 GB of memory and a CPU clock speed of 2.10 GHz. The source strength was solved point-by-point and its frequency domain extension was implemented through programming. Using conventional methods, the total time for solving the matrix equation point-by-point and calculating the sound field in the 50 Hz to 600 Hz frequency band was approximately 810 seconds. Using the method of this application, the total time for broadband sound field calculation, performing point-by-point calculation in the 200 Hz-400 Hz range and frequency domain extension in other frequency bands, was approximately 240 seconds.

[0172] Figure 8 This is a schematic diagram showing the sound pressure level (SPL) distribution of a two-dimensional sound field at 100Hz calculated using the method provided in this application, for comparison. Figure 9 This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 100Hz, calculated using conventional methods (i.e., solving the matrix equations point by point). Figure 10This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 500Hz, calculated using the method provided in this application. Figure 11 This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 500Hz, calculated using conventional methods. Figure 12 This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 600Hz, calculated using the method provided in this application. Figure 13 This is a schematic diagram of the sound pressure level distribution of a two-dimensional sound field at 600Hz, calculated using conventional methods. It can be seen that in the frequency extension region, the sound field calculation results obtained using the fast calculation method of this application agree well with those obtained using conventional matrix solving methods. The mean square error (MSE) is used to statistically analyze the error of the method of this application compared to the conventional method; the errors at 100Hz, 500Hz, and 600Hz are 0.17 dB, 2.79 dB, and 4.78 dB, respectively.

[0173] Figure 14 The results of sound pressure level calculations at a fixed receiving point (50 m, 500 m) in the frequency range of 50 Hz to 600 Hz are presented. Figure 14 It can be seen that the sound pressure level variation curve at the receiving point obtained by using the method provided in this application is in good agreement with that obtained by using conventional matrix calculation method, and the mean square error (MSE) is 1.54dB.

[0174] The specific embodiments of this application have been described in detail above. For those skilled in the art, several improvements and modifications can be made to this application without departing from the principle of this application, and these improvements and modifications also fall within the protection scope of the claims of this application.

Claims

1. A fast calculation method for broadband sound field at complex undulating boundaries in shallow seas, characterized in that, Includes the following steps: S1. Establish an equivalent source model for the actual sound source. The equivalent source model includes three equivalent source groups for the sea surface reflection field, seabed transmission field, and seabed reflection field, which are used to represent the actual sound source, respectively. The number of virtual sound sources included in the equivalent source group is determined based on the upper frequency limit of the target frequency band and the severity of boundary undulations in the target shallow sea area. S2, extract several discrete frequency points from the target frequency band to generate a combination of discrete frequency points; S3, sequentially extract discrete frequency points from the discrete frequency point combination as the first discrete frequency point, and for each extracted first discrete frequency point, solve... The dimensional matrix equation determines the source intensity of each virtual sound source in the equivalent source model at the first discrete frequency point, until the trend of the source intensity of each virtual sound source changing with the first discrete frequency point satisfies the frequency domain extension condition of the source intensity, where, The total number of virtual sound sources; S4. For each virtual sound source, perform frequency domain extension of the source strength based on the source strength at each first discrete frequency point to determine the source strength of the virtual sound source at each second discrete frequency point, wherein the second discrete frequency point is the discrete frequency point other than the first discrete frequency point in the combination of discrete frequency points. S5, based on the source strength of each virtual sound source at each discrete frequency point in the discrete frequency point combination, calculate the broadband sound field of the target shallow sea area. Step S3 includes the following steps: S31, Obtain discrete frequency point combinations ,in and , , These are the lower and upper limits of the target frequency band, respectively. S32, from all unsolved problems in the discrete frequency point combination In the discrete frequency points of the dimensional matrix equation, the first discrete frequency point is extracted in sequence and marked as the first discrete frequency point; S33, Solve for the first discrete frequency point marked in step S32. The 3D matrix equation is used to obtain the source intensity of each virtual sound source at the first discrete frequency point and extract its real and imaginary parts. S34, update the curves of the real and imaginary source strengths of each virtual sound source as a function of the first discrete frequency point; S35, determine whether the trend of the real part and imaginary part of the source strength of each virtual sound source changing with the first discrete frequency point satisfies the frequency domain extension condition of the source strength. If not, return to step S32. If so, mark all discrete frequency points in the discrete frequency point combination that are not marked as the first discrete frequency point as the second discrete frequency point, and proceed to step S4. The frequency domain extension of the source strength for each virtual sound source in step S4 includes the following steps: S41, Extract the source strength of the virtual sound source at each first discrete frequency point. And based on its actual part and the virtual part Construct source strong real part vector The source strong imaginary part vector ,in, The total number of the first discrete frequency points; S42, respectively for and Perform dynamic mode decomposition operation to obtain The corresponding eigenvalue matrix eigenmode matrix sum coefficient vector ,as well as The corresponding eigenvalue matrix eigenmode matrix sum coefficient vector ; S43, determine the virtual sound source at each of the second discrete frequency points based on the following formula. Source strength : ， in, for of Power of 1 for of Power of 1 This indicates retrieving the real part of the element within the square brackets. This indicates the expression to retrieve the first element within the curly braces. There are elements, among which... The length of each modal eigenvector in the eigenmode matrix is ​​given. It is the imaginary unit.

2. The method for rapid calculation of broadband sound field at complex undulating boundaries in shallow seas according to claim 1, characterized in that, The equivalent source set includes a sea surface equivalent source set for representing sea surface reflections. Used to represent the first equivalent seabed sound source group for seabed transmission. And a second group of equivalent seabed sound sources used to represent seabed reflections. ,in, This is achieved by discretely generating several virtual sound sources on the sea surface and then translating these virtual sound sources upwards along the normal direction from the sea surface. The position vectors of each virtual sound source in the middle are ; and This is achieved by generating two sets of discrete virtual sound sources on an undulating seabed, translating one set of virtual sound sources upwards along the normal direction from the seabed, and translating the other set of virtual sound sources downwards along the normal direction from the seabed. and The position vectors of each virtual sound source are respectively and .

3. The method for rapid calculation of broadband sound field at complex undulating boundaries in shallow seas according to claim 2, characterized in that, , and The spacing between each virtual sound source is of to ,in To determine the upper frequency limit of the target frequency band The lower limit of wavelength in seawater.

4. The method for rapid calculation of broadband sound field at complex undulating boundaries in shallow seas according to claim 1, characterized in that, For each first discrete frequency point, its corresponding The dimensional matrix equation is: ， in, , Equivalent sound source groups , The Green's function matrix of sound pressure at each virtual sound source at each sea surface node. , , Equivalent sound source groups , , The Green's function matrix of the sound pressure of each virtual sound source at each seabed node. , , Equivalent sound source groups , , The displacement Green's function matrix of each virtual sound source at each seabed node. , These are the position vectors of each surface node and seabed node, respectively. Let the wave number in the water be the first discrete frequency point. Let be the wave number in the seabed medium at the first discrete frequency point.

5. The method for rapid calculation of broadband sound field at complex undulating boundaries in shallow seas according to claim 1, characterized in that, The frequency domain extension condition is: The curves of the real and imaginary parts of the source strength of each virtual sound source change with frequency exhibit a periodic repeating structure, and the number of repetitions reaches a preset value.

6. The method for rapid calculation of broadband sound field at complex undulating boundaries in shallow seas according to claim 1, characterized in that, In step S42, the following steps are performed to obtain... , and : A1, using a length of The windows are for the sequence respectively and Perform a matrix 2D transformation operation to obtain the first real part state matrix. Second real part state matrix ; A2, yes Perform truncation rank Truncated SVD decomposition: ， in, , , They are respectively The left singular vector matrix, the singular value matrix, and the right singular vector matrix. for The conjugate transpose of ; A3, Construct a low-dimensional mapping matrix of the real part based on the following formula. : ， in, for The conjugate transpose of the matrix. for The inverse matrix; A4. Solve the eigenvalue problem shown in the following equation to obtain the eigenvalue matrix. and the low-dimensional eigenvector matrix of the real part : ； A5, determined based on the following formula eigenmode matrix : ； A6, the coefficient vector is determined by solving the following formula. : 。 7. The method for rapid calculation of broadband sound field at complex undulating boundaries in shallow seas according to claim 6, characterized in that, In step S42, the following steps are performed to obtain... , and : B1, using a length of The windows are for the sequence respectively and Perform a matrix 2D transformation operation to obtain the first imaginary state matrix. Second imaginary state matrix ; B2, for Perform truncation rank Truncated SVD decomposition: ， in, , , They are respectively The left singular vector matrix, the singular value matrix, and the right singular vector matrix. for The conjugate transpose of ; B3, Construct the low-dimensional mapping matrix of the imaginary part based on the following formula. : ， in, for The conjugate transpose of the matrix. for The inverse matrix; B4. Solve the eigenvalue problem shown in the following equation to obtain the eigenvalue matrix. and the low-dimensional eigenvector matrix of the imaginary part : ； B5 is determined based on the following formula. eigenmode matrix : ； B6, the coefficient vector is determined by solving the following equation. : 。

Images