CN117992707B - Method for predicting low-frequency electric wave propagation characteristics in complex path based on integral equation - Google Patents
Method for predicting low-frequency electric wave propagation characteristics in complex path based on integral equation Download PDFInfo
- Publication number
- CN117992707B CN117992707B CN202410381119.1A CN202410381119A CN117992707B CN 117992707 B CN117992707 B CN 117992707B CN 202410381119 A CN202410381119 A CN 202410381119A CN 117992707 B CN117992707 B CN 117992707B
- Authority
- CN
- China
- Prior art keywords
- integral
- formula
- equation
- point
- integration
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 110
- 238000004364 calculation method Methods 0.000 claims abstract description 35
- 230000005684 electric field Effects 0.000 claims abstract description 23
- 230000010354 integration Effects 0.000 claims description 47
- 230000005284 excitation Effects 0.000 claims description 7
- 238000006467 substitution reaction Methods 0.000 claims description 5
- 230000035699 permeability Effects 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 2
- 239000004020 conductor Substances 0.000 claims description 2
- 238000012885 constant function Methods 0.000 claims description 2
- 238000005070 sampling Methods 0.000 claims description 2
- 230000001131 transforming effect Effects 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 16
- 238000010586 diagram Methods 0.000 description 5
- 238000004088 simulation Methods 0.000 description 4
- 239000000243 solution Substances 0.000 description 4
- 230000005428 wave function Effects 0.000 description 2
- 102100029469 WD repeat and HMG-box DNA-binding protein 1 Human genes 0.000 description 1
- 101710097421 WD repeat and HMG-box DNA-binding protein 1 Proteins 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000012088 reference solution Substances 0.000 description 1
- HOWHQWFXSLOJEF-MGZLOUMQSA-N systemin Chemical compound NCCCC[C@H](N)C(=O)N[C@@H](CCSC)C(=O)N[C@@H](CCC(N)=O)C(=O)N[C@@H]([C@@H](C)O)C(=O)N[C@@H](CC(O)=O)C(=O)OC(=O)[C@@H]1CCCN1C(=O)[C@H]1N(C(=O)[C@H](CC(O)=O)NC(=O)[C@H](CCCN=C(N)N)NC(=O)[C@H](CCCCN)NC(=O)[C@H](CO)NC(=O)[C@H]2N(CCC2)C(=O)[C@H]2N(CCC2)C(=O)[C@H](CCCCN)NC(=O)[C@H](CO)NC(=O)[C@H](CCC(N)=O)NC(=O)[C@@H](NC(=O)[C@H](C)N)C(C)C)CCC1 HOWHQWFXSLOJEF-MGZLOUMQSA-N 0.000 description 1
- 108010050014 systemin Proteins 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Landscapes
- Complex Calculations (AREA)
Abstract
The invention belongs to the technical field of low-frequency wave propagation, and particularly discloses a method for predicting low-frequency wave propagation characteristics in a complex path based on an integral equation. The method comprises the following steps: inputting a model file; decomposing an integral term in the electric wave attenuation function at the receiving point P into three sections I 1、I2、I3 by utilizing a sectional integral method; solving the integral of I 1 by adopting a Gaussian integral formula; solving the integral of I 2 based on a simpson formula of quadratic interpolation polynomial approximation and a method of cubic interpolation polynomial approximation; solving the integral of I 3 by adopting quadratic interpolation polynomial approximation and a method of transforming integral limits; and combining the three-section integral result to obtain an electric wave attenuation function, and calculating the total electric field and the secondary time delay of electric wave propagation. The method improves the calculation efficiency, ensures the calculation precision, realizes the modeling of the long-distance complex path by utilizing the non-uniform grid, and expands the applicability and flexibility of the integral equation method.
Description
Technical Field
The invention belongs to the technical field of low-frequency wave propagation, and particularly relates to a method for predicting low-frequency wave propagation characteristics in a complex path based on an integral equation.
Background
The low-frequency electric wave is widely used in telecommunication and navigation services by virtue of its advantages such as long wavelength, reduced attenuation, and stable propagation characteristics. In a long-distance complex propagation path, irregular and uneven terrain environments may seriously affect propagation characteristics of low-frequency electric waves. Currently, the integral equation method is one of the important methods for predicting the propagation of low frequency electric waves in such environments. However, most of the existing numerical techniques solve the integral equation based on uniform grids, and for long-distance complex path problems, smaller grid sizes are required to ensure good calculation accuracy, but the adoption of uniform small grids on the whole propagation path leads to a reduction in calculation efficiency. In addition, the general linear approximation method (such as a trapezoidal formula) for solving the integral equation has the defect of lower precision.
Therefore, for the long-distance complex path problem, a numerical method needs to be proposed to solve the integral equation so that the calculation accuracy and the calculation efficiency of the low-frequency electric wave propagation characteristics can be ensured at the same time.
Disclosure of Invention
The invention aims to provide a method for predicting low-frequency wave propagation characteristics in a complex path based on an integral equation, which aims to solve the problems of low calculation efficiency and large occupied computer memory in the traditional method for predicting the low-frequency wave propagation characteristics of the complex path by adopting a fine-grid integral equation method with uniform grid size, and also can ensure the calculation accuracy of the low-frequency wave propagation characteristics.
In order to achieve the above purpose, the invention adopts the following technical scheme:
the method for predicting the low-frequency electric wave propagation characteristics in the complex path based on the integral equation comprises the following steps:
step 1, inputting a model file, wherein the content of the input model file comprises grid parameters of a calculation area, excitation source parameters, electric parameters of an electric wave propagation path and electric parameters of a free space;
step 2, utilizing a sectional integration method to attenuate the electric wave function at the receiving point P The integral term contained in (a) is decomposed into the following three segments: integration involving singular point 0Principal integralIncluding singular pointsAnd unknown integralIntegral of (2)；
Step 3, calculating by adopting a 5-point Legend-Gaussian integration method and combining a flat ground formulaIs a combination of the integration of (2);
Step 4. For the composition containing Of the integral intervalsIntegration is calculated in two cases, namely:
When (when) When the number is even, a simpson formula based on quadratic interpolation polynomial approximation is adopted to solveIs a combination of the integration of (2);
When (when) When odd, forFront partThe quadratic interpolation polynomial approximation-based simpson formula is adopted for solving the integral intervalsThe integration of the interval is solved by adopting a cubic interpolation polynomial approximation method in the last 3 integral intervals;
Step 5, solving by adopting a quadratic interpolation polynomial approximation and transformation integral limit method Is a combination of the integration of (2);
Step 6. By integration 、、Solving the electric wave attenuation function at the receiving point PAnd obtaining the total electric field in the whole calculation area and the secondary time delay generated by the propagation of the low-frequency electric wave on the complex path through the electric wave attenuation function.
The invention has the following advantages:
As described above, the present invention describes a method for predicting low frequency wave propagation characteristics in a complex path based on integral equations, which derives a general numerical solution of the integral equation of a hybrid grid (i.e., using both coarse and fine grids) based on the simpson formula product of quadratic interpolation polynomial approximation, and improves the conventional integral equation method using only fine or coarse grids. Compared with a trapezoidal formula product method by means of linear approximation, the method provided by the invention has higher calculation accuracy in predicting the low-frequency wave propagation characteristics in a complex path. In addition, the method of the invention models the long-distance complex path by utilizing the mixed grid technology, compared with the integral equation method under the traditional uniform grid, the method expands the applicability and flexibility of the integral equation method, and improves the calculation efficiency of predicting the low-frequency wave propagation characteristics in the complex path under the condition of not losing the precision.
Drawings
FIG. 1 is a flow chart of a method for predicting propagation characteristics of low frequency electric waves in a complex path based on an integral equation in the present invention;
FIG. 2 is a schematic diagram of a method for predicting propagation characteristics of low frequency electric waves in a complex path based on an integral equation in the present invention;
FIG. 3 is a schematic diagram showing the comparison of the electric field intensity of the improved integral equation method of the present invention with the original integral equation method (coarse grid and fine grid) and the time domain finite difference method on Gao Sishan paths with heights of 1km and 2km, respectively, in embodiment 1;
FIG. 4 is a schematic diagram showing the comparison of the second time delay of the improved integral equation method of the present invention with the original integral equation method (coarse grid and fine grid) and the time domain finite difference method on Gao Sishan paths with heights of 1km and 2km respectively in embodiment 1;
FIG. 5 is a schematic diagram showing the comparison of the electric field intensity of the improved integral equation method of the present invention with the original integral equation method (coarse grid and fine grid) and the time domain finite difference method on the triangular mountain path with the height of 1km and 2km respectively in embodiment 2;
fig. 6 is a schematic diagram showing the comparison of the second time delay of the improved integral equation method of the present invention with the original integral equation method (coarse grid and fine grid) and the time domain finite difference method on the paths of the triangular mountains with the heights of 1km and 2km respectively in embodiment 2.
Detailed Description
Noun interpretation:
emission point: the position in the two-dimensional coordinate system where the low-frequency radio wave signal is transmitted is defined as the origin O of coordinates in the present invention.
Receiving points: the position in the two-dimensional coordinate system at which the low-frequency radio wave signal is received is defined as point P in the present invention.
Integration point: the position of the integral I, which in the present invention is defined as point Q, needs to be calculated in a two-dimensional coordinate system.
The invention is described in further detail below with reference to the attached drawings and detailed description:
The embodiment of the invention discloses a method for predicting low-frequency wave propagation characteristics in a complex path based on an integral equation. Meanwhile, the method adopts the non-uniform grid to model the long-distance complex path, and compared with the traditional uniform grid scheme, the method can simultaneously consider the calculation precision and the calculation efficiency when predicting the low-frequency wave propagation characteristics in the complex path.
The method of the invention has the general flow: firstly, utilizing a sectional integration method to attenuate the electric wave at a receiving point PThe integral term in (3) is decomposed into three sections I 1、I2、I3; secondly, solving the integral of I 1 by adopting a Gaussian integral formula; then, solving the integral of I 2 by adopting a simpson formula based on quadratic interpolation polynomial approximation and a cubic interpolation polynomial approximation method; then, solving the integral of I 3 by adopting quadratic interpolation polynomial approximation and a method of transforming integral limits; finally, the electric wave attenuation function is obtained by combining the integral result of I 1、I2、I3 And utilizeThe total electric field and the secondary time delay of the wave propagation are calculated.
Based on the above inventive concept, the method for predicting the propagation characteristics of the low-frequency electric wave in the complex path based on the integral equation according to the present invention will be described in detail with reference to the accompanying drawings, and as shown in fig. 1 and 2, the method of the present invention comprises the following steps:
Step 1, inputting a model file, wherein the content of the input model file comprises grid parameters of a calculation area, excitation source parameters, electric parameters of an electric wave propagation path and electric parameters of a free space.
In a two-dimensional rectangular coordinate systemIn defining the length of the calculation region asAnd willDiscrete intoA plurality of grids, each grid having a length ofAnd satisfies the following relationship: (1)
Wherein x represents the lateral distance, z represents the longitudinal height, Representation ofFirst in the directionThe integration length of the individual integration points,，And (2) andInitial value of (1)Subscript ofRepresent the firstAnd integral subintervals.
The mixed grids are adopted for different terrains on the low-frequency wave propagation path, namely grids with different sizes are adopted for different terrains, coarse grids are adopted when the propagation path is a flat ground, and fine grids are adopted when the propagation path is a fluctuation mountain.
The mixed grids are adopted for different terrains on the low-frequency wave propagation path, namely, grids with different sizes are adopted for different terrains, and when the propagation path is a flat ground or a fluctuant mountain, coarse grids or fine grids are adopted respectively.
For example, when the propagation path is a flat ground, a coarse mesh is used, the mesh length of which can be set toWhen the propagation path is a undulating mountain, a fine mesh is used, the mesh length of which can be set to=100m。
The excitation source is a vertical electric dipole, and the current and charge distance are respectively I,The electric parameters of the electric wave propagation path include the relative dielectric constant epsilon r and the electric conductivity sigma of the earth surface; the electrical parameters of free space include the air permittivity epsilon 0 and the air permeability mu 0.
Step 2, utilizing a sectional integration method to attenuate the electric wave function at the receiving point PThe integral term contained in (a) is decomposed into the following three segments: integration involving singular point 0Principal integralIncluding singular pointsAnd unknown integralIntegral of (2)。
An electric wave attenuation function represented by Hufford integral equation at a receiving point PGiven by equation (2):
(2)
Wherein P is the position of the receiving point, and its coordinates are Q is the position of the surface area minute point of the terrain, and the coordinates areS represents the distance coordinate of the Q point in the x direction in a two-dimensional coordinate system (x, z),Represents the terrain height, W (Q) represents the electric wave attenuation function of the point Q,Is the distance from the emission point O to the reception point P, i represents imaginary units, λ is the wavelength of free space; The definition is as follows:
(3)
Wherein, Is the distance from the emission point O to the integration point Q,Is the distance from the integration point Q to the receiving point P, k 0 is the wavenumber of free space, and. Δ g is the ground surface impedance on the radio wave propagation path, and is obtained by substituting the electrical parameter of the radio wave propagation path into equation (4).
(4)
Wherein,Is the relative dielectric constant of the ground,Is the ground conductivity.
(5)
Wherein ζ' (·) represents the first derivative of the terrain height; the integral term in the electric wave attenuation function W (P) is defined as I, i.e.:
(6)
the integral equation given by equation (2) is AndThere are two singular points, so the integration I is decomposed into three integrated sums by a segmented integration method, and the numerical solution of I is expressed as follows in combination with the discrete grid defined in the step 1:
(7)
Wherein, 、、As shown in the formulas (8) to (10), respectively.
(8)
(9)
(10)
Wherein,、、Respectively represent、、Is a constant value.
Step 3, calculating by adopting a 5-point Legend-Gaussian integration method and combining a flat ground formulaIs a function of the integral of (a).
For the integral I 1 where the singular point s=0, a gaussian integral solution is used.
Make variable substitution, orderWhereinIs an introduced substitution variable, s is substituted into formula (8), thenIs converted into (a) integration interval：
(11)
At this time, I 1 integrates over the integration intervalNo longer contains any singularities.
The integration of equation (11) is calculated using a 5-point Legend-Gauss product, as shown in equation (12).
(12)
Wherein the method comprises the steps of，AndRespectively, represent the integration interval asWeight factor and sampling point in the schlempe-gaussian integration formula, k=1, 2,3,4,5.
Table 15 Legend-GaussAndIs of the value of (2)
In formula (12)The method is calculated by adopting a flat ground formula, and comprises the following specific steps:
I. Calculation using flat ground equation (13) Electric field strength at：
(13)
Wherein ω is the angular frequency of the low frequency electric wave and has ω=2pi f, f is the frequency of the low frequency electric wave, μ 0 is the air permeability, I is the current of the vertical electric dipole, dl is the charge spacing of the vertical electric dipole, h is the ground clearance of the vertical electric dipole on the z axis,Representation pointsThe terrain elevation at which the position is located,Wavenumber of ground，Is the dielectric constant of air; d 1 and d 2 are defined as:
(14)
(15)
in formula (13) The intermediate parameter is represented and given by equation (16).
(16)
F (P 2) is the Fresnel integral, defined as:
(17)
II. utilizing electric field strength And a wave attenuation functionThe relation between them, equation (18), is calculatedIs a value of (2);
(18)
Wherein, The electric field intensity when the ground is regarded as a conductive plane is defined as:
(19)。
Step 4, based on the Simpson formula, adopting quadratic interpolation polynomial approximation and cubic interpolation polynomial approximation to solve And (5) integrating.
In particular, fromTo the point ofA kind of electronic deviceTotal integration is given byThe integral intervals are calculated in the following two cases:
When (when) When the number is even, a simpson formula based on quadratic interpolation polynomial approximation is adopted to solveIs a function of the integral of (a).
When (when)When odd, forFront of (2)The quadratic interpolation polynomial approximation-based simpson formula is adopted for solving the integral intervalsAnd finally solving the integral of the interval by adopting a cubic interpolation polynomial approximation method in the 3 integral subintervals.
I. when (when)In the case of even number, theThe integral intervals are arranged asGroup jCovering two consecutive integral intervals; defining the integrated function in equation (9) as:
(20)
the simpson formula is used for formula (9), and there are:
(21)
In group j, the integrand Fitting by adopting a quadratic interpolation polynomial:
(22)
Wherein, Is an interpolation nodeConstant function at, and。
Base function of interpolationGiven by equation (23).
(23)
Substituting formula (22) into formula (21) to obtain:
(24)
Wherein:
(25)
(26)
(27)
Wherein, 、Superscript of A, representing the size of the 2j and 2j+1 th gridsRepresents the j-th group，。
II. WhenWhen the number is odd, the integration interval is calculated by adopting a Simpson formula based on quadratic interpolation polynomial approximationI 2 integral of (1), at this timeHas n-5 integral intervals and is arranged asA group. For the followingThe I 2 integral of the interval is then solved approximately using a cubic polynomial.
Thus, I 2 in equation (9) is split into two parts:
(28)
The integral can be solved by using the quadratic polynomial approximation of equation (24) 。
For integrationIs calculated using a cubic interpolation polynomial approximation, so willThe substitution is as follows:
(29)
Definition of the definition ThenThe expression of (2) is as follows:
(30)
Substituting formula (29) The method comprises the following steps:
(31)
Wherein:
(32)
(33)
(34)
(35)
。
Wherein, 、AndThe sizes of the n-3, n-2 and n-1 grids are indicated, respectively.
Step 5, solving by adopting a quadratic interpolation polynomial approximation and transformation integral limit methodIs a function of the integral of (a).
I 3 in equation (10) contains not onlySingular points at, and also including unknown functionsSince the Gaussian integration method cannot be adopted, the method is to first useSeparating.
Order the (36)
When (when)When the method is used, the following steps are included:
(37)
where ζ "(. Cndot.) represents the second derivative of the terrain elevation.
For a pair ofIn the integration intervalPerforming quadratic interpolation polynomial fitting to obtain:
(38)
Wherein, Indicating the size of the n-1 th grid,Indicating the size of the nth grid. To eliminate the singular point in the I 3 integrand, the following variables are used instead:
(39)
Where y represents a substituted variable, and substituting formula (39) into formula (10) yields:
(40)
Substituting equation (38) into equation (40) and directly integrating:
(41)
Wherein:
(42)
(43)
Step 6. By integration 、、Solving the electric wave attenuation function at the receiving point PAnd the total electric field E z in the whole calculation area and the secondary time delay t w generated by the propagation of the low-frequency electric wave on the complex path are obtained through the electric wave attenuation function.
Substituting the formula (7) in the step 2, the formula (12) in the step 3, the formulas (24), (28) in the step 4 and the formula (41) in the step 5 into the formula (2) to obtain a recurrence formula:
(44)
After the initial solution is calculated by the flat ground formula in the step 3, the iterative calculation of the subsequent point sequence is performed 。
Obtaining an electric wave attenuation function W (P) at each point in the calculation region by using the formula (44), substituting the electric field intensity given in step 3And a wave attenuation functionThe relation between them, equation (18), yields the total electric field over the whole calculation region.
Secondary time delay generated by low-frequency wave propagation on complex pathThe calculation formula of (2) is as follows: ; wherein, For the phase of the electric field as it propagates over a complex path,Is the phase of the electric field as it propagates on the good conductor path.
The method can realize modeling of the long-distance complex path by utilizing the mixed grid technology, improves the calculation efficiency of predicting the low-frequency wave propagation characteristic on the complex path by using the integral equation method, ensures the calculation precision, and expands the applicability of the integral equation method in predicting the low-frequency wave propagation characteristic on the complex path.
In addition, in order to verify the effectiveness of the method of the present invention, two specific examples are given below.
Example 1: and (5) predicting the ground field intensity and the secondary time delay of the Gaussian mountain landform.
The radiation power of the excitation source is 1kw, and the frequency of the signal source is 100kHz.
As shown in fig. 3 and 4, the calculated area sizes of the path 1 and the path 2 are bothAnd the relative dielectric constants of the two paths are epsilon=13, and the conductivities are sigma=3× -3 S/m.
Path 1 and path 2 have a maximum height h=1 km and h=2km Gao Sishan at x c =50 km, respectively, the mountain width L is 5km, the height function is: ; wherein, Is a parameter for controlling the width L of a real mountain,。
In example 1, three different grids were used for simulation by the integral equation method, respectively: (1) Grid sizeA fine grid integral equation method of 100 m; (2) Grid sizeA coarse grid integral equation method of 1000 m; (3) The improved integral equation method of the coarse grid (the grid size is 1000 m) and the fine grid (the grid size is 100 m) are mixed, wherein the improved integral equation method adopted by the invention applies the fine grid around the mountain, the rest space uses the coarse grid, and in order to avoid approximate overfitting of a polynomial, a grid smooth transition region with the length of 5km is used at the boundary of the fine and coarse grids.
Fig. 3 and 4 are graphs showing the comparison of the electric field intensity E z and the quadratic time delay t w of the improved integral equation method according to the present invention with the conventional integral equation method (fine grid and coarse grid) and the time domain finite difference method on two gaussian mountain paths (for distinguishing curves clearly, the electric field intensity and the secondary time delay of path 1 are shifted by-20 dB and-1 mus, respectively).
For better visualization, the differences between each of the integral equation method and the time domain finite difference method, E z and t w, are plotted in fig. 3 and 4 (as in the lower part of fig. 3 and 4). The difference is defined asWhereinRepresenting the numerical results of each integral equation method calculation,Is a reference solution obtained using a finite difference method. As can be seen from fig. 3 and 4, the three integration equation schemes provide good results in calibrating the amplitude according to the time domain finite difference method. Whereas for phase, the integral equation method with fine grid and the improved integral equation method with mixed grid have higher prediction accuracy than the integral equation method with coarse grid, and the difference is particularly obvious in the case of higher Gao Sishan.
Example 2: and predicting the ground field intensity and the secondary time delay of the triangular mountain landform.
As shown in fig. 5 and 6, the parameters in path 3 and path 4 are the same as those used in path 1 and path 2 in example 1, respectively, except for the shape of the mountain. The height function of the triangular mountains employed by path 3 and path 4 is as follows:
；
Wherein, 。
The three integral equation methods employed in example 2 are also the same as those described in example 1.
As can be seen from fig. 5 and 6, the integral equation method with fine meshes and the modified integral equation method are very identical in both amplitude and phase to the result of the finite difference method in the time domain. However, in the case of terrain with significant discontinuities, the integral equation with coarse grid is not well-matched to the result of the time-domain finite difference method due to poor discretization.
Table 2 comparison of the calculation costs of the original integral equation method using the fine mesh and the coarse mesh, respectively, with the improved integral equation method using the hybrid mesh according to the present invention
Table 2 lists the memory requirements and computation time required for each of the integral equation methods simulations.
Since the mountain model used by paths 1 through 4 has the same width and position, the grid discretization on each path is the same, which results in nearly the same computational cost on each path.
Thus, path 1 is chosen as an example to show the computational consumption of each integral equation method in a simulation.
As shown in Table 2, under the same calculation precision, compared with a fine-grid integral equation method, the simulation time can be reduced by 8 times and the memory requirement can be reduced by 6 times by adopting the improved integral equation method. In addition, the improved integral equation method provided by the invention has higher calculation accuracy than the coarse grid integral equation method under the condition of not increasing too much calculation cost.
The foregoing description is, of course, merely illustrative of preferred embodiments of the present invention, and it should be understood that the present invention is not limited to the above-described embodiments, but is intended to cover all modifications, equivalents and alternatives falling within the spirit and scope of the present invention as defined by the appended claims.
Claims (1)
1. A method for predicting the propagation characteristics of low-frequency electric waves in a complex path based on an integral equation is characterized in that,
The method comprises the following steps:
Step 1, inputting a model file, wherein the content of the input model file comprises grid parameters of a calculation area, excitation source parameters, electric parameters of an electric wave propagation path and electric parameters of a free space;
Step 2, decomposing an integral term contained in the electric wave attenuation function W (P) at the receiving point P into the following three sections by utilizing a sectional integral method: an integral I 1 containing the singular point 0, a main integral I 2, and an integral I 3 containing the singular point x n and the unknown integral W (x n);
Step 3, calculating the integral of I 1 by adopting a 5-point Legend-Gaussian integral method and combining a flat ground formula;
step 4, for I 2 integration with n-2 integration sub-intervals, the following two cases are calculated, namely:
When n-2 is even, solving the integral of I 2 by adopting a Simpson formula based on quadratic interpolation polynomial approximation;
When n-2 is an odd number, solving a simpson formula based on quadratic interpolation polynomial approximation for n-5 integral intervals before I 2, and solving integral of intervals for the last 3 integral intervals of I 2 by adopting a method of cubic interpolation polynomial approximation;
Step 5, solving the integral of I 3 by adopting a quadratic interpolation polynomial approximation and transformation integral limit method;
step 6, solving an electric wave attenuation function W (P) at a receiving point P through an integral I 1、I2、I3, and obtaining a total electric field in the whole calculation area and a secondary time delay generated by low-frequency electric wave propagation on a complex path through the electric wave attenuation function;
The step1 specifically comprises the following steps:
In a two-dimensional rectangular coordinate system (x, z), a length of a calculation region is defined as x T, and x T is discretized into n grids, each of which has a length of Δx t, and satisfies the following relationship:
Wherein x represents a lateral distance, z represents a longitudinal height, x l represents an integration length of the first integration point in the x direction, l=1, 2,..n, t=1, 2,..l, and an initial value x 0 =0 of x l, and subscript t represents the t-th integration subinterval;
The method comprises the steps of adopting mixed grids aiming at different terrains on a low-frequency electric wave propagation path, namely adopting grids with different sizes for different terrains, and adopting relatively thicker or relatively thinner grids when the propagation path is a flat ground or a fluctuant mountain;
The excitation source is a vertical electric dipole, the current and charge distance of the excitation source are I, dl respectively, and the electric parameters of the electric wave propagation path comprise the relative dielectric constant epsilon r and the electric conductivity sigma of the earth surface; the electrical parameters of free space include air permittivity epsilon 0 and air permeability mu 0;
The step 2 specifically comprises the following steps:
The electric wave attenuation function W (P) represented by Hufford integral equation at the receiving point P is given by formula (2):
Wherein P is the position of a receiving point, the coordinate of the P is (x T,ξ(xT)), Q is the position of a topographic surface area dividing point, the coordinate of the Q is (s, ζ (s)), s represents the distance coordinate of the Q point in the two-dimensional coordinate system (x, z) along the x direction, ζ (·) represents the topographic height, W (Q) represents the electric wave attenuation function of the Q point, r 0 is the distance from the transmitting point O to the receiving point P, i represents an imaginary unit, and λ is the wavelength of free space; f (x T, s) is defined as:
where r 1 is the distance from the emission point O to the integration point Q, r 2 is the distance from the integration point Q to the reception point P, k 0 is the wave number of free space, and k 0 =2pi/λ;
δ g is the ground surface impedance on the radio wave propagation path, obtained by substituting the electrical parameter of the radio wave propagation path into equation (4);
Wherein ε r is the ground relative permittivity and σ is the ground conductivity;
A＝-ξ(xT)+ξ(s)+ξ'(s)·(xT-s) (5)
Wherein ζ' (·) represents the first derivative of the terrain height; the integral term in the electric wave attenuation function W (P) is defined as I, i.e.:
The integral equation given by equation (2) has two singular points at s=0 and s=x T, so the integral I is decomposed into the form of the sum of three integrals using a piecewise integration method, and the numerical solution of I is expressed in combination with the discrete grid defined in step 1:
Wherein, the definition of I 1(xn)、I2(xn)、I3(xn) is shown in the formulas (8) to (10);
Wherein, [0, x 1]、(x1,xn-1]、(xn-1,xn ] each represent the integration interval of I 1(xn)、I2(xn)、I3(xn);
The step 3 specifically comprises the following steps:
For the integral I 1 with the singular point s=0, solving by using a Gaussian integral method; make variable substitution, order Where y is the introduced substitution variable, substituting s into equation (8), the integration interval of I 1 is converted to [ -1,1];
At this time, the I 1 integral no longer contains any singular point in the integral interval [ -1,1 ];
the integration of equation (11) is calculated by 5-point Legend-Gauss product, as shown in equation (12);
Wherein the method comprises the steps of Lambda K and y K represent the weight factors and sampling points in the Legendre-Gaussian integration formula when the integration interval is [ -1,1], K=1, 2,3,4,5, respectively; w (s K) is calculated by adopting a flat ground formula, and the specific method is as follows:
I. The electric field strength E z at s K is calculated using flat ground equation (13):
Wherein ω is an angular frequency of the low-frequency electric wave and ω=2pi f, f being a frequency of the low-frequency electric wave;
mu 0 is air permeability, I is current of the vertical electric dipole, dl is charge interval of the vertical electric dipole, h is ground clearance of the vertical electric dipole on a z axis, and xi (s K) represents terrain height at a point s K;
k g is the wavenumber of the ground and Wherein ε 0 is the air dielectric constant;
d 1 and d 2 are defined as:
p 2 in equation (13) represents an intermediate parameter, given by equation (16);
F (P 2) is the Fresnel integral, defined as:
II, calculating the value of W (s K) by using a formula (18) which is the relation between the electric field intensity E z and the electric wave attenuation function W (P);
Ez＝E0W(P) (18)
wherein, E 0 is the electric field strength when the ground is regarded as a conductive plane, and is defined as:
the step 4 specifically comprises the following steps:
I. When n-2 is an even number, the n-2 integral intervals are arranged in j= (n-2)/2 groups, and the J-th group (x 2j-1,x2j,x2j+1) covers two consecutive integral intervals; defining the integrated function in equation (9) as:
the simpson formula is used for formula (9), and there are:
in group j, the multiplicative function f(s) is fitted with a quadratic interpolation polynomial:
f(s)＝y2j-1l2j-1(s)+y2jl2j(s)+y2j+1l2j+1(s) (22)
Wherein y m＝f(xm) is a constant function at interpolation node x m, and m=2j-1, 2j,2j+1;
The interpolated basis function l m(s) is given by equation (23);
Substituting formula (22) into formula (21) to obtain:
Wherein:
Wherein Δχ 2j、Δx2j+1 represents the size of the 2j and 2j+1 th grids, the superscript [ J ] of a represents the J-th group (x 2j-1,x2j,x2j+1), j=1, 2,;
II, when n-2 is an odd number, calculating I 2 integration with an integration interval of (x 1,xn-4) by adopting a Simpson formula based on quadratic interpolation polynomial approximation, wherein (x 1,xn-4) has n-5 integration subintervals and is arranged into J' = (n-5)/2 groups;
thus, I 2 in equation (9) is split into two parts:
The integral can be solved by using the quadratic polynomial approximation of equation (24)
For integrationCalculated using a cubic interpolation polynomial approximation, f(s) is replaced by:
f(s)＝yn-4ln-4(s)+yn-3ln-3(s)+yn-2ln-2(s)+yn-1ln-1(s) (29)
Definition m' =n-4, n-3,..n-1, then the expression of l m'(s) is as follows:
Substituting formula (29) The method comprises the following steps:
Wherein:
Wherein K' =Δx n-3+Δxn-2+Δxn-1;
Deltax n-3、Δxn-2 and Deltax n-1 represent the sizes of the n-3, n-2 and n-1 grids, respectively; the step 5 specifically comprises the following steps:
Order the
When s=x n, there are:
where ζ "(. Cndot.) represents the second derivative of terrain elevation;
Performing quadratic interpolation polynomial fitting on W(s) g (x T, s) in an integration interval (x n-1,xn) to obtain:
Wherein Δx n-1 represents the size of the n-1 th mesh, and Δx n represents the size of the n-th mesh;
To eliminate the singular point in the I 3 integrand, the following variables are used instead:
Where y represents a substituted variable, and substituting formula (39) into formula (10) yields:
Substituting equation (38) into equation (40) and directly integrating:
Wherein:
the step 6 specifically comprises the following steps:
Substituting the formula (7) in the step 2, the formula (12) in the step 3, the formulas (24), (28) in the step 4 and the formula (41) in the step 5 into the formula (2) to obtain a recurrence formula:
after the initial solution is calculated by the flat ground formula in the step 3, the W (x n) is calculated by iteration of the subsequent point sequence;
obtaining an electric wave attenuation function W (P) at each point in the calculation region by using a formula (44), and substituting the electric field intensity E z given in the step 3 and the electric wave attenuation function W (P) into a relation, namely a formula (18), to obtain a total electric field in the whole calculation region;
The calculation formula of the secondary time delay t w generated by the propagation of the low-frequency electric wave on the complex path is as follows: Wherein/> For the phase of the electric field as it propagates over a complex path,/>Is the phase of the electric field as it propagates on the good conductor path.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410381119.1A CN117992707B (en) | 2024-04-01 | 2024-04-01 | Method for predicting low-frequency electric wave propagation characteristics in complex path based on integral equation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202410381119.1A CN117992707B (en) | 2024-04-01 | 2024-04-01 | Method for predicting low-frequency electric wave propagation characteristics in complex path based on integral equation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117992707A CN117992707A (en) | 2024-05-07 |
| CN117992707B true CN117992707B (en) | 2024-06-14 |
Family
ID=90889225
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202410381119.1A Active CN117992707B (en) | 2024-04-01 | 2024-04-01 | Method for predicting low-frequency electric wave propagation characteristics in complex path based on integral equation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117992707B (en) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105184070A (en) * | 2015-08-31 | 2015-12-23 | 华南理工大学 | Piecewise fitting method for calculating volt-second characteristic curve based on voltage integration method |
| CN105740204A (en) * | 2016-03-14 | 2016-07-06 | 西安理工大学 | Low-frequency-band earth conductivity rapid inversion method under irregular terrain |
Family Cites Families (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| ITTO20030675A1 (en) * | 2003-09-05 | 2005-03-06 | Torino Politecnico | PROCEDURE FOR THE NUMERICAL MODELING OF SIZES |
| CN101937026B (en) * | 2010-07-30 | 2013-04-17 | 西安理工大学 | Method for predicating ground wave propagation attenuation factor at high precision |
| CN107341284B (en) * | 2017-05-18 | 2020-11-17 | 西安理工大学 | Bidirectional parabolic equation method for predicting low-frequency radio wave propagation characteristics with high precision |
| CN107966679A (en) * | 2017-11-16 | 2018-04-27 | 中国电子科技集团公司第二十研究所 | Very low frequency navigation Real-time Network based on cloud framework is formatted radio wave propagation modification method |
| CN116756575B (en) * | 2023-08-17 | 2023-11-03 | 山东科技大学 | Non-invasive load decomposition method based on BGAIN-DD network |
| CN116758104B (en) * | 2023-08-21 | 2023-11-14 | 山东科技大学 | Multi-instance portrait matting method based on improved GCNet |
| CN116992192B (en) * | 2023-09-28 | 2023-12-12 | 山东科技大学 | Low-frequency electric wave propagation method in sea ice mixing path based on parabolic side Cheng Yuce |
-
2024
- 2024-04-01 CN CN202410381119.1A patent/CN117992707B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105184070A (en) * | 2015-08-31 | 2015-12-23 | 华南理工大学 | Piecewise fitting method for calculating volt-second characteristic curve based on voltage integration method |
| CN105740204A (en) * | 2016-03-14 | 2016-07-06 | 西安理工大学 | Low-frequency-band earth conductivity rapid inversion method under irregular terrain |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117992707A (en) | 2024-05-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Apaydin et al. | Numerical investigations of and path loss predictions for surface wave propagation over sea paths including hilly island transitions | |
| Apaydin et al. | The split-step-fourier and finite-element-based parabolic-equation propagation-prediction tools: Canonical tests, systematic comparisons, and calibration | |
| CN112949134B (en) | Earth-well transient electromagnetic inversion method based on non-structural finite element method | |
| CN107341284B (en) | Bidirectional parabolic equation method for predicting low-frequency radio wave propagation characteristics with high precision | |
| CN109684740B (en) | Electromagnetic multi-scale calculation method based on mixed grid and time step | |
| CN112865883A (en) | Channel modeling method combining ray tracing method and time domain finite difference method | |
| CN111537987A (en) | Genetic algorithm-based incoherent scattering radar data inversion method | |
| CN112733364A (en) | Foil strip cloud scattering rapid calculation method based on impedance matrix blocking | |
| CN114707458A (en) | Sensor deployment method suitable for monitoring partial discharge of switch cabinet | |
| CN117992707B (en) | Method for predicting low-frequency electric wave propagation characteristics in complex path based on integral equation | |
| CN111929348A (en) | Method for calculating surface potential of direct current grounding electrode in horizontal multilayer structure earth environment | |
| CN116992192B (en) | Low-frequency electric wave propagation method in sea ice mixing path based on parabolic side Cheng Yuce | |
| Ozyalcin et al. | A novel TLM-based time-domain wave propagator | |
| Baba et al. | The effect of localized ionospheric perturbations on subionospheric VLF propagation on the basis of finite element method | |
| CN109614732B (en) | Electromagnetic compatibility modeling method and device for object | |
| CN116430126B (en) | Electromagnetic background cognition-based electromagnetic silence target detection method and device and computer equipment | |
| Marchenko et al. | Finite-difference time-domain analysis of ELF radio wave propagation in the spherical Earth–ionosphere waveguide and its validation based on analytical solutions | |
| Yilmaz et al. | A single-boundary implicit and FFT-accelerated time-domain finite element-boundary integral solver | |
| Sewell et al. | Extracting modal field profiles from 3D unstructured transmission line modelling meshes for use as sources and observers | |
| Omar et al. | Fast corona discharge solver for precipitators using multi‐grid methods on fine grids | |
| Yee et al. | Impedance boundary condition simulation in the FDTD/FVTD hybrid | |
| CN118013172B (en) | Low-frequency ground wave propagation characteristic prediction method for solving integral equation based on trapezoidal formula | |
| Nutaro et al. | An event driven, simplified TLM method for predicting path-loss in cluttered environments | |
| Kato et al. | FDTD simulations of LEMP propagation in the earth-ionosphere waveguide using different lightning models | |
| Aydemir et al. | Genetic approach to the minimization of the coupling between aircraft antennas |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |