Modeling method for simulating propagation of electric waves in irregular terrain

Abstract

The invention discloses a modeling method for simulating propagation of electric waves in irregular terrains, and relates to the field of radio. The method comprises the following steps: constructing a terrain-related first parabolic equation about radio wave propagation; introducing an auxiliary function, and determining a second parabolic equation; writing the second parabolic equation into a solution in an exponential function form, and performing rational approximation on an exponential term containing a topographic factor and an exponential term containing a pseudo differential factor; solving an intermediate variable and a field function subjected to terrain translation based on a finite difference method; and restoring the field function in combination with the auxiliary function and the terrain profile. Compared with the prior art, the method is high in calculation precision and especially suitable for predicting radio wave propagation loss on irregular terrains.

Description

Technical Field

[0001] The present invention relates to the field of radio technology, and more particularly to a modeling method for simulating radio wave propagation on irregular terrain. Background Art

[0002] The problem of radio wave propagation over irregular terrain has always been a difficult issue in the field of radio wave propagation. The efficient and accurate solution to this problem has attracted widespread attention from scholars and industry. Parabolic equation methods are an important approach to solving this problem. Currently, the split-step Fourier transform method and the finite difference method can be used to solve radio wave propagation over irregular terrain. However, neither method is currently effective in addressing this problem. The former can only solve symmetric parabolic equations and has low accuracy, while the latter cannot handle high-order parabolic equations and has low computational accuracy. Summary of the Invention

[0003] In order to overcome the low precision defect of the above-mentioned prior art, the present invention provides a modeling method for simulating the propagation of radio waves on irregular terrain.

[0004] In order to solve the above technical problems, the technical solutions of the present invention are as follows:

[0005] In a first aspect, a modeling method for simulating radio wave propagation over irregular terrain includes:

[0006] Based on the radio wave propagation characteristics, environmental characteristics, and the topographic profile of the irregular terrain, the first parabolic equation of radio wave propagation related to terrain is constructed;

[0007] introducing an auxiliary function into the first parabolic equation to determine a second parabolic equation suitable for solving by a finite difference method;

[0008] The second parabolic equation is expressed as a solution in the form of an exponential function, and a rational approximation is made to the exponential term containing the terrain factor and the exponential term containing the pseudo-differential factor;

[0009] Based on the finite difference method, the intermediate variables and the field functions after terrain translation are solved;

[0010] The auxiliary function and the terrain profile are combined to restore the field function.

[0011] In a second aspect, a computer program product comprises a computer program or computer executable instructions, wherein when the computer program or computer executable instructions are executed by a processor, the method of the first aspect is implemented.

[0012] Compared with the prior art, the beneficial effects of the technical solution of the present invention are:

[0013] The present invention discloses a modeling method for simulating radio wave propagation over irregular terrain. This method solves the problem of radio wave propagation under irregular terrain conditions by using high-order parabolic equations that carry terrain factors. The method includes constructing a first parabolic equation for radio wave propagation, introducing an auxiliary function to construct a second parabolic equation, writing the second parabolic equation as an exponential function and approximating the solution form to obtain a high-order parabolic equation suitable for irregular terrain. The method then uses the finite difference method and introduces intermediate variables to solve the field function after terrain translation. Finally, the auxiliary function and terrain profile are combined to restore the field function. Compared to existing technologies, the present invention provides higher accuracy and is particularly suitable for predicting radio wave propagation loss over irregular terrain. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 Schematic diagram of the modeling method described in Example 1 of the present application;

[0015] Figure 2 Schematic diagram of the linearized irregular terrain profile and spatial grid division in Example 2 of the present application;

[0016] Figure 3 This is a pseudo-color image of the field intensity of the initial field simulated in Example 2 of the present application;

[0017] Figure 4 This is a comparison diagram of the field calculated in Example 2 of the present application and the geometric diffraction theory. DETAILED DESCRIPTION

[0018] The terms "first", "second" etc. in the specification and claims of the present application and the above-mentioned drawings are used to distinguish similar objects, and are not necessarily used to describe a specific order or sequential order. It should be understood that the terms used in this way can be interchangeable in appropriate circumstances, and this is merely a way of distinguishing the objects of the same attribute when describing them in the embodiments of the present application. In addition, the terms "including" and "having" and any of their variations are intended to cover non-exclusive inclusions, so that the process, method, system, product or equipment comprising a series of units need not be limited to those units, but may include other units that are not clearly listed or inherent to these processes, methods, products or equipment. The term "determine" widely covers various actions, may include obtaining, calculating, computing, processing, deriving, investigating, searching (for example, searching in a table, a database or other data structure), ascertaining and similar actions, may also include receiving (for example, receiving information), accessing (for example, accessing data in a memory) and similar actions, may also include generating, creating, establishing and similar actions, and parsing, selecting, selecting and similar actions etc. The relevant definitions of other terms will be provided in the following description.

[0019] It should be noted that when an element is considered to be "connected" to another element, it can be directly connected to the other element or connected to the other element through an intervening element. In addition, the "connection" in the following embodiments should be understood as "electrical connection", "communication connection", etc., if there is transmission of electrical signals or data between the connected objects.

[0020] The accompanying drawings are for illustrative purposes only and are not to be construed as limiting this patent;

[0021] In order to better illustrate this embodiment, some parts in the drawings may be omitted, enlarged, or reduced, and do not represent the actual product size;

[0022] It is understandable to those skilled in the art that some well-known structures and descriptions thereof may be omitted in the drawings.

[0023] The technical solution of the present invention is further described below with reference to the accompanying drawings and embodiments.

[0024] Example 1

[0025] This embodiment proposes a modeling method for simulating the propagation of radio waves in irregular terrain. Figure 1 ,include:

[0026] S1. Constructing the first parabolic equation related to the terrain for radio wave propagation based on radio wave propagation characteristics, environmental characteristics, and the terrain profile of the irregular terrain;

[0027] S2. Introducing an auxiliary function into the first parabolic equation to determine a second parabolic equation suitable for solving by the finite difference method;

[0028] S3. Rewrite the second parabolic equation as a solution in the form of an exponential function, and make rational approximations to the exponential terms containing terrain factors and the exponential terms containing pseudo-differential factors;

[0029] S4. Based on the finite difference method, solve the intermediate variables and the field function after terrain translation;

[0030] S5. Restore the field function by combining the auxiliary function and the terrain profile.

[0031] In some preferred embodiments, the radio wave propagation characteristics include radio wave propagation direction, simulation frequency and wave number, and the environmental characteristics include components of the magnetic field or electric field and atmospheric refractive index;

[0032] The constructing of the first parabolic equation related to the terrain for radio wave propagation is specifically constructing a parabolic equation with a terrain slope, including:

[0033] Assume that the horizontal axis is used to establish a vw rectangular coordinate system. The irregular terrain function in the vw coordinate system is T(v). Next, the xz rectangular coordinate system is established with the height of the slope boundary of each irregular surface segment as the horizontal boundary z = 0. The conversion relationship between the vw and xz rectangular coordinate systems is:

[0034]

[0035] In the vw rectangular coordinate system, the wave equation is:

[0036]

[0037] Where w represents the height, v represents the propagation direction, φ represents the horizontally or vertically polarized magnetic or electric field component, n represents the atmospheric refractive index, and k represents the wave number. Assuming that the slope S of each terrain segment is constant, vw is converted to the xz coordinate system based on (2-1) as follows:

[0038]

[0039] ψ represents φ after coordinate transformation.

[0040] Substituting (2-3) and (2-4) into (2-2) and decomposing them, we obtain the following first parabolic equation:

[0041]

[0042] Where S represents the slope of each section of terrain, which is a constant.

[0043] In some optional embodiments, step S2 includes: introducing an auxiliary function to obtain a second parabolic equation suitable for solving by a finite difference method, wherein the expression of the auxiliary function is as follows:

[0044] ψ(x,z)=e ikx u(x,z) (2-6)

[0045] Substituting (2-6) into (2-5) we get

[0046]

[0047] Equation (2-7) is the second parabolic equation suitable for solving by the finite difference method; m is the atmospheric refractive index after n is transformed into a coordinate system.

[0048] It should be noted that the introduction of the auxiliary function can significantly reduce the truncation error of the difference format when calculating small-angle plane wave components.

[0049] Furthermore, in step S3, the expression of the solution of the exponential function form of the second parabolic equation is:

[0050]

[0051] in,

[0052]

[0053] Δx is the step size of the horizontal grid. It should be understood that Z1 and Z2 are used to simplify the expression of formula (2-8).

[0054] Furthermore, in step S3, a rational approximation is performed on the exponential term containing the pseudo differential factor, and its expression is:

[0055]

[0056] Furthermore, in step S3, a rational approximation is made to the exponential term containing the terrain factor, that is, the exponential term with S is approximated as follows:

[0057]

[0058] Furthermore, the finite difference method is used to solve the intermediate variables and the field function after terrain translation, including:

[0059] Use u p represents the numerical calculation result of u at the position pΔx, and similarly ψ p represents the numerical calculation result of ψ at the position pΔx; u p (z j ) represents the numerical calculation result of u(pΔx,jΔz), and similarly ψ p (z j ) is also the numerical calculation result of ψ(pΔx,jΔz). (Written in full as ), represents an intermediate variable, and the subscript s represents an ordinal number (e.g., The subscript s in the The ordinal number of the grid), Δx represents the horizontal grid step, and Δz represents the vertical grid step. The solution process is as follows:

[0060]

[0061] In particular, when g = 1, equation (2-13) becomes

[0062]

[0063] When g = 2, formula (2-13) becomes

[0064]

[0065] When l=1, formula (2-13) becomes

[0066]

[0067] When l=2, formula (2-13) becomes

[0068]

[0069] The intermediate process of formula (2-13)-(2-17) can be written as follows:

[0070]

[0071] Then use difference instead of differential, which is implemented as follows:

[0072]

[0073] Where E represents the translation operation, that is, E j u(0)=u(z j ) (It should be noted that the absence of superscripts in equations (2-23) and (2-24) indicates that u p-1 、u p E has the same effect on v and w. χ、 is the difference format parameter, which can be determined by combining series expansion and undetermined coefficients.

[0074] For v and w, the same parameters can be used for difference instead of differentiation. It should be noted that the order of solving v and w can be swapped, that is, the positions of d and b, c and a, w and v, and Z1 and Z2 can be swapped.

[0075] Substituting equations (2-23) and (2-24) into equations (2-18)-(2-22) yields a system of linear equations. However, the number of unknowns in the resulting system is greater than the number of equations. Therefore, it is necessary to introduce impedance boundary conditions and expand the number of equations to achieve the same number of equations as the number of unknowns. The impedance boundary condition is expressed as follows:

[0076]

[0077] Here, η is related to the material of the boundary.

[0078] Since (2-25) and (2-27) actually only add one equation group, they may not necessarily meet the requirements for solving the equation group. Therefore, high-order impedance boundaries are introduced, and their expressions are as follows:

[0079]

[0080] For equations (2-28)-(2-30), the discretization method is based on the single-sided difference format. The specific expressions are as follows:

[0081]

[0082] Combining (2-13)-(2-33), we can p-1 Solution and the terrain-translated field function u p+1 .

[0083] Furthermore, in step S6, the field function is restored by combining the auxiliary function, and its expression is:

[0084] ψ p (z j )=e ikpΔx u p (z j ) (2-34)

[0085] Furthermore, in step S6, the field function is restored in combination with the terrain profile, and its expression is:

[0086] φ p (z j +T p )=ψ p (z j ) (2-35)

[0087] Here φ p represents the restored field function, T p Represents the height of the piecewise linear terrain at position pΔx.

[0088] Furthermore, to prevent the reflection of the upper boundary, a window function is added to absorb the boundary, including:

[0089] At height j a Δz is followed by a window function, [j a Δz,j max The field function in the interval Δz] gradually becomes 0 to avoid reflection at the boundary.

[0090] In some specific implementations, a window function such as a Hanning window may be used.

[0091] Furthermore, to prevent reflection of the upper boundary, a PML may be introduced in step S3, including:

[0092] The complex coordinate axis is used instead of the original real coordinate axis to realize the attenuation of the field function. The expression of the complex coordinate axis is as follows:

[0093]

[0094] Where, represents complex coordinates;

[0095] In complex coordinates, it needs to be written in the following form:

[0096]

[0097] In complex coordinates, it needs to be written in the following form:

[0098]

[0099] Among them, the value of γ(z) changes gradually (note that this is the step S4, so the field function has not been translated yet), starting from 0 and gradually increasing to avoid reflection introduced by coordinate discontinuity.

[0100] This embodiment also provides a computer-readable storage medium, on which is stored at least one instruction, at least one program, code set or instruction set, and the at least one instruction, at least one program, code set or instruction set is loaded and executed by a processor, so that the processor performs some or all steps of the method provided in the aforementioned embodiment.

[0101] It is understood that the storage medium may be transient or non-transient. Exemplarily, the storage medium includes, but is not limited to, a USB flash drive, a mobile hard drive, a read-only memory (ROM), a random access memory (RAM), a magnetic disk, or an optical disk, among other media capable of storing program code.

[0102] Exemplarily, the processor may be a central processing unit (CPU), a microprocessor (MPU), a digital signal processor (DSP), or a field programmable gate array (FPGA).

[0103] Exemplarily, the read-only memory includes but is not limited to MASK ROM, PROM, EPROM, EEPROM, Flash, etc.

[0104] Exemplarily, the random access memory includes but is not limited to DRAM, SRAM, SDRAM, DDR SDRAM, etc.

[0105] In some examples, a computer program product is provided, which can be implemented in hardware, software, or a combination thereof. As a non-limiting example, the computer program product can be embodied as the storage medium, or as a software product, such as an SDK (Software Development Kit).

[0106] As a non-limiting example, a computer program product is provided, comprising a computer program or computer-executable instructions stored in a computer-readable storage medium. A processor of an electronic device reads the computer program or computer-executable instructions from the computer-readable storage medium and executes the computer-executable instructions, causing the electronic device to perform some or all of the steps of the method described in the embodiments of the present application.

[0107] In some examples, a computer program is provided, comprising a computer-readable code. When the computer-readable code is run in a computer device, a processor in the computer device executes the code to implement part or all of the steps in the method.

[0108] This embodiment also proposes an electronic device, including a memory and a processor, wherein the memory stores at least one instruction, at least one program, code set or instruction set, and when the processor executes the at least one instruction, at least one program, code set or instruction set, it implements part or all of the steps of the method described in the above embodiment.

[0109] In some examples, a hardware entity of the electronic device is provided, including: a processor, a memory and a communication interface; wherein the processor generally controls the overall operation of the electronic device; the communication interface is used to enable the electronic device to communicate with other terminals or servers through a network; the memory is configured to store instructions and applications executable by the processor, and can also cache data to be processed or processed by the processor and various modules in the electronic device (including but not limited to image data, audio data, voice communication data and video communication data), which can be implemented by flash memory (FLASH), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM) or random access memory (RAM).

[0110] A processor may include one or more processing elements. Thus, a processor may include one or more integrated circuits (ICs) configured to perform the functions of the processor. Furthermore, each integrated circuit may include circuits (e.g., a first circuit, a second circuit, and other circuits) configured to perform the functions of the processor.

[0111] Furthermore, data may be transmitted between the processor, the communication interface and the memory via a bus, which may include any number of interconnected buses and bridges, connecting various circuits of one or more processors and memories.

[0112] Example 2

[0113] To facilitate those skilled in the art to implement this application, this embodiment provides a more specific simulation implementation process based on Example 1, including:

[0114] S1: Construct the first parabolic equation related to the terrain for radio wave propagation based on the direction of radio wave propagation, propagation direction, components of the magnetic or electric field, atmospheric refractive index, simulation frequency, wave number, and profile of irregular terrain.

[0115] S2: Introduce an auxiliary function into the first parabolic equation to obtain the second parabolic equation suitable for solving by the finite difference method.

[0116] S3: Write the solution of the second parabolic equation in the form of an exponential function.

[0117] S4: Make rational approximations to the exponential terms containing terrain factors and the exponential terms containing pseudo-differential factors.

[0118] S5: Use the finite difference method, combined with approximations and impedance boundary conditions, to solve for the intermediate variables and the translated field functions.

[0119] S6: Combine auxiliary functions and terrain profiles to restore field functions.

[0120] S7: To prevent reflection at the upper boundary, add a window function to absorb the boundary or introduce PML in step S5.

[0121] In this embodiment, step S1 includes: constructing a first parabolic equation with terrain slope, the process of which is as follows:

[0122] See Figure 2 Assume that the horizontal axis is used to establish the vw rectangular coordinate system. The irregular terrain function in the vw coordinate system is T(v). The xz rectangular coordinate system is established with the height of the slope boundary of each irregular surface as the horizontal boundary z = 0. The conversion relationship between the vw and xz rectangular coordinate systems is:

[0123]

[0124] In the vw rectangular coordinate system, the wave equation is:

[0125]

[0126] Where w represents the height, v represents the propagation direction, ψ represents the horizontally or vertically polarized magnetic or electric field component, n represents the atmospheric refractive index, and k represents the wave number. Assuming that the slope S of each terrain segment is constant, vw is converted to the xz coordinate system based on (3-1) as follows:

[0127]

[0128] ψ represents φ after coordinate transformation.

[0129] Substituting (3-3) and (3-4) into (3-2) and performing decomposition, we obtain the following first parabolic equation:

[0130]

[0131] In this embodiment, step S2 includes: introducing an auxiliary function to obtain a second parabolic equation suitable for solving by the finite difference method. The expression of the auxiliary function is as follows:

[0132] ψ(x,z)=e ikx u(x,z) (3-6)

[0133] In this embodiment, in step S2, the second parabolic equation applicable to the finite difference method for solution has the following derivation process and expression:

[0134] Substituting (3-6) into (3-5) we get

[0135]

[0136] Equation (3-7) is the second parabolic equation suitable for solving by the finite difference method.

[0137] In this embodiment, in step S3, the solution in the form of an exponential function is expressed as follows:

[0138]

[0139] in,

[0140]

[0141] Δx is the step size of the horizontal grid.

[0142] In this embodiment, in step S4, an approximate transformation is performed on the exponential term with the terrain slope S, and its expression is as follows:

[0143]

[0144] In this embodiment, in step S4, a rational approximation is performed on the exponential term containing the pseudo-differential term, and its expression is as follows:

[0145]

[0146] Assume that the simulation frequency is 300 MHz, the speed of light is 3e8 m / s, the horizontal grid step size Δx = 0.25 m, the vertical grid step size Δz = 0.1 m, and the atmospheric refractive index is 1. The parameters of a and b are shown in Table 1:

[0147] Table 1 Numerical values of parameters a and b

[0148] a b 0.254527097037178+0.102770902139200i 0.891170206488095-0.0242267631058184i 0.888333285237185+0.0139499091667711i 0.610500586375455-0.0884391409569107i 0.0236884063317687+0.105731792581913i 0.284488919118403-0.161489102867564i 0.598097595061600+0.0519099433553410i 0.0458018794189215-0.165041847832471i 0 -0.0673152077331475-0.0718387613914577i

[0149] In this embodiment, step S5 includes: using the finite difference method, combined with the approximate terms and the impedance boundary conditions, to solve the intermediate variables and the translated field function, the process is as follows:

[0150] Use u p (z j ) represents the numerical calculation result of u(pΔx,jΔz), ψ p (z j ) represents the numerical calculation result of ψ(pΔx,jΔz), and the intermediate variable is The solution process is as follows:

[0151]

[0152] The intermediate process in (3-13) can be written as follows:

[0153]

[0154] The last step of (3-13) is written as follows:

[0155]

[0156] Then, we use the difference to replace the partial derivatives of equations (3-13) to (3-16) with respect to z. The implementation method is as follows:

[0157]

[0158] The same method is used to discretize v and w. Substituting (3-17) and (3-18) into equations (3-15) to (3-16) yields a set of linear equations.

[0159] In step S5, the parameters of the impedance boundary condition are as follows:

[0160]

[0161] η=α-S (3-21)

[0162]

[0163] Among them, ε r is the relative complex permittivity of the lower boundary and satisfies the following relationship:

[0164] ε r =ε g +i60σλ (3-23)

[0165] Among them, ε g represents the relative permittivity of the lower boundary, σ represents the conductivity of the lower boundary, and λ represents the wavelength of the radio wave.

[0166] In this embodiment, it is assumed that the polarization direction of the radio wave is horizontal polarization, σ=6.17e7, ε g =1.

[0167] In this embodiment, a single-side differential format is used to discretize the impedance boundary, as follows:

[0168]

[0169] v p (z0) and u(z0) can be expressed by v(z1) and u(z1), then (3-24) and (3-26) can be used to supplement the number of equations so that v and u can be solved.

[0170] In this embodiment, in step S6, the field function is restored by combining the auxiliary function, and its expression is as follows:

[0171] ψ p (z j )=e ikpΔx u p (z j ) (3-27)

[0172] In this embodiment, in step S6, the field function is restored in combination with the terrain profile, and its expression is as follows:

[0173] φ p (z j +T p )=ψ p (z j ) (3-28)

[0174] Where, φ p represents the restored field function.

[0175] In this embodiment, step S7 includes: preventing reflection at the upper boundary, using a Hanning window as an absorption window, which is expressed as follows:

[0176]

[0177] In this embodiment, Z0=1000m, Z max =2000m.

[0178] In this embodiment, the initial field used is a Levy source, such as Figure 3 As shown, its expression is as follows:

[0179]

[0180] Where β is the half-power beam angle, which in this case is 1°. θ is the source's tilt angle, which in this case is aimed at the top of the wedge, so θ = 1.7905°. The wedge's parameters are as follows: a height of 300 m, a slope of 30°, and a location of 8 km.

[0181] Figure 4 In this embodiment, Figure 3 From the modeling results of the field strength at an altitude of 300 m, it can be seen that the calculation accuracy of the method described in this embodiment is higher than that of the geometric diffraction theory.

[0182] It can be understood that the options in the above embodiment 1 are also applicable to this embodiment, so they will not be described again here.

[0183] The same or similar reference numerals correspond to the same or similar components;

[0184] The terms used in the drawings to describe positional relationships are for illustrative purposes only and are not to be construed as limiting the present application.

[0185] It should be noted that, unless there is any conflict, the embodiments and features in the embodiments of this application can be combined with each other.

[0186] In different specific implementations, the method or system described in this application can be implemented in software, hardware or a combination thereof. In addition, the order of the steps of the method can be changed, and various elements can be added, reordered, combined, omitted, modified, etc.

[0187] Obviously, the above embodiments of the present application are merely examples for clearly illustrating the present application, and are not intended to limit the implementation methods of the present application, and are not intended to limit the present application. For those skilled in the art, other different forms of changes or modifications can be made based on the above description. Each discrete structural / functional module or unit can be integrated together to form an independent part, or each module can exist alone, or two or more modules can be integrated to form an independent part, and the structure and function of the discrete components can be implemented as a combined structure or component. It is not necessary and impossible to enumerate all the implementation methods here. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present application should be included in the scope of protection of the claims of the present application.

Claims

1. A modeling method for simulating radio wave propagation in irregular terrain, characterized in that: include: Based on the radio wave propagation characteristics, environmental characteristics and the terrain profile of the irregular terrain, the first parabolic equation related to the terrain of radio wave propagation is constructed; introducing an auxiliary function into the first parabolic equation to determine a second parabolic equation suitable for solving by a finite difference method; The second parabolic equation is expressed as a solution in the form of an exponential function, and a rational approximation is made to the exponential term containing the terrain factor and the exponential term containing the pseudo-differential factor; Based on the finite difference method, the intermediate variables and the field functions after terrain translation are solved; The auxiliary function and the terrain profile are combined to restore the field function.

2. A modeling method for simulating radio wave propagation in irregular terrain according to claim 1, characterized in that: The radio wave propagation characteristics include the radio wave propagation direction, simulation frequency and wave number, and the environmental characteristics include the components of the magnetic field or electric field and the atmospheric refractive index; The constructing of a first parabolic equation related to terrain for radio wave propagation includes: With the horizontal axis as the abscissa, a vw rectangular coordinate system is established. The irregular terrain function in the vw coordinate system is T(v). Next, the xz rectangular coordinate system is established with the height of the slope boundary of each irregular terrain segment as the horizontal boundary z = 0. The conversion relationship between the vw and xz angular coordinate systems is: In the vw rectangular coordinate system, the wave equation is: Where w represents the altitude, v represents the propagation direction, φ represents the horizontally or vertically polarized component of the magnetic or electric field, n represents the atmospheric refractive index, and k represents the wave number; Convert the vw coordinate system to the xz coordinate system: ψ represents φ after coordinate transformation; Substituting into the wave equation and decomposing it, we get the first parabolic equation: Where S represents the slope of each section of terrain, which is a constant.

3. A modeling method for simulating radio wave propagation in irregular terrain according to claim 2, characterized in that: The expression of the auxiliary function is: ψ(x,z)=e ikx u(x,z) (1-6) Substituting into the first parabolic equation, the expression of the second parabolic equation is: Where m is the atmospheric refractive index after the coordinate transformation of n.

4. A modeling method for simulating radio wave propagation in irregular terrain according to claim 3, characterized in that: The expression of the solution of the exponential function form of the second parabolic equation is: in, Where Z1 and Z2 are used to simplify the expression of formula (1-8); Making rational approximations for the terms containing pseudo-differential terms, we have: Where Δx represents the grid step in the horizontal direction; The rational approximation is also made for the term containing the terrain slope exponent:

5. The modeling method for simulating radio wave propagation in irregular terrain according to claim 4, characterized in that: The method of solving the intermediate variables and the field function after terrain translation based on the finite difference method includes: Use u p represents the numerical calculation result of u at the position pΔx, and similarly ψ p represents the numerical calculation result of ψ at the position pΔx; u p (z j ) represents the numerical calculation result of u(pΔx,jΔz), and similarly ψ p (z j ) is also the numerical calculation result of ψ(pΔx,jΔz); is an intermediate variable, the subscript s represents the ordinal number; Δx represents the horizontal grid step, and Δz represents the vertical grid step; The solution process is as follows: In particular, when g = 1, equation (1-13) becomes When g = 2, formula (1-13) becomes When l=1, formula (1-13) becomes When l=2, formula (1-13) becomes Among them, the intermediate process of formula (1-13)-(1-17) is written as follows: Then, we use the difference to replace the partial derivative of z in the solution process, and we have: Where E represents the translation operation, that is, E j u(0)=u(z j ); E has the same effect on v and w; β, χ、 is the difference format coefficient; for v and w, the same parameters are used for difference instead of differentiation; Substituting equations (1-23) and (1-24) into equations (1-18)-(1-22) yields a system of linear equations. However, the number of unknowns in the resulting system is greater than the number of equations. Therefore, the impedance boundary condition is introduced to expand the number of equations to achieve the condition where the number of equations is the same as the number of unknowns. The expression for the impedance boundary condition is as follows: Where η represents a parameter related to the material of the boundary; Introducing a higher-order impedance boundary, the expression is as follows: For equations (1-28)-(1-30), the impedance boundary is discretized using a single-sided difference format. The specific expressions are as follows: Combined with (1-13)-(1-33), based on u p-1 Solution and the terrain-translated field function u p+1 .

6. A modeling method for simulating radio wave propagation in irregular terrain according to claim 5, characterized in that: The combining of the auxiliary function and the terrain profile to restore the field function includes: Combined with the auxiliary function, u p Restore to ψ p , whose expression is: The field function is restored by combining the terrain profile, and its expression is: φ p (z j +T p )=ψ p (z j ) (1-35) Where, φ p represents the restored field function, T p Represents the height of the piecewise linear terrain at position pΔx.

7. A modeling method for simulating radio wave propagation in irregular terrain according to claim 6, characterized in that: To prevent reflection of the upper boundary, the method further includes: adding a window function absorbing boundary, or introducing a PML before solving the intermediate variable and the field function after terrain translation.

8. The modeling method for simulating radio wave propagation in irregular terrain according to claim 7, characterized in that: The adding of the window function boundary comprises: At height j a Δz is followed by a window function, [j a Δz,j max The field function in the interval Δz] gradually becomes 0 to avoid reflection at the boundary.

9. The modeling method for simulating radio wave propagation in irregular terrain according to claim 7, characterized in that: The introducing of PML comprises: The complex coordinate axis is used instead of the original real coordinate axis to realize the attenuation of the field function. The expression of the complex coordinate axis is as follows: Where, represents complex coordinates; In complex coordinates, it needs to be written in the following form: In complex coordinates, it needs to be written in the following form: The value of γ(z) changes gradually from 0 to increase gradually to avoid reflection introduced by coordinate discontinuity.

10. A computer program product comprising a computer program or computer executable instructions, characterized in that When the computer program or computer executable instructions are executed by a processor, the method according to any one of claims 1 to 9 is implemented.

Images