A method for predicting underwater very low frequency electromagnetic wave propagation characteristics

Abstract

This invention belongs to the field of electromagnetic wave propagation technology and discloses a method for predicting the propagation characteristics of very low frequency (VLF) electromagnetic waves underwater. This method uses the 3D-HIE-FDTD method to numerically model the propagation process of VLF electromagnetic waves originating underwater, enabling relatively accurate prediction of electromagnetic wave propagation characteristics in complex underwater environments. Furthermore, this invention employs a PCR solver to solve the tridiagonal linear equations generated during the implicit update process in parallel, effectively improving solution efficiency and shortening computation time. This meets the practical needs for high-precision and rapid prediction of propagation characteristics in underwater navigation and communication scenarios. This invention not only solves the problems of CFL stability limitations and excessive computational overhead in traditional 3D-FDTD methods for underwater electromagnetic wave propagation but also overcomes the shortcomings of the 3D-HIE-FDTD method, such as low parallelism and difficulty in efficient acceleration and expansion on GPUs when using a serial Thomas solver.

Description

Technical Field

[0001] This invention belongs to the field of electromagnetic wave propagation technology, specifically relating to a method for predicting the propagation characteristics of very low frequency electromagnetic waves underwater. Background Technology

[0002] Very Low Frequency (VLF) electromagnetic waves refer to radio waves with frequencies ranging from 3 kHz to 30 kHz and wavelengths of approximately 10 km to 100 km. VLF electromagnetic waves are characterized by their long wavelength, relatively low propagation loss, good phase stability, and certain penetration ability through seawater. Therefore, they have significant application value in fields such as long-distance communication, navigation, timing, and underwater information transmission. Especially in marine environments, VLF electromagnetic waves can penetrate seawater to a certain depth, providing communication and navigation signal support for submarines, unmanned underwater vehicles, and other underwater equipment, thus possessing unique advantages in the fields of underwater navigation and communication.

[0003] Most existing systems employ land-based transmission, relying on large ground-based transmitters and large-sized transmitting antennas to achieve long-distance radiation. However, these systems generally suffer from problems such as bulky antennas, complex transmitting devices, high construction and maintenance costs, poor deployment flexibility, and insufficient concealment. Currently, for the propagation problem of VLF electromagnetic waves originating underwater, especially for propagation modeling and rapid prediction for underwater navigation and communication applications, there is a lack of mature theoretical models and numerical methods that balance accuracy and efficiency. Due to the characteristics of this type of problem, such as low-frequency long waves, strong attenuation, multi-boundary coupling, and long-distance propagation, its modeling and solution difficulty is significantly higher than that of conventional electromagnetic field problems.

[0004] The Finite-Difference Time-Domain (FDTD) method is widely used in numerical analysis of electromagnetic propagation because it can directly solve Maxwell's equations in the time domain and handles complex media and boundary problems well. However, in underwater VLF navigation and communication problems, due to the significant attenuation in the seawater medium, a small mesh size is usually required for fine discretization to ensure the accuracy of the propagation process and field distribution calculations. Under three-dimensional modeling conditions, if the traditional explicit FDTD method is used, the total number of meshes increases sharply as the mesh size decreases; at the same time, the time step must also satisfy the strict Courant-Friedrichs-Lewy (CFL) stability condition, leading to a significant increase in the number of computation steps. This results in problems such as long computation time, large storage overhead, and low solution efficiency, making it difficult to meet the practical needs of underwater navigation for rapid prediction of propagation characteristics.

[0005] To address the aforementioned issues, the Hybrid Implicit-Explicit Finite-Difference Time-Domain (HIE-FDTD) method provides an effective approach for electromagnetic wave propagation calculations under unidirectional fine-grid conditions. This method introduces an implicit solution mechanism into the update of the electric field components in the x and y directions, effectively relaxing the time step limitations of traditional explicit FDTD methods while maintaining high solution accuracy, thereby reducing the computational burden caused by fine-grid meshing. Particularly in the 3D underwater VLF propagation problem, the HIE-FDTD method is more suitable than the traditional FDTD method for handling complex scenarios with large computational domains, fine meshes, and high stability requirements.

[0006] However, the HIE-FDTD method typically requires solving a tridiagonal linear system of equations during the implicit update process, and existing methods mostly use the Thomas solver for computation. While the Thomas solver is simple to implement, its recursive computation process has strong seriality and insufficient parallel computing capabilities. Especially in large-scale 3D problems and GPU parallel implementations, it is difficult to fully utilize the computational advantages of high-performance hardware, thus limiting further improvements in the overall solution efficiency of the HIE-FDTD method. Summary of the Invention

[0007] The purpose of this invention is to propose a method for predicting the propagation characteristics of very low frequency electromagnetic waves underwater. This method can accurately predict the propagation characteristics of electromagnetic waves in complex underwater environments, while effectively improving the solution efficiency and shortening the calculation time, thereby meeting the practical needs for high-precision and rapid prediction of propagation characteristics in underwater navigation and communication scenarios.

[0008] To achieve the above objectives, the present invention adopts the following technical solution:

[0009] A method for predicting the propagation characteristics of very low frequency electromagnetic waves underwater includes the following steps:

[0010] Step 1. Input the model file and initialize it. The input model file includes the mesh parameters of the computational domain, the excitation source parameters, the electrical parameters of the electromagnetic wave propagation path, the electrical parameters of the free space, the CFS-PML parameters of the absorbing boundary, and the time settings.

[0011] Step 2. Update the entire compute region. Electric field components in the direction Magnetic field components And auxiliary components;

[0012] Step 3. Construct implicit Directional electric field components and Directional electric field components The resulting tridiagonal linear equation system;

[0013] Step 4. Use the PCR solver to solve the tridiagonal linear equations constructed in Step 3 in parallel to obtain the updated electric field components at the current time step. , The value;

[0014] Step 5. Add incentive sources to Electric field components in the direction superior;

[0015] Step 6. Update the entire compute region. and Magnetic field components in direction , And auxiliary components;

[0016] Step 7. Update the current runtime. The value is ;

[0017] Determine if the current running time step is equal to the preset running time step; if yes, proceed to step 8; otherwise, proceed to step 2.

[0018] Step 8. Extract the calculation region Electric field components in the direction The peak field strength is output.

[0019] Furthermore, based on the above-mentioned method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves, this invention also proposes a computer device, which includes a memory and one or more processors.

[0020] The memory stores executable code, and when the processor executes the executable code, it implements the steps of the underwater very low frequency electromagnetic wave propagation characteristic prediction method described above.

[0021] Furthermore, based on the above-mentioned method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves, this invention also proposes a computer-readable storage medium on which a program is stored; when the program is executed by a processor, it is used to implement the steps of the above-mentioned method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves.

[0022] The present invention has the following advantages:

[0023] As described above, this invention relates to a method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves. This method is the first to use the 3D-HIE-FDTD method to numerically model the propagation process of VLF electromagnetic waves with the source located underwater, thereby accurately characterizing the propagation characteristics of electromagnetic waves in complex marine environments. It solves the problems of limited time steps, low computational efficiency, and large storage overhead caused by fine mesh discretization in traditional explicit FDTD methods for underwater VLF propagation problems.

[0024] Meanwhile, this invention addresses the need to solve a tridiagonal linear equation system during the implicit update process of the 3D-HIE-FDTD method. It employs the PCR solver for parallel solving, overcoming the shortcomings of the traditional Thomas solver, which suffers from strong seriality and insufficient parallelism. This effectively improves the computational efficiency of the implicit solution process and further shortens the overall simulation time. While ensuring computational accuracy, this invention is better suited to application scenarios with large computational regions, fine meshes, and complex boundaries in 3D underwater VLF propagation. It provides an efficient and feasible technical solution for high-precision and rapid prediction of the propagation characteristics of very low frequency electromagnetic waves in underwater navigation, making it more suitable for engineering applications and widespread adoption. Attached Figure Description

[0025] Figure 1 This is a flowchart of the underwater very low frequency electromagnetic wave propagation characteristic prediction method in an embodiment of the present invention.

[0026] Figure 2 This is a schematic diagram illustrating the principle of the underwater very low frequency electromagnetic wave propagation characteristic prediction method in an embodiment of the present invention.

[0027] Figure 3 To compare the electric field components along a smooth sea surface using the method of this invention with conventional methods A schematic diagram of the prediction results.

[0028] Figure 4 To apply the method of this invention to the electric field components when there are wave fluctuations along the propagation path. A schematic diagram of the prediction results.

[0029] in, Figure 4 In the diagram, (a) represents the electric field components along the propagation path. Curve of variation with distance, Figure 4 (b) in the figure represents the wave height variation curve along the propagation path. Detailed Implementation

[0030] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0031] Example 1

[0032] This embodiment specifically discloses a three-dimensional HIE-FDTD (3D-HIE-FDTD) underwater very low frequency electromagnetic wave propagation characteristic prediction method based on a parallel cyclic reduction (PCR) solver. This method uses an improved 3D-HIE-FDTD algorithm to predict the impact of propagation path changes caused by wave fluctuations on the field strength. It not only solves the problem of excessive computational overhead in underwater electromagnetic wave propagation due to CFL stability limitations of traditional 3D-FDTD methods, but also overcomes the shortcomings of low parallelism and difficulty in efficient acceleration and expansion on GPUs when using a serial Thomas solver in the 3D-HIE-FDTD algorithm. It can provide theoretical basis and technical support for improving the service accuracy of underwater very low frequency radio navigation and timing systems.

[0033] The general idea of ​​the method of this invention is as follows: First, input the model file; second, update the electric field components of the explicit part of the entire computational domain. and magnetic field components And auxiliary components; then construct implicit... and Directional electric field components , The system of tridiagonal linear equations, including auxiliary components, is solved in parallel using a PCR solver; then the entire computational domain is updated. and Magnetic field components in direction and And auxiliary components; finally, repeat the above process until the termination condition is met, and extract the data from the specified observation point or path when the calculation ends. Electric field components in the direction The peak field strength data is obtained and the prediction results are output.

[0034] This invention employs the 3D-HIE-FDTD method to numerically model the propagation process of very low frequency electromagnetic waves originating underwater, enabling accurate prediction of electromagnetic wave propagation characteristics in complex underwater environments. Simultaneously, by using a PCR solver to solve the tridiagonal linear equations generated during the implicit update process in parallel, it effectively improves solution efficiency and shortens computation time, thus meeting the practical needs for high-precision and rapid prediction of propagation characteristics in underwater navigation and communication scenarios.

[0035] Based on the above inventive concept, the underwater very low frequency electromagnetic wave propagation characteristic prediction method proposed in this invention will be described in detail below with reference to the accompanying drawings. Figure 1 and Figure 2 As shown, the underwater very low frequency electromagnetic wave propagation characteristic prediction method includes the following steps:

[0036] Step 1. Input the model file and initialize it. The input model file includes the mesh parameters of the computational domain, the excitation source parameters, the electrical parameters of the electromagnetic wave propagation path, the electrical parameters of free space, the CFS-PML parameters of the absorbing boundary, and the time settings.

[0037] In this embodiment, step 1 specifically includes:

[0038] In a three-dimensional rectangular coordinate system In the definition, the computational region size is as follows: ,in Horizontal coordinate Number of grids in the direction, Vertical axis Number of grids in the direction, height coordinates Number of grid cells in the direction; , , The meshing step size in each direction is set to... , , .

[0039] The excitation source parameters include the excitation source location, excitation source loading method, signal frequency, and signal amplitude.

[0040] Among them, the location of the excitation source is determined by , , The starting and ending coordinates in three directions are determined; the excitation source loading methods include hard source loading and soft source loading; in the source model, the excitation signal adopts a sine wave, and the excitation source is equivalent to a horizontal electric dipole, with its current I and dipole length dl preset.

[0041] The electrical parameters of the radio wave propagation path are the same as the electrical parameters of the sea level, including the relative permittivity of the seawater. and conductivity .

[0042] The electrical parameters of free space include the air dielectric constant. and air permeability .

[0043] The CFS-PML parameter of the absorbing boundary includes the number of absorbing boundary layers. Conductivity distribution parameters Elongation coefficient and attenuation coefficient , , and The distribution along the absorption boundary region varies according to a preset pattern, used to jointly regulate the impedance matching characteristics, attenuation intensity, and absorption stability of electromagnetic waves within the boundary region, thereby reducing non-physical reflections at the computational boundary. Specifically:

[0044] (1)

[0045] (2)

[0046] (3)

[0047] in, Indicates the order of parameter distribution; , Indicates the operating wavelength; Indicates the coordinate direction variable. The coordinates represent the starting position of the absorption layer; Indicates the thickness of the absorption boundary; It is an integer. It is a constant.

[0048] Time settings include the current running time step. Time step and preset running time steps .

[0049] Parameters initialized to 0 include:

[0050] , , Electric field components in the direction , , .

[0051] , , Magnetic field components in direction , , .

[0052] Electric field auxiliary component and magnetic field auxiliary components ,in .

[0053] Electric field component calculation coefficients , Magnetic field component calculation coefficients , Its calculation expression is:

[0054] (4)

[0055] (5)

[0056] (6)

[0057] (7)

[0058] in, = ; This represents the relative permeability.

[0059] Parameters of auxiliary variables in the CFS-PML framework , , ,in .

[0060] Current running time step .

[0061] This step, by initializing the medium parameters, boundary parameters, and update coefficients, provides a consistent numerical basis for subsequent implicit and explicit hybrid updates and provides a prerequisite for accurately reflecting the impact of complex propagation paths on the field strength. Pre-setting the absorbing boundary CFS-PML parameters can effectively suppress the interference of truncated boundary reflections on the calculation results and improve the stability and reliability of large-area propagation calculations.

[0062] Step 2. Update the entire compute region. Electric field components in the direction Magnetic field components And auxiliary components.

[0063] In this embodiment, step 2 specifically includes:

[0064] Update the entire computing area Electric field components in the direction Magnetic field components and electric field auxiliary components and Magnetic field auxiliary component and The specific calculation formula is as follows:

[0065] (8)

[0066] (9)

[0067] (10)

[0068] (11)

[0069] (12)

[0070] (13)

[0071] in, , , For three-dimensional discrete space indexing, express At time step ,space The value at that location, Indicates time difference, , , Indicates spatial difference, Including electric field components , , Magnetic field components , , Electric field auxiliary component and magnetic field auxiliary components .

[0072] express In the x-direction The values ​​at each grid point express In the y-direction The values ​​at each grid point express In the z-direction The values ​​at each grid point; Including parameters of auxiliary variables under the CFS-PML framework , , .

[0073] This step first completes the explicit field component update, allowing the 3D-HIE-FDTD method to retain the advantages of time-domain recursion and discrete accuracy of the 3D-FDTD method, while leaving only the field updates in the necessary directions to subsequent implicit solution processing. Compared to the fully explicit 3D-FDTD method, this invention reduces the computational burden caused by the strict CFL condition limiting the time step under small grid conditions through a combination of explicit and implicit update methods, thus providing a more efficient numerical solution basis for large-scale, fine-grid underwater very low frequency propagation problems.

[0074] Step 3. Construct implicit Directional electric field components and Directional electric field components The resulting tridiagonal linear equation system.

[0075] In this embodiment, step 3 specifically includes:

[0076] Constructing implicit Directional electric field components The resulting tridiagonal linear equation system is as follows:

[0077] (14)

[0078] Constructing implicit Directional electric field components The resulting tridiagonal linear equation system is as follows:

[0079] (15)

[0080] Constructing implicit and Direction with respect to auxiliary components , , , The tridiagonal linear equation system is as follows:

[0081] (16)

[0082] (17)

[0083] (18)

[0084] (19)

[0085] The method of this invention transforms the implicit coupling update into a standard tridiagonal linear equation system, which can convert the time-progression problem, which was originally constrained by stability conditions, into an algebraic problem with regular structure and suitable for parallel processing. This approach not only maintains the modeling capability of the 3D-HIE-FDTD method in complex cross-media environments, but also lays the foundation for the subsequent introduction of the PCR solver to achieve efficient parallel acceleration.

[0086] Step 4. Use the PCR solver to solve the tridiagonal linear equations constructed in Step 3 in parallel to obtain the updated electric field components at the current time step. , The value.

[0087] In this embodiment, step 4 specifically includes:

[0088] Step 4.1. Establish the coefficient representation of the tridiagonal linear equation system.

[0089] Acquisition scale is tridiagonal linear equations , No. The equations are written as follows:

[0090] (20)

[0091] in, , .

[0092] The number of equations; The index number of the equation; Let be the vector of unknowns to be solved. , , They represent the first The, the The, the One unknown quantity; The coefficient is the lower diagonal coefficient, corresponding to coefficient ; The main diagonal coefficients correspond to coefficient ; The coefficient of the upper diagonal is the corresponding coefficient. coefficient ; Given the terms on the right-hand side, written in vector form as follows: .

[0093] and If it does not exist, the boundary is handled in the following equivalent way:

[0094] right The equation is: , equivalent to .

[0095] right The equation is: , equivalent to .

[0096] Step 4.2. Initialize PCR iteration level and span parameters.

[0097] Set the number of specification layers .

[0098] Define the span for each layer , and the first The system of equations for the layer is denoted as ;in, , , , They represent the first After the first layer specification The coefficients of the lower diagonal, the coefficients of the main diagonal, the coefficients of the upper diagonal, and the known terms on the right-hand side of each equation.

[0099] The initial layer satisfies: , , , .

[0100] In the After the layer reduction was completed, the original system was divided into There are three uncoupled subsystems; the number of unknowns in each subsystem is... Overall Completed within the specified layer.

[0101] Step 4.3. Parallel reduction and generation of half-size subsystems.

[0102] The first After the first layer specification The equation is linearly combined with the equations on its left and right sides to obtain the equation of the second degree. New coefficients for the layer:

[0103] (twenty one)

[0104] (twenty two)

[0105] (twenty three)

[0106] (twenty four)

[0107] (25)

[0108] in, , , , They represent the first After the first layer specification The lower diagonal coefficients, main diagonal coefficients, upper diagonal coefficients, and known terms on the right-hand side of each equation; , , , They represent the first Layer-specification and the first The lower diagonal coefficient, main diagonal coefficient, upper diagonal coefficient, and known terms on the right side of the equations adjacent to the left side of the equations; , , , They represent the first Layer-specification and the first The lower diagonal coefficient, main diagonal coefficient, upper diagonal coefficient, and known terms on the right side of the equations adjacent to each other on the right. , Indicates the first Elimination coefficients used in layer reduction to eliminate the influence of left and right adjacent unknowns.

[0109] After completing the first After the layer update, two independent tridiagonal sub-equations are obtained, each of which is about half the size of the original equation set.

[0110] Then, taking each sub-equation as the processing object, the parallel cyclic reduction process shown in formulas (21) to (25) is used to update and calculate the coefficients of the next layer, so that the coupling span of the equations doubles and is further split into more smaller sub-equation sets; the parallel cyclic reduction process is repeated on the sub-equation sets until each subsystem is reduced to a 2×2 tridiagonal subsystem containing only two unknowns.

[0111] Step 4.4. Solve all 2×2 tridiagonal subsystems in parallel.

[0112] Only two unknowns remain , At that time, solve the 2×2 linear equation system in parallel:

[0113] (26)

[0114] in, and They represent the first After the first layer specification The equation and the first The main diagonal coefficients of the equations; Indicates the first After the first layer specification The upper diagonal coefficients of the equations; Indicates the first After the first layer specification The lower diagonal coefficients of the equations; and They represent the first After the first layer specification The equation and the first The right-hand side of the equation has known terms.

[0115] The unknowns are solved synchronously and in parallel for all 2×2 tridiagonal subsystems, yielding the final solution vector. Each component does not require a back-substitution stage.

[0116] Step 4.5. Backfill and combine the unknowns obtained from all 2×2 tridiagonal subsystems according to their original index positions to complete the solution of the tridiagonal equations and obtain the updated electric field components at the current time step. , The value.

[0117] Existing HIE-FDTD implicit updates often employ the Thomas algorithm to solve tridiagonal equations. However, the Thomas algorithm suffers from significant serial characteristics and poor scalability on parallel hardware platforms such as GPUs. In contrast, this invention uses a PCR solver to transform the original tridiagonal system into a multi-layered parallel reduction problem, significantly improving the solver's parallelism and mitigating the efficiency bottleneck caused by serial recursion. This allows for more full utilization of the computational power of parallel computing platforms.

[0118] Step 5. Add incentive sources to Electric field components in the direction superior.

[0119] In this embodiment, step 5 specifically includes:

[0120] The added field source, i.e., the excitation source, is a sinusoidal signal, consisting of electric field components. Excitation, current waveform of the field source Represented as:

[0121] (27)

[0122] in, .

[0123] when hour, ;when hour, .

[0124] Step 6. Update the entire compute region. and Magnetic field components in direction , And auxiliary components.

[0125] In this embodiment, step 6 specifically includes:

[0126] Update the entire computing region Magnetic field components in direction , Magnetic field components in direction and auxiliary components , , , The specific calculation formula is as follows:

[0127] (28)

[0128] (29)

[0129] (30)

[0130] (31)

[0131] (32)

[0132] (33)

[0133] The method of this invention completes the unified update of the magnetic field components in step 6, forming a complete iterative closed loop of explicit field update - implicit equation construction and parallel solution - magnetic field update. This closed loop not only retains the stable discrete structure of the Yee-type time-domain advancement method in electromagnetic propagation solution, but also reduces the time advancement pressure of the traditional 3D-FDTD method in small grid and large area scenes by combining implicit and explicit methods, thus ensuring both solution accuracy and computational efficiency.

[0134] Step 7. Update the current runtime. The value is Then determine whether the current running time step is equal to the preset running time step; if yes, go to step 8; otherwise, go to step 2.

[0135] Step 8. Extract the calculation region Electric field components in the direction The peak field strength is output.

[0136] In this embodiment, step 8 specifically includes:

[0137] The peak radiated power at the envelope is time The calculation formula is:

[0138] (34)

[0139] in, Indicates the level of the extracted receiving point signal. Directional electric field components Peak value It is a constant value.

[0140] Extract the calculation area using formula (34) Electric field components in the direction The peak field strength is output.

[0141] This step converts the numerical calculation results into propagation characteristic indicators that can be directly used for engineering analysis. Since the entire solution process is completed within a three-dimensional cross-medium environment, considering the influence of ocean wave fluctuations and employing an efficient implicit parallel solution framework, the final field strength results can more realistically reflect the electromagnetic wave propagation laws in complex underwater propagation scenarios and provide stronger technical support for system service accuracy analysis.

[0142] In addition, to verify the effectiveness of the method proposed in this invention, the following specific experiments are also provided:

[0143] Experiment 1 shows the predicted field strength along a path on a flat sea surface.

[0144] The size of the computation region is: The traditional 3D-FDTD mesh generation parameters are: , The time steps are 3m, 9m, 15m, and 21m respectively. for The 3D-HIE-FDTD mesh generation parameters of the method of this invention are as follows: , The time steps are 3m, 9m, 15m, and 21m respectively. for The relative permittivity of seawater electrical conductivity .

[0145] Figure 3 The method of this invention and the traditional 3D-FDTD method are compared on the propagation path of the electric field component on a smooth sea surface with different electrical parameters. Prediction results. (By...) Figure 3 It can be seen that the errors of the two methods are very small, proving that the method of the present invention is suitable for field strength prediction.

[0146] Experiment 2 shows the predicted field strength along the propagation path of ocean wave undulations.

[0147] The size of the computation region is: The mesh size is , Time step for The relative permittivity of seawater electrical conductivity .

[0148] The wave undulation model calculates the wave height, wave number, wavelength, period, and angular frequency of the waves based on the input wind speed, and establishes a radial propagation distance model centered on a preset source point. Then, by combining the exponential decay function and the cosine wave function, the instantaneous undulation height of the sea surface at each location on the two-dimensional sea area grid is obtained. Finally, the continuous wave height is discretized and mapped to the seawater filling termination layer in the three-dimensional grid, realizing the construction of marine medium parameters with dynamic wave boundaries.

[0149] Figure 4 The solid black line within (a) in the diagram represents the conditions of undulating sea surface waves. Electric field distribution, the red dashed line represents the electric field distribution under flat sea surface conditions. Electric field distribution; Figure 4 Figure (b) shows the corresponding wave height distribution. Figure 4 It can be seen that, compared with a smooth sea surface, the electric field distribution under undulating sea surface conditions exhibits more obvious local fluctuations and non-uniform variations, especially in the propagation region after the main peak, where the electric field curve shows continuous oscillations. This result indicates that undulating waves alter the sea surface boundary conditions and the reflection, scattering, and interference relationships during electromagnetic wave propagation, thus causing a fluctuation in electric field amplitude of approximately 14 dB compared to a smooth sea surface.

[0150] The method of this invention overcomes the shortcomings of traditional 3D-FDTD methods, which are limited by the stability condition of CFL and have excessive computational overhead in large-area cross-media propagation modeling. For the tridiagonal matrix equation system formed after implicit direction discretization, it uses the PCR solver to replace the serial Thomas solver for parallel acceleration, thereby improving the parallelism and computational efficiency of the implicit solution process and enhancing the acceleration and scalability of the algorithm on the GPU platform.

[0151] Example 2

[0152] This embodiment 2 describes a computer device that includes a memory and one or more processors.

[0153] The memory stores executable code, which, when executed by the processor, is used to implement the steps of the underwater very low frequency electromagnetic wave propagation characteristic prediction method in Embodiment 1 above.

[0154] In this embodiment, the computer device can be any device or apparatus with data processing capabilities, and will not be described in detail here.

[0155] Example 3

[0156] This embodiment 3 describes a computer-readable storage medium storing a program that, when executed by a processor, implements the steps of a method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves.

[0157] The computer-readable storage medium can be an internal storage unit of any device or apparatus with data processing capabilities, such as a hard disk or memory, or an external storage device of any device with data processing capabilities, such as a plug-in hard disk, smart media card (SMC), SD card, flash card, etc.

[0158] Of course, the above description is only a preferred embodiment of the present invention. The present invention is not limited to the above-described embodiments. It should be noted that any equivalent substitutions or obvious modifications made by those skilled in the art under the guidance of this specification fall within the scope of this specification and should be protected by the present invention.

Claims

1. A method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves, characterized in that, Includes the following steps: Step 1. Input the model file and initialize it. The input model file includes the mesh parameters of the computational domain, the excitation source parameters, the electrical parameters of the electromagnetic wave propagation path, the electrical parameters of the free space, the CFS-PML parameters of the absorbing boundary, and the time settings. Step 2. Update the entire compute region. Electric field components in the direction Magnetic field components And auxiliary components; Step 3. Construct implicit Directional electric field components and Directional electric field components The resulting tridiagonal linear equation system; Step 4. Use the PCR solver to solve the tridiagonal linear equations constructed in Step 3 in parallel to obtain the updated electric field components at the current time step. , The value; Step 5. Add incentive sources to Electric field components in the direction superior; Step 6. Update the entire compute region. and Magnetic field components in direction , And auxiliary components; Step 7. Update the current runtime. The value is ; Determine if the current running time step is equal to the preset running time step; if so, proceed to step 8. Otherwise, proceed to step 2; Step 8. Extract the calculation region Electric field components in the direction The peak field strength is output.

2. The method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves according to claim 1, characterized in that, Step 1 specifically involves: In a three-dimensional rectangular coordinate system In the definition, the computational region size is as follows: ,in Horizontal coordinate Number of grids in the direction, Vertical axis Number of grids in the direction, height coordinates Number of grid cells in the direction; , , The meshing step size in each direction is set to... , , ; Excitation source parameters include excitation source location, excitation source loading method, signal frequency, and signal amplitude; Among them, the location of the excitation source is determined by , , The starting and ending coordinates in three directions are determined; the excitation source loading methods include hard source loading and soft source loading; in the source model, the excitation signal adopts a sine wave, and the excitation source is equivalent to a horizontal electric dipole, with its current I and dipole length dl preset; The electrical parameters of the radio wave propagation path are the same as the electrical parameters of the sea level, including the relative permittivity of the seawater. and conductivity ; The electrical parameters of free space include the air dielectric constant. and air permeability ; The CFS-PML parameter of the absorbing boundary includes the number of absorbing boundary layers. Conductivity distribution parameters Elongation coefficient and attenuation coefficient ,in: (1) (2) (3) in, Indicates the order of parameter distribution; , Indicates the operating wavelength; Indicates the coordinate direction variable. The coordinates represent the starting position of the absorption layer; Indicates the thickness of the absorption boundary; It is an integer. It is a constant; Time settings include the current running time step. Time step and preset running time steps ; Parameters initialized to 0 include: , , Electric field components in the direction , , ; , , Magnetic field components in direction , , ; Electric field auxiliary component and magnetic field auxiliary components ,in ; Electric field component calculation coefficients , Magnetic field component calculation coefficients , Its calculation expression is: (4) (5) (6) (7) in, = ; Indicates relative permeability; Parameters of auxiliary variables in the CFS-PML framework , , ,in ; Current running time step .

3. The method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves according to claim 2, characterized in that, Step 2 specifically involves: Update the entire computing area Electric field components in the direction Magnetic field components and electric field auxiliary components and Magnetic field auxiliary component and The specific calculation formula is as follows: (8) (9) (10) (11) (12) (13) in, , , For three-dimensional discrete space indexing, express At time step ,space The value at that location, Indicates time difference, , , Indicates spatial difference, Including electric field components , , Magnetic field components , , Electric field auxiliary component and magnetic field auxiliary components ; express In the x-direction The values ​​at each grid point express In the y-direction The values ​​at each grid point express In the z-direction The values ​​at each grid point; Including parameters of auxiliary variables under the CFS-PML framework , , .

4. The method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves according to claim 3, characterized in that, Step 3 specifically involves: Constructing implicit Directional electric field components The resulting tridiagonal linear equation system is as follows: (14) Constructing implicit Directional electric field components The resulting tridiagonal linear equation system is as follows: (15) Constructing implicit and Direction with respect to auxiliary components , , , The tridiagonal linear equation system is as follows: (16) (17) (18) (19)。 5. The method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves according to claim 4, characterized in that, Step 4 specifically involves: Step 4.

1. Establish the coefficient representation of the tridiagonal linear equation system; Acquisition scale is tridiagonal linear equations , No. The equations are written as follows: (20) in, , ; The number of equations; This is the index number of the equation; Let be the vector of unknowns to be solved. , , They represent the first The, the The, the One unknown quantity; The coefficient is the lower diagonal coefficient, corresponding to coefficient ; The main diagonal coefficients correspond to coefficient ; The coefficient of the upper diagonal corresponds to coefficient ; Given the terms on the right-hand side, written in vector form as follows: ; and If it does not exist, the boundary is handled in the following equivalent way: The equation is: , equivalent to ; The equation is: , equivalent to ; Step 4.

2. Initialize PCR iteration levels and span parameters; Set the number of specification layers ; Define the span for each layer , and will the The system of equations for the layer is denoted as ;in, , , , They represent the first After the first layer specification The lower diagonal coefficients, main diagonal coefficients, upper diagonal coefficients, and known terms on the right-hand side of each equation; The initial layer satisfies , , , ; In the After the layer reduction was completed, the original system was divided into There are three uncoupled subsystems; the number of unknowns in each subsystem is... Overall Completed within the specified layer; Step 4.

3. Parallel reduction and generation of half-size subsystems; The first After the first layer specification The equation is linearly combined with the equations on its left and right sides to obtain the equation of the second degree. New coefficients for the layer: (21) (22) (23) (24) (25) in, , , , They represent the first After the first layer specification The lower diagonal coefficients, main diagonal coefficients, upper diagonal coefficients, and known terms on the right-hand side of each equation; , , , They represent the first Layer-specification and the first The lower diagonal coefficient, main diagonal coefficient, upper diagonal coefficient, and known terms on the right side of the equations adjacent to the left side of the equations; , , , They represent the first Layer-specification and the first The lower diagonal coefficient, main diagonal coefficient, upper diagonal coefficient, and known terms on the right side of the equations adjacent to each other on the right. , Indicates the first Elimination coefficients used in layer reduction to eliminate the influence of left and right adjacent unknowns; After completing the first After the layer update, two independent tridiagonal sub-equations are obtained; Then, taking each sub-equation system as the processing object, the parallel cyclic reduction process shown in formulas (21) to (25) is used to update and calculate the coefficients of the next level. The parallel cyclic reduction is repeated on the sub-equation system until each subsystem is reduced to a 2×2 tridiagonal subsystem containing only two unknowns. Step 4.

4. Solve all 2×2 tridiagonal subsystems in parallel; Only two unknowns remain , At that time, solve the 2×2 linear equation system in parallel: (26) in, and They represent the first After the first layer specification The equation and the first The main diagonal coefficients of the equations; Indicates the first After the first layer specification The upper diagonal coefficients of the equations; Indicates the first After the first layer specification The lower diagonal coefficients of the equations; and They represent the first After the first layer specification The equation and the first The right-hand side of each equation has known terms; The unknowns are solved synchronously and in parallel for all 2×2 tridiagonal subsystems, yielding the final solution vector. Each component; Step 4.

5. Backfill and combine the unknowns obtained from all 2×2 tridiagonal subsystems according to their original index positions to complete the solution of the tridiagonal equations and obtain the updated electric field components at the current time step. , The value.

6. The method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves according to claim 5, characterized in that, Step 5 specifically involves: The added field source, i.e., the excitation source, is a sinusoidal signal, consisting of electric field components. Excitation, current waveform of the field source Represented as: (27) in, ; when hour, ;when hour, .

7. The method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves according to claim 6, characterized in that, Step 6 specifically involves: Update the entire computing region Magnetic field components in direction , Magnetic field components in direction and auxiliary components , , , The specific calculation formula is as follows: (28) (29) (30) (31) (32) (33)。 8. The method for predicting the propagation characteristics of underwater very low frequency electromagnetic waves according to claim 7, characterized in that, Step 8 specifically involves: Reduced to the peak envelope radiated power is time The calculation formula is: (34) in, Indicates the level of the extracted receiving point signal. Directional electric field components Peak value It is a constant value; Extract the calculation area using formula (34) Electric field components in the direction The peak field strength is output.

9. A computer device comprising a memory and one or more processors, wherein the memory stores executable code, characterized in that, When the processor executes the executable code, it implements the steps of the underwater very low frequency electromagnetic wave propagation characteristic prediction method as described in any one of claims 1 to 8.

10. A computer-readable storage medium having a program stored thereon, characterized in that, When the program is executed by the processor, it implements the steps of the underwater very low frequency electromagnetic wave propagation characteristic prediction method as described in any one of claims 1 to 8.

Images