Prediction method of low-frequency ground wave propagation characteristics over undulating terrain based on FD algorithm

Abstract

The present invention belongs to the field of radio wave propagation technology, and specifically discloses a method for predicting low-frequency ground wave propagation characteristics in undulating terrain paths based on the FD algorithm. The method first establishes a 3DPE model under irregular terrain, then constructs a discretized matrix equation, then merges boundary impedance conditions based on a unilateral approximation method, then iteratively solves the field matrix, and finally uses the field results in the entire calculation area to calculate the magnetic field intensity and quadratic delay of radio wave propagation. The method of the present invention effectively solves the problem of inaccurate characterization of low-frequency ground wave propagation characteristics in lateral irregular terrain modeling using the 2DPE method, and at the same time overcomes the limitation of the traditional 3DPE model in insufficient low-frequency ground wave propagation prediction accuracy in complex undulating terrain environments. It can provide theoretical guidance for the performance optimization of low-frequency radio navigation and timing systems, and further improve the accuracy of navigation positioning and timing services in complex terrain areas.

Description

Technical Field

[0001] The present invention belongs to the technical field of radio wave propagation, and in particular relates to a method for predicting low-frequency ground wave propagation characteristics of a path over undulating terrain based on an FD algorithm. Background Art

[0002] Among the key means currently used to obtain precise time and location information, low-frequency radio navigation / timing systems have become an important component of the Positioning, Navigation and Timing (PNT) system due to their outstanding advantages such as low propagation loss, strong anti-interference ability, and long propagation distance.

[0003] However, the existing two-dimensional propagation delay prediction theory cannot fully reflect the interaction between low-frequency ground waves and complex terrain environments. Among them, the parabolic equation (PE) method, which can achieve both high precision and fast speed, has problems such as unstable Brewster angle and inaccurate terrain modeling in its three-dimensional extension, which leads to inaccurate prediction results of low-frequency ground wave propagation characteristics.

[0004] Accurately modeling complex terrain and comprehensively considering the impact of changes in earth conductivity and undulating terrain on low-frequency ground waves along the propagation path effectively address the problem that two-dimensional models cannot account for phenomena such as reflection and diffraction caused by lateral wave propagation. At the same time, solving the three-dimensional parabolic equation (3DPE) to more accurately predict the propagation delay of low-frequency ground waves when propagating over irregular terrain is a key issue in improving the accuracy of low-frequency wireless navigation / timing systems.

[0005] Therefore, it is necessary to propose a method that can accurately predict the propagation characteristics of low-frequency ground waves in paths over undulating terrain. Summary of the Invention

[0006] The purpose of the present invention is to propose a method for predicting the propagation characteristics of low-frequency ground waves on undulating terrain paths based on the FD algorithm. This method uses a finite difference (FD) algorithm to solve the 3DPE (3D-Conform Map Model Parabolic Equation, 3D-CMM-PE) that introduces a shift transformation model. This method can accurately predict the propagation characteristics of low-frequency ground waves on undulating terrain.

[0007] In order to achieve the above-mentioned purpose, the present invention adopts the following technical solutions:

[0008] The method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm includes the following steps:

[0009] Step 1. Input the model file. The content of the input model file includes the grid parameters of the calculation area, the electrical parameters of the radio wave propagation path, and the excitation source parameters.

[0010] Step 2. Calculate the initial field using the flat ground formula and use the initial field as the source field for 3DPE.

[0011] Step 3. Establish a 3DPE model under irregular terrain;

[0012] Step 4. Construct the discretized matrix equation for solving the 3DPE model;

[0013] Step 5. Merge the boundary impedance conditions of the 3DPE model based on the unilateral approximation method;

[0014] Step 6. Iteratively solve the field matrix based on the discretized matrix equation and the boundary impedance conditions of the 3DPE model to obtain the field distribution results in the entire calculation area;

[0015] Step 7. Using the field results in the entire calculation area, calculate the magnetic field strength of the low-frequency ground wave propagating along the path of the undulating terrain. And the secondary delay generated .

[0016] In addition, based on the above-mentioned method for predicting low-frequency ground wave propagation characteristics of a path over undulating terrain based on the FD algorithm, the present invention also proposes a computer device, which includes a memory and one or more processors;

[0017] The memory stores executable codes, and when the processor executes the executable codes, it is used to implement the steps of the above-mentioned method for predicting low-frequency ground wave propagation characteristics of a path over undulating terrain based on the FD algorithm.

[0018] The present invention has the following advantages:

[0019] As described above, the present invention describes a method for predicting low-frequency ground wave propagation characteristics over undulating terrain paths based on an FD algorithm. This method is an FD algorithm capable of solving 3D-CMM-PE. It introduces a displacement transformation model (CMM) suitable for diffraction propagation into terrain modeling, establishing a 3D-CMM-PE form for complex terrain, reducing the impact of dispersion errors while predicting the lateral diffraction effects caused by the interaction between irregular terrain and low-frequency ground waves. The present method not only effectively solves the problem of inaccurate characterization of low-frequency ground wave propagation characteristics in modeling laterally irregular terrain using the two-dimensional parabolic equation (2DPE) method, but also overcomes the limitations of the traditional 3DPE model, which lacks sufficient accuracy in predicting low-frequency ground wave propagation in complex undulating terrain environments, thereby improving the accuracy of predicting low-frequency ground wave propagation characteristics over undulating terrain. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 This is a flow chart of the method for predicting low-frequency ground wave propagation characteristics over undulating terrain paths based on the FD algorithm in the present invention.

[0021] Figure 2 Schematic diagram of the principle of the method for predicting low-frequency ground wave propagation characteristics of undulating terrain paths based on the FD algorithm in the present invention.

[0022] Figure 3 The magnetic field intensity predicted on a plane by the method of the present invention and the traditional method is Schematic diagram of the results.

[0023] Figure 4 In order to use the method of the present invention and the traditional method Quadratic delay predicted on the plane Schematic diagram of the results.

[0024] Figure 5 is the low-frequency ground wave magnetic field intensity when there are two Gaussian mountains on the propagation path in the present invention Schematic diagram of the prediction results.

[0025] Figure 6 is the secondary delay of the low-frequency ground wave when there are two Gaussian mountains on the propagation path in the present invention. Schematic diagram of the prediction results.

[0026] Figure 7 for Figure 4 A partial enlarged view of point A in the middle. DETAILED DESCRIPTION

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

[0028] The present invention provides a method for predicting the propagation characteristics of low-frequency ground waves over undulating terrain based on the FD algorithm. This method improves the problem that the 2DPE method cannot consider the impact of lateral irregular terrain on low-frequency ground waves, and overcomes the problem that the traditional 3DPE method is insufficient in solving and predicting the propagation characteristics of low-frequency ground waves over undulating terrain. The overall process is as follows:

[0029] First, a 3DPE model under irregular terrain is established. Then, a discretized matrix equation is constructed. Then, the boundary impedance condition is set based on the unilateral approximation method. Then, the field matrix is ​​solved iteratively. Finally, the magnetic field intensity of the radio wave propagation is obtained using the field results in the entire calculation area. and quadratic delay .

[0030] Based on the above invention concept, the method for predicting low-frequency ground wave propagation characteristics of undulating terrain paths based on the FD algorithm proposed by the present invention is described in detail below with reference to the accompanying drawings. Figure 1 As shown, the method of the present invention comprises the following steps:

[0031] Step 1. Input the model file. The content of the input model file includes the grid parameters of the calculation area, the electrical parameters of the radio wave propagation path, and the excitation source parameters.

[0032] The grid parameters of the simulation area, i.e. the calculation area: in the three-dimensional rectangular coordinate system In the calculation area, the size is defined as , where the number of grids is set to 、 、 , N x is the number of grids in the x-direction of the horizontal coordinate, N y is the number of grids in the y-direction of the longitudinal coordinate, N z is the number of grids in the height coordinate z direction; 、 、 is the meshing step size. Figure 2 As shown, O represents the location of the origin of the calculation area coordinates. Indicates the coordinate origin position after coordinate transformation.

[0033] The simulation path is the electrical parameters of the radio wave propagation path: the electrical parameters of the radio wave propagation path include the relative dielectric constant of free space and conductivity , relative dielectric constant in mountainous areas and conductivity .

[0034] Source parameters: The excitation source is a vertical electric dipole, and the current of the vertical electric dipole is , the charge separation of the vertical electric dipole is , the vertical electric dipole is at a height of .

[0035] Step 2. Use the flat ground formula to calculate the initial field and use the initial field as the source field of 3DPE.

[0036] In the rectangular coordinate system Next, let the time harmonic factor be ,in is the imaginary unit, represents the angular velocity, Represents time. Use the flat ground formula to solve the initial field of 3DPE, the ground wave magnetic field component Expressed as:

[0037] (1)

[0038] in, represents the great circle distance; and are the wave numbers in vacuum and on the ground respectively; From source point to observation point The straight-line distance at is the straight-line distance from the source's mirror image to the observation point; is the Fresnel integral, is an intermediate variable.

[0039] Great Circle Distance Expressed as:

[0040] .

[0041] Fresnel integral Defined as:

[0042] .

[0043] Intermediate variables Expressed as:

[0044] .

[0045] in, is the observation point location, Indicates the vertical height coordinate.

[0046] For the preset initial distance According to formula (1), the initial field of low-frequency ground wave is obtained for:

[0047] (2)

[0048] in, express The ground wave magnetic field component at As the source field of 3DPE.

[0049] Step 3. Build a 3DPE model under irregular terrain.

[0050] The three-dimensional parabola equation on a flat surface is expressed as follows:

[0051] (3)

[0052] in, represents the transverse Laplacian operator, ; is the target field to be solved, is the refractive index of the medium, . is the dielectric constant of the propagation medium, such as when propagating in free space, i.e. air. . is the relative magnetic permeability of the propagation medium, such as when propagating in free space, i.e. air. , Indicates the relative magnetic permeability when propagating in air.

[0053] The coordinate transformation model CMM is used to transform the original undulating terrain into a flat ground, and the height and distance variables are set as follows:

[0054] (4)

[0055] in, is the terrain function; after coordinate transformation, the spatial step size satisfies , , .

[0056] Substituting formula (4) into formula (3) to simplify the equation, we can obtain the 3DPE model under irregular terrain as shown in formula (5):

[0057] (5)

[0058] To simplify the formula, the terrain function Simplified to .

[0059] The 3DPE model under irregular terrain in formula (5) takes into account the complex terrain conditions with obvious lateral terrain changes encountered when dealing with three-dimensional irregular terrain, retains the partial derivative of the terrain in the y direction in the 3DPE model, and ensures the prediction accuracy of the low-frequency ground wave propagation characteristics.

[0060] Step 4. Construct the discretized matrix equation for solving the 3DPE model.

[0061] First, the calculation area is along the propagation direction Discretize the field into points, the waveguide consists of a continuous The horizontal plane is composed of Indicates the Points in The coordinates of the direction, where ; For the calculation area in the horizontal direction The field is discretized into Points, Indicates the Points in The coordinates of the direction, where ; For the calculation area at height The field is discretized into Points, Indicates the Points in The coordinates of the direction, where .

[0062] Next, let the grid points at The value is ,definition , as shown in formula (6):

[0063] (6)

[0064] in, Representation plane middle Where l represents the columns of the matrix and m represents the rows of the matrix.

[0065] According to the CMM-3DPE form derived from formula (5), at the midpoint Place The first derivative of the direction is written as:

[0066] (7)

[0067] in, is the grid point at value.

[0068] First, for The first derivative of the direction is expressed as forward difference and averaging as follows:

[0069] (8)

[0070] in, is the grid point at value, is the grid point at value, is the grid point at value, is a grid point at value.

[0071] Secondly, and The second-order derivative of the direction is expressed in the form of central difference and averaging:

[0072] (9)

[0073] (10)

[0074] in, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value.

[0075] Then, about and The central difference of the directional second-order mixed partial derivative is expressed as:

[0076] (11)

[0077] in, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value.

[0078] Finally, substitute formula (7) to formula (11) into formula (5) to obtain:

[0079] (12)

[0080] Simplifying formula (12), we get:

[0081] (13)

[0082] in:

[0083] (14)

[0084] (15)

[0085] (16)

[0086] (17)

[0087] (18)

[0088] (19)

[0089] By replacing formula (13) with a matrix form, we obtain the discretized matrix equation shown in formula (20):

[0090] (20)

[0091] in, Representation plane middle The matrix, 、 、 、 、 、 Represents the coefficient matrix; matrix is the field matrix The result of removing the last column and inserting a column of zero vectors to the left of the first column; the matrix is the field matrix The result of removing the first column and inserting a column of zero vectors to the right of the last column; the matrix is the field matrix The result of removing the last column and inserting a column of zero vectors to the left of the first column; the matrix is the field matrix Result of dropping the first column and inserting a column of zero vectors to the right of the last column.

[0092] The coefficient matrices are:

[0093] (21a)

[0094] (21b)

[0095] (22a)

[0096] (22b)

[0097] (23a)

[0098] (23b)

[0099] matrix Expressed as:

[0100] (twenty four)

[0101] matrix Expressed as:

[0102] (25)

[0103] It is worth noting that when the method of the present invention is used to solve 3DPE, Absorption boundaries need to be set in the a and z directions to simulate the wave propagation to infinity and eliminate the reflection caused by artificial truncation. The absorbing boundary settings are as follows:

[0104] (26)

[0105] (27)

[0106] in, and is the window function, Indicates the coordinate system along The maximum distance propagated in the positive direction of the axis, Indicates the maximum distance propagated along the positive direction of the z-axis in the coordinate system.

[0107] Then 3DPE In-plane absorbing boundary attenuation function for:

[0108] (28)

[0109] Starting from the initial distance, in each matrix calculation, the field distribution calculated by 3DPE needs to be multiplied by .

[0110] Step 5. Incorporate the boundary impedance conditions of the 3DPE model based on the one-side approximation method.

[0111] The impedance boundary condition in the transformed coordinate system is:

[0112] (29)

[0113] in, express The field value at represents the ground surface impedance, .

[0114] On the boundary surface The finite difference expression of the second-order partial derivative of Equation (29) is obtained by using the one-sided approximation, which is approximately expressed as:

[0115] (30)

[0116] in, express The field value at .

[0117] Using central difference approximation at and use the boundary conditions to calculate the boundary at Partial derivatives, we get:

[0118] (31)

[0119] in, express The field value at .

[0120] Substituting formula (28) and formula (31) into the 3D-CMM-PE of formula (28), we obtain:

[0121] (32)

[0122] in, is a grid point at value, Grid Points at value, Grid Points at value, Grid Points at value; is a grid point at value, Grid Points at value, Grid Points at value, Grid Points at value.

[0123] 、 、 、 、 They are:

[0124] (33)

[0125] (34)

[0126] (35)

[0127] (36)

[0128] (37)

[0129] Combining formula (31) into the corresponding coefficient matrix in matrix equation (20) yields:

[0130] (38a)

[0131] (38b)

[0132] (39a)

[0133] (39b)

[0134] because and The form and and The form is the same, so the coefficient matrix constant.

[0135] In step 5 of the method of the present invention, the second-order derivative on the boundary surface is expressed by a unilateral approximation, and the second-order derivative is expressed by an impedance boundary condition. Then, the PE equation, i.e., the CMM-PE equation, which introduces a shift transformation model, is substituted into the equation to obtain an equation that satisfies both the impedance boundary condition and the CMM-PE. This simplifies the amount of calculation while ensuring the calculation accuracy of the 3DPE model.

[0136] Step 6. Based on the discretized matrix equation and the boundary impedance conditions of the 3DPE model, the field matrix is ​​iteratively solved to obtain the field distribution results in the entire calculation area.

[0137] First, let the result on the right side of formula (20) be the matrix ,Right now:

[0138] (40)

[0139] Next, let the matrix for:

[0140] (41)

[0141] Then formula (20) is expressed as:

[0142] (42)

[0143] Let the initial solution , and name the intermediate solution after s iterations as .

[0144] Let the matrix ,in for The intermediate solution after iterations, the matrix is the field matrix The result of removing the last column and inserting a column of zero vectors to the left of the first column; the matrix is the field matrix Result of dropping the first column and inserting a column of zero vectors to the right of the last column.

[0145] The coefficient matrix corresponding to each column field of the field matrix and Are different, multiply each column field by a different coefficient matrix and Average, and solve the field matrix column by column , so the Column This is done by replacing the Shift to the right to calculate.

[0146] (43)

[0147] in, is the coefficient matrix corresponding to the field value of the lth column of the field matrix, Represented as a matrix The first column of is a matrix The first column of is a matrix The first column of .

[0148] middle field is done by The combination obtained, Expressed as:

[0149] (44)

[0150] If the convergence conditions are met, then:

[0151] (45)

[0152] in, is the convergence decision constant, which can generally be taken .

[0153] Solve the intermediate solution by formula (44) In this case The iterative process is to determine the current Does it satisfy formula (45)? If not, let , and return to formula (43) and continue to refine the solution , until formula (45) is satisfied, the calculation is completed, and the field distribution results in the entire calculation area are obtained. .

[0154] Step 7. Using the field results in the entire calculation area, calculate the magnetic field strength of the low-frequency ground wave propagating along the path of the undulating terrain. And the secondary delay generated .

[0155] The calculated value in step 6 Substitute into formula (2) to find all magnetic fields in the calculation area :

[0156] .

[0157] The magnetic field at each location in the calculation area Extract the real and imaginary parts and calculate the phase :

[0158] (46)

[0159] in, represents the operation of extracting the real part, Represents the operation of extracting the imaginary part. It represents the accumulated phase delay of electromagnetic waves in the propagation path.

[0160] When the radio wave propagation path is set to an undulating terrain path, the relative dielectric constant and conductivity are and , the phase of the low-frequency ground wave propagating on the undulating terrain path is obtained by formula (46): .

[0161] When the wave propagation path is set to a good conductor path, the relative dielectric constant and conductivity are and , the phase of the low-frequency ground wave propagating on the good conductor path is obtained by formula (46): .

[0162] Calculate the secondary delay of low-frequency ground waves propagating over undulating terrain The formula is:

[0163] (47).

[0164] In addition, in order to verify the effectiveness of the method proposed in the present invention, the following specific experiments are given:

[0165] Experiment 1 compares the magnetic field intensity and quadratic time delay prediction results of the 3D-CMM-PE method based on the unilateral approximate boundary merging scheme, namely the method of the present invention, the three-dimensional finite-difference time-domain method (3D Finite-Difference Time-Domain, 3D-FDTD) method, and the two-dimensional parabolic equation method based on the shift transformation model (2D-Conform Map Model Parabolic Equation, 2D-CMM-PE).

[0166] In the 3D-CMM-PE calculation model, the calculation area is { , , }, the grid size is defined as m. In the 3D-FDTD calculation model, the calculation area is { , , }, the grid size is defined as m, time step The ground electrical parameters are , S / m, the mountain range is a Gaussian structure, and its height function in the rectangular coordinate system The form is as follows: .

[0167] Among them, H is the height of Gauss Mountain, which is km, and Gauss Mountain and The center of the direction, the value is km, km, and To control the mountain and Parameters of the direction width, the width control parameters are km.

[0168] Figure 3 、 Figure 4 and Figure 7The present invention is based on the 3D-CMM-PE method of the unilateral approximate boundary merging scheme, the 2D parabolic equation 2D-CMM-PE method based on the shift transformation model and the 3D finite difference time domain 3D-FDTD method in Magnetic field strength at the surface on a plane and quadratic delay The prediction results. Figure 3 、 Figure 4 and Figure 7 It can be seen that the field strength and quadratic time delay calculated by the three methods are basically the same. However, the 2D-CMM-PE method cannot predict the influence of lateral terrain on the propagation of low-frequency ground waves, resulting in a significant difference between its results behind the mountain and those of the other two three-dimensional methods.

[0169] Experiment 2 shows the prediction results of the magnetic field strength and secondary time delay of two mountain ranges on the propagation of low-frequency ground waves.

[0170] In this experiment, the calculation area is set the same as in Experiment 1, and the low-frequency ground wave propagation characteristics are predicted by the method proposed in this invention. Terrain height function in the propagation path in the rectangular coordinate system The form is:

[0171] .

[0172] Among them, the height of the two peaks and Both are 1.5km, mountain width control parameters km, km. Fixed the first peak at and The center positions of the directions are km, km; the second peak is The center position of the direction is also km.

[0173] Figure 5 and Figure 6 For two Gaussian peaks, when the position of the first mountain is fixed, the second mountain is The center positions of the directions are 8km on the low-frequency ground wave propagation characteristics. Figure 5 and Figure 6 Can see the magnetic field strength and quadratic delay There is a significant superposition between the two mountains, which is caused by the reflection effect of low-frequency ground wave transverse waves between the two mountains.

[0174] The method of the present invention overcomes the problem that the 2DPE method cannot consider the influence of lateral irregular terrain on low-frequency ground waves, and overcomes the shortcomings of the traditional 3DPE method in solving and predicting the accuracy of low-frequency ground wave propagation characteristics in undulating terrain. The FD algorithm is used to solve 3D-CMM-PE, which not only reduces the influence of dispersion error, but also can predict the lateral diffraction effect caused by the interaction between irregular terrain and low-frequency ground waves, thereby improving the prediction accuracy of low-frequency ground wave propagation characteristics in undulating terrain.

[0175] Example 2

[0176] This embodiment 2 describes a computer device, which includes a memory and one or more processors.

[0177] The memory stores executable codes, which, when executed by the processor, are used to implement the steps of the method for predicting low-frequency ground wave propagation characteristics of a path over undulating terrain based on the FD algorithm in the first embodiment.

[0178] In this embodiment, the computer device is any device or apparatus with data processing capability, which will not be described in detail here.

[0179] Of course, the above description is only a preferred embodiment of the present invention, and the present invention is not limited to the above-mentioned embodiments. It should be noted that all equivalent substitutions and obvious deformation forms made by any technician familiar with this field under the guidance of this specification fall within the substantive scope of this specification and should be protected by the present invention.

Claims

1. A method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm, characterized by: The steps include: Step 1. Input the model file. The content of the input model file includes the grid parameters of the calculation area, the electrical parameters of the radio wave propagation path, and the excitation source parameters. Step 2. Calculate the initial field using the flat ground formula and use the initial field as the source field for 3DPE. Step 3. Establish a 3DPE model under irregular terrain; Step 4. Construct the discretized matrix equation for solving the 3DPE model; For the calculation area along the propagation direction Discretize the field into points, the waveguide consists of a continuous The horizontal plane is composed of Indicates the Points in The coordinates of the direction, where ; For the calculation area in the horizontal direction The field is discretized into Points, Indicates the Points in The coordinates of the direction, where ; For the calculation area at height The field is discretized into Points, Indicates the Points in The coordinates of the direction, where ; The discretized matrix equation is: ； in, Representation plane middle The matrix, Representation plane middle The matrix, is the target field to be solved; 、 、 、 、 、 Represents the coefficient matrix; matrix is the field matrix The result of removing the last column and inserting a column of zero vectors to the left of the first column; the matrix is the field matrix The result of removing the first column and inserting a column of zero vectors to the right of the last column; the matrix is the field matrix The result of removing the last column and inserting a column of zero vectors to the left of the first column; the matrix is the field matrix The result of removing the first column and inserting a column of zero vectors to the right of the last column; Step 5. Merge the boundary impedance conditions of the 3DPE model based on the unilateral approximation method; The boundary impedance condition is: ； in, is a grid point at value, Grid Points at value, Grid Points at value, Grid Points at value; is a grid point at value, Grid Points at value, Grid Points at value, Grid Points at value; 、 、 、 、 They are: (33) (34) (35) (36) (37) in, 、 、 is the spatial step length after coordinate transformation; is the terrain function, is the ground surface impedance, is the wave number of vacuum, is an imaginary unit; Expressed as: ； Step 6. Iteratively solve the field matrix based on the discretized matrix equation and the boundary impedance conditions of the 3DPE model to obtain the field distribution results in the entire calculation area; The discretized matrix equation is simplified to: ； in, , ; Let the initial solution , and name the intermediate solution after s iterations as ; Let the matrix ,in for The intermediate solution after iterations, the matrix is the field matrix The result of removing the last column and inserting a column of zero vectors to the left of the first column; the matrix is the field matrix The result of removing the first column and inserting a column of zero vectors to the right of the last column; ； in, is the coefficient matrix corresponding to the field value of the lth column of the field matrix, Represented as a matrix The first column of is a matrix The first column of is a matrix Column l of middle field is done by The combination obtained, Expressed as: ； The convergence conditions are: ； in, is the convergence decision constant; Step 7. Using the field results in the entire calculation area, calculate the magnetic field strength of the low-frequency ground wave propagating along the path of the undulating terrain. And the secondary delay generated .

2. The method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm according to claim 1, characterized in that: The step 1 is specifically as follows: In the three-dimensional rectangular coordinate system In the calculation area, the size is defined as , where N x is the number of grids in the x-direction of the horizontal coordinate, N y is the number of grids in the y-direction of the longitudinal coordinate, N z is the number of grids in the height coordinate z direction; 、 、 is the meshing step size; The electrical parameters of the radio wave propagation path include the relative dielectric constant of free space and conductivity , relative dielectric constant in mountainous areas and conductivity ; The excitation source is a vertical electric dipole with a current of , the charge distance is , the height from the ground is .

3. The method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm according to claim 2, characterized in that: The step 2 is specifically as follows: Cartesian coordinate system Next, let the time harmonic factor be ,in is the angular velocity, For time; Use the flat ground formula to solve the initial field of 3DPE and the ground wave magnetic field component Expressed as: (1) in, represents the great circle distance; is the wave number on the ground; From source point to observation point The straight-line distance at is the straight-line distance from the source's mirror image to the observation point; is the Fresnel integral, is an intermediate variable; Great Circle Distance Expressed as: ； Fresnel integral Defined as: ； Intermediate variables Expressed as: ； in, is the observation point location, Indicates the vertical height coordinate; For the preset initial distance According to formula (1), the initial field of low-frequency ground wave is obtained for: (2) in, for The ground wave magnetic field component at As the source field of 3DPE.

4. The method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm according to claim 3, characterized in that: The step 3 is specifically as follows: The three-dimensional parabolic equation on a flat surface is expressed as shown in formula (3): (3) in, represents the transverse Laplacian operator, ; is the refractive index of the medium, , is the dielectric constant of the medium, Represents relative magnetic permeability; The coordinate transformation model CMM is used to transform the original undulating terrain into a flat ground, and the height and distance variables are set as follows: (4) in, is the terrain function; after coordinate transformation, the spatial step size satisfies , , ; Substituting formula (4) into formula (3) to simplify the equation, we can obtain the 3DPE model under irregular terrain as shown in formula (5): (5) Among them, the terrain function Simplified to .

5. The method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm according to claim 4, characterized in that: The step 4 is specifically as follows: Let the grid point at The value is ,definition , as shown in formula (6): (6) According to formula (5), at the midpoint Place The first derivative of the direction is written as: (7) for The first-order derivative of the direction is expressed in the form of forward difference and averaging: (8) in, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value; right and The second-order derivative of the direction is expressed in the form of central difference and averaging: (9) (10) in, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value; about and The central difference of the directional second-order mixed partial derivative is expressed as: (11) in, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value, is a grid point at value; Substituting formula (7) to formula (11) into formula (5), we obtain: (12) Simplifying formula (12), we get: (13) in: (14) (15) (16) (17) (18) (19) Substituting formula (13) into matrix form, we obtain the discretized matrix equation shown in formula (20): (20) The coefficient matrices are: (21a) (21b) (22a) (22b) (23a) (23b) matrix Expressed as: (24) matrix Expressed as: (25)。 6. The method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm according to claim 5, characterized in that: In step 4, And set absorbing boundaries in the z direction: (26) (27) in, and is the window function, Indicates the coordinate system along The maximum distance propagated in the positive direction of the axis, Indicates the maximum distance propagated along the positive direction of the z-axis in the coordinate system; Then 3DPE In-plane absorbing boundary attenuation function for: (28)。 7. The method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm according to claim 6, characterized in that: The step 5 is specifically as follows: The impedance boundary condition in the transformed coordinate system is: (29) in, express Field value at ; On the boundary surface Based on the one-sided approximation method, the finite difference expression of the second-order partial derivative of formula (29) is approximately expressed as: (30) in, express Field value at ; Using central difference approximation at and use the boundary conditions to calculate the boundary at Partial derivatives, we get: (31) in, express Field value at ; Substituting formula (28) and formula (31) into formula (5), we obtain: (32) Combining formula (31) into the corresponding coefficient matrix in matrix equation (20) yields: (38a) (38b) (39a) (39b) because and The form and and The form is the same, so the coefficient matrix constant.

8. The method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm according to claim 7, characterized in that: The step 6 is specifically as follows: Let the result on the right side of formula (20) be the matrix ,Right now: (40) Let the matrix for: (41) Then formula (20) is expressed as: (42) Let the initial solution , and name the intermediate solution after s iterations as ; Let the matrix ; (43) middle field is done by The combination obtained, Expressed as: (44) If the convergence conditions are met, then: (45) Solve it by formula (44) , ; Determine the current Does it satisfy formula (45)? If not, let , and return to formula (43) and continue to calculate , until formula (45) is satisfied, and the field distribution results in the entire calculation area are obtained. .

9. The method for predicting low-frequency ground wave propagation characteristics over undulating terrain based on the FD algorithm according to claim 8, characterized in that: The step 7 is specifically as follows: The result calculated in step 6 Substitute into formula (2) to find all magnetic fields in the calculation area : ； The magnetic field at each location in the calculation area Extract the real and imaginary parts and calculate the phase : (46) in, represents the operation of extracting the real part, Indicates the operation of extracting the imaginary part; Represents the cumulative phase delay of the electromagnetic wave in the propagation path; When the radio wave propagation path is set to an undulating terrain path, the relative dielectric constant and conductivity are and , the phase of the low-frequency ground wave propagating on the undulating terrain path is obtained by formula (46): ; When the wave propagation path is set to a good conductor path, the relative dielectric constant and conductivity are and , the phase of the low-frequency ground wave propagating on the good conductor path is obtained by formula (46): ; Calculate the secondary delay of low-frequency ground waves propagating over undulating terrain The formula is: (47)。 10. 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, the steps of the method for predicting low-frequency ground wave propagation characteristics of a path over undulating terrain based on the FD algorithm are implemented as described in any one of claims 1 to 9.

Images