A method for calculating tropospheric scatter transmission loss and propagation delay

Abstract

A method for calculating tropospheric scatter transmission loss and propagation delay is as follows: Obtain meteorological data on the scatter propagation path through a numerical meteorological model; Calculate the electromagnetic wave propagation path in the main axis direction of the transmitting antenna and the electromagnetic wave propagation path in the main axis direction of the receiving antenna based on the meteorological data; Dissect the transmitting beam and the receiving beam, and intercept the path from the transmitting antenna to the scatter point and then to the receiving antenna from the ray paths of the obtained transmitting sub-beams and the ray paths of the receiving sub-beams, which is the tropospheric scatter propagation path of the scatter sub-beams; Calculate the volume of the common scatterers of each sub-beam; Calculate the received power of the scatter sub-beams and the propagation time of the scatter sub-beams; Calculate the tropospheric scatter transmission loss and propagation delay. The present invention can solve the problem that the existing tropospheric scatter model cannot accurately reflect the influence of atmospheric environment changes in different regions and at different times on the transmission loss and propagation delay, and can calculate the tropospheric scatter loss and propagation delay more accurately.

Description

Technical Field

[0001] The present invention belongs to the technical field of tropospheric scatter communication, and particularly relates to a calculation method for tropospheric scatter transmission loss and transmission delay. Background Art

[0002] The tropospheric scatter communication link has the deficiencies of large transmission loss, serious multipath effect, and being significantly affected by the atmospheric environment. Therefore, in order to improve the system performance of the application of the tropospheric scatter link, it is necessary to accurately master the transmission loss characteristics and propagation delay characteristics of the tropospheric scatter link. In addition to the semi-empirical estimation model, various loss prediction models have been generated based on the scatter propagation mechanism, such as the tropospheric scatter transmission loss calculation method based on the generalized scatter cross-section proposed by Zhang Minggao, the tropospheric scatter transmission loss calculation method based on the parabolic equation method proposed by Li Lei, the tropospheric scatter channel model based on rays proposed by Ergin Dinc, etc.

[0003] There is also an obvious multipath effect in tropospheric scatter. The electromagnetic wave reaches the receiver through different paths, and there are propagation delay and time delay spread in the received signal. The existing models have insufficient characterization ability of meteorological characteristics. Usually, it is assumed that the electromagnetic wave propagates at the speed of light in a straight line. For example, Sunde derived a rough calculation formula for the average time delay difference based on the symmetric link and the smooth spherical surface starting from the link geometry; Bello derived a calculation method for the scattered time delay power spectrum based on the scattering cross-section theory of Boor-Gordon and the two-dimensional plane hypothesis starting from the geometry of the scattering link on the basis of the turbulent incoherent scatter propagation mechanism; Zhang Minggao derived the normalized delay power spectrum based on the generalized scatter cross-section theory combined with the geometry of the scattering link; Ergin Dinc's group delay model for the tropospheric scatter channel based on rays, etc.

[0004] The above models have played an important role in assisting the design of the tropospheric scatter communication link and the determination of the parameters of the scatter system equipment. However, with the continuous expansion of the application of tropospheric scatter, the time synchronization system and the over-the-horizon detection system based on the tropospheric scatter link have put forward higher requirements for the prediction and analysis capabilities of accurately predicting and analyzing the transmission loss and propagation delay under different meteorological conditions. In the existing transmission loss and propagation delay estimation models, the meteorological parameters generally use the empirical fitting values or the statistical result averages of the meteorological parameters to describe, and cannot give the specific meteorological environment information in different regions and at different times; moreover, it is assumed that the electromagnetic wave propagates in a straight line at a constant speed of light in the atmosphere during the calculation process, and it is impossible to effectively reflect the path bending, propagation speed change, and atmospheric absorption attenuation that occur when the electromagnetic wave propagates in the complex and changeable atmospheric environment. There is still room for improvement in the calculation accuracy of the tropospheric scatter transmission loss and propagation delay. Summary of the Invention

[0005] The object of the present invention is to provide a method for calculating tropospheric scatter transmission loss and propagation delay based on a numerical weather model, which can calculate the tropospheric scatter transmission loss and propagation delay more accurately.

[0006] To achieve the above object, the present invention adopts the following technical solution:

[0007] A method for calculating tropospheric scatter transmission loss and propagation delay, comprising the following steps:

[0008] S1. Obtain meteorological data on the scatter propagation path through a numerical weather model; the numerical weather model stores meteorological data in a grid form, with a grid interval in the horizontal direction being the longitude and latitude resolution, and the vertical grid being an isobaric surface. Based on the numerical weather model, interpolation or extrapolation methods are used to obtain meteorological data at a specified location, and the meteorological data includes temperature, specific humidity, air pressure, and wind speed;

[0009] S2. Calculate the electromagnetic wave propagation path in the main axis direction of the transmitting antenna and the electromagnetic wave propagation path in the main axis direction of the receiving antenna based on the meteorological data obtained through the numerical weather model;

[0010] S3. Perform cross-sectional and longitudinal dissections and path rotations on the transmitting beam to obtain the ray paths of each transmitting sub-beam, perform cross-sectional dissection, path rotation, and ray alignment on the receiving beam to obtain the ray paths of the receiving sub-beams within the effective range, and intercept the path from the transmitting antenna to the scatter point and then to the receiving antenna from the obtained ray paths of the transmitting sub-beams and the ray paths of the receiving sub-beams within the effective range, which is the tropospheric scatter propagation path of the scatter sub-beam;

[0011] S4. Calculate the volume of the common scatterer of each sub-beam, where the common scatterer of the sub-beam is the intersecting part of the transmitting sub-beam and the receiving sub-beam;

[0012] S5. Calculate the received power of the scatter sub-beam and the propagation time of the scatter sub-beam;

[0013] S6. Calculate the tropospheric scatter transmission loss and propagation delay;

[0014] The tropospheric scatter transmission loss is: where P r_ij is the received power of the scatter sub-beam with the serial number subscript ij, and P t is the transmitted power of the transmitting antenna;

[0015] Normalize the received power of all scatter sub-beams and arrange the propagation times of all scatter sub-beams in ascending order to obtain the normalized delay power spectrum of the tropospheric scatter link, which is the tropospheric scatter propagation delay.

[0016] More specifically, the step of obtaining meteorological data at a specified location based on a numerical weather model in step S1 is as follows:

[0017] S101. Convert the geodetic height H at the specified location g to the geopotential height wherein, is the latitude at the specified location, is the normal gravity of the earth ellipsoid surface, γ 45 ° is the standard value of the gravitational acceleration at 45° latitude, is the local earth radius;

[0018] S102. Determine the corresponding interpolation method and extrapolation method by comparing the relative position relationship between the interpolation point and the geopotential height of the isobaric surface of the numerical weather model, including the following three cases:

[0019] (1) When the geopotential height of the interpolation point is inside the isobaric surface, determine 8 grid points m1 to m8 adjacent to the interpolation point W in the numerical weather model through the longitude, latitude and geopotential height of the interpolation point W, and interpolate to the geopotential height of the interpolation point for 4 groups of two vertically opposite grid points with the same longitude and latitude: vertically, linear interpolation is used for temperature, specific humidity and wind speed, and exponential model interpolation is used for air pressure; horizontally, bilinear interpolation is used for temperature, specific humidity, wind speed and air pressure;

[0020] (2) When the geopotential height of the interpolation point is below the bottom layer of the isobaric surface, determine 8 grid points m1 to m8 adjacent to the interpolation point W through the longitude, latitude and geopotential height of the interpolation point W. For the 4 grid points of the upper isobaric surface among the 8 grid points, bilinear interpolation method is used horizontally to obtain the meteorological parameters at the corresponding position of the interpolation point on the isobaric surface, and the following methods are used respectively to extrapolate and calculate the meteorological parameters at the corresponding position of the interpolation point on the isobaric surface vertically: exponential model is used for air pressure; temperature T j = T v_j - 0.0065h g_j , T v_j is the temperature obtained by horizontal bilinear interpolation on the isobaric surface; specific humidity q j = q v_j , q v_j is the specific humidity obtained by horizontal bilinear interpolation on the isobaric surface; power law interpolation is used for wind speed;

[0021] (3) When the potential height of the interpolation point is above the highest layer of the isobaric surface, the eight grid points m1 to m8 adjacent to the interpolation point W are determined by the longitude and latitude of the interpolation point W and the potential height. For the four grid points on the lower isobaric surface among the eight grid points, the bilinear interpolation method is used in the horizontal direction to obtain the meteorological parameters at the corresponding positions of the interpolation points on the isobaric surface. In the vertical direction, the meteorological parameters at the corresponding positions of the interpolation points on the isobaric surface are extrapolated using the following methods: the specific humidity is zero; the air pressure is extrapolated using the exponential model based on the horizontal bilinear interpolation; the temperature and wind speed components are obtained using the CIRA86 international reference atmospheric model.

[0022] More specifically, in step S2, a ray tracing method is used to calculate the electromagnetic wave propagation path in the main axis direction of the transmitting antenna, and the steps are as follows:

[0023] S201, taking the location of the transmitting antenna as the coordinate origin, establish a local rectangular coordinate system xyz based on the right-hand spiral criterion, the x-axis of the coordinate system is the projection of the initial direction of the ray on the tangent plane of the earth's surface, the z-axis is the direction vector from the center of the earth to the location of the transmitter, α0 and β0 are the azimuth and elevation angles of the main axis of the transmitting antenna respectively;

[0024] In this local rectangular coordinate system xyz, the electromagnetic wave propagation path satisfies the linear differential equation system Where n is the atmospheric refractive index and N is the atmospheric refractive index;

[0025] The above linear differential equations have the following initial boundary conditions: Under the constraints and the height limit of the troposphere top, the propagation path of electromagnetic waves in the main axis direction of the transmitting antenna can be obtained by iteratively solving the meteorological data obtained from the numerical meteorological model through the numerical method;

[0026] For the receiving antenna, the same steps are used to calculate with the receiving antenna as the starting point to obtain the electromagnetic wave propagation path in the direction of the antenna's main axis.

[0027] More specifically, in step S3, the transmit beam is split by using a rectangular cross-section splitting method, and the transmit beam is split into K equal parts at equal angle intervals dw along the horizontal axis and the vertical axis of the transmit beam cross section, respectively, and the transmit beam is decomposed into a plurality of closely arranged and non-overlapping regular tetrahedral transmit sub-beams;

[0028] The electromagnetic wave propagation path l in the main axis direction of the transmitting antenna transmit_ray_main Convert from the local rectangular coordinate system xyz to the spherical coordinate system BLH, and then rotate the initial pointing direction of the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system BLH to the initial pointing direction of each transmitting sub-beam ray. The ray path l of the transmitting sub-beam with the serial number subscript ij transmit_ray_ij = l transmit_ray_main(B, L, H) + (dw·i, dw·j, 0), l transmit_ray_main (B, L, H) is the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system BLH, and i, j = -K / 2, -K / 2 + 1, …, K / 2 - 1, K / 2.

[0029] More specifically, in step S3, the receiving beam is dissected by the transverse dissection method. The receiving beam is transversely dissected into K equal parts at equal angular intervals dw in the cross-section of the receiving beam, obtaining several receiving sub-beams with the same azimuth angle but different elevation angles; the electromagnetic wave propagation path l of the main axis direction of the transmitting antenna transmit_ray_main is converted from the local rectangular coordinate system xyz to the spherical coordinate system BLH, and then the initial direction of the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system BLH is rotated to the initial direction of each transmitting sub-beam ray. The ray path l of the transmitting sub-beam with the serial number subscript ij transmit_ray_ij = l transmit_ray_main (B, L, H) + (dw·i, dw·j, 0), l transmit_ray_main (B, L, H) is the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system BLH, i = -K / 2, -K / 2 + 1, …, K / 2 - 1, K / 2; for each transmitting sub-beam, calculate the minimum distance point between the receiving sub-beam ray and the transmitting sub-beam ray and the minimum distance between these two minimum distance points, obtain the angle between these two minimum distance points with respect to the transmitting point of the receiving antenna, rotate the azimuth angle of the receiving sub-beam ray according to the size of the angle until the minimum distance between the two minimum distance points is less than the set threshold, calculate the angle between the initial azimuth angle of the receiving sub-beam ray and the azimuth angle of the main axis of the receiving antenna at this time, and discard the receiving sub-beam rays whose azimuth angle exceeds the beam width, obtaining the ray paths of the receiving sub-beams within the effective range.

[0030] More specifically, in step S4, the volume of the sub-beam common scatterer is calculated through the following steps:

[0031] In the plane formed by the scatter point S, the transmitting point T, and the receiving point R, starting from the scatter point S, make sub-beam ray path tangent vectors in the directions of the transmitting point T and the receiving point R respectively. The lengths of the tangent vectors are respectively set as the actual ray path lengths from the transmitting point T to the scatter point S and from the scatter point S to the receiving point R. The end points T' and R' of the tangent vectors are the virtual transmitting point and the virtual receiving point;

[0032] In the plane formed by the virtual emission point T', the scattering point S, and the virtual reception point R', with the virtual emission point T' as the origin o', the direction from the virtual emission point T' to the scattering point S as the x'-axis, and the direction from the center of the earth to the virtual emission point T' as the y'-axis, a three-dimensional rectangular coordinate system x'y'z' is constructed based on the right-hand screw rule. In this three-dimensional rectangular coordinate system x'y'z', T'S and SR' are the virtual straight propagation paths of the electromagnetic wave. Rotate T'S and SR' around the origin o' and the virtual reception point R' by 0.5dw and -0.5dw respectively in the x'o'y' plane to obtain the upper and lower boundaries of the transmitted sub-beam and the received sub-beam, and further obtain the coordinates of the intersection points P1 - P4 of the upper and lower boundaries of the transmitted sub-beam and the received sub-beam in x'y'z'; the coordinates of each vertex of the sub-beam common scatterer can be obtained through the coordinates of P1 - P4 in x'y'z'; divide the sub-beam common scatterer into several trapezoidal prisms along the x'-axis direction, calculate the volume of each trapezoidal prism, and the sum of the volumes of all trapezoidal prisms is the volume of the sub-beam common scatterer.

[0033] More specifically, in step S5, the steps of calculating the received power of the scattered sub-beam using the bistatic radar equation are as follows:

[0034] Calculate the received power of each trapezoidal prism. The received power of the l-th trapezoidal prism in the scattered sub-beam is: Where P t_ij is the transmitted power of the scattered sub-beam with the serial number subscript ij, G t , G r are the transmitting antenna gain and the receiving antenna gain respectively, g t , g r are the directivity function of the transmitting antenna and the directivity function of the receiving antenna respectively, ξ is the wavelength of the electromagnetic wave, dV ij_l is the volume of the l-th trapezoidal prism, R ij_l is the distance from the center point of the l-th trapezoidal prism to the receiving antenna, S ij_l is the distance from the center point of the l-th trapezoidal prism to the transmitting antenna, σ ij_l is the scattering cross-section of the l-th trapezoidal prism;

[0035] The received power P r_ij of the scattered sub-beam with the serial number subscript ij is the sum of the received powers of all trapezoidal prisms in this scattered sub-beam, P r_ij = ∑P r_ij_l .

[0036] More specifically, in step S5, the refractive index integration method is used to calculate the propagation time τ ij , In the formula, S t is the ray path from the emission point to the scattering point, S ris the ray path from the scattering point to the receiving point, n is the atmospheric refraction index on the ray path, c is the speed of light, and ds is the curve element.

[0037] Preferably, in step S6, first based on the transmission loss L of the sub-beam ij correct the received power of the scattered sub-beam, and the received power of the corrected scattered sub-beam is: P in the formula r_ij is the received power of the scattered sub-beam with the serial number subscript ij, and P t_ij is the transmitted power of the scattered sub-beam with the serial number subscript ij, is the atmospheric absorption loss of the scattered sub-beam with the serial number subscript ij, and then according to the received power P of the corrected scattered sub-beam r ' _ij calculate the tropospheric scatter transmission loss, and the tropospheric scatter transmission loss

[0038] Preferably, according to the received power P of the corrected scattered sub-beam r ' _ij calculate the tropospheric scatter transmission delay.

[0039] As can be seen from the above technical solutions, the present invention applies the numerical weather model to the calculation of tropospheric scatter transmission loss and propagation delay, enabling it to have a more detailed and accurate weather analysis ability than existing models; applies the ray tracing method to the calculation of tropospheric scatter paths, and for the first time considers the influence of path bending and delay on transmission loss and propagation delay in the tropospheric scatter model, making the calculation results of propagation delay and transmission loss more accurate and reliable; proposes a method for tropospheric scatter beam dissection and ray alignment based on the turbulent incoherent scattering mechanism, realizing an accurate simulation of the tropospheric scatter process; gives a new method for accurately calculating the volume of the common scatterer, improving the accuracy of transmission loss calculation. Applying the present invention to the tropospheric scatter link can analyze the variation characteristics of transmission loss and propagation delay under different meteorological environments in different regions, providing a reference for the parameter design of tropospheric scatter equipment and the performance analysis of scatter links. Description of the Drawings

[0040] Figure 1 is the flow chart of the method of the present invention;

[0041] Figure 2 is the schematic diagram of the interpolation calculation of meteorological data of the numerical weather model;

[0042] Figure 3 is the schematic diagram of the three-dimensional ray tracing method in the local rectangular coordinate system;

[0043] Figure 4 is the schematic diagram of the emission beam dissection;

[0044] Figure 5 It is a schematic diagram of receiving beam dissection;

[0045] Figure 6 It is a schematic diagram of the ray alignment method;

[0046] Figure 7 It is a schematic diagram of the sub-beam scattering common body;

[0047] Figure 8 It is a schematic diagram of the virtual emission point and the virtual reception point;

[0048] Figure 9 It is a schematic diagram of the subdivided trapezoid of the scattering common body;

[0049] Figure 10 It is the transmission loss of link 2063 under different elevation angle conditions by using the method of the present invention and the methods of ITU-R P.617-2 and ITU-R P.617-5;

[0050] Figure 11 It is the transmission loss of link 2307 under different elevation angle conditions by using the method of the present invention and the methods of ITU-R P.617-2 and ITU-R P.617-5;

[0051] Figure 12 It is a schematic diagram of the geographical location of link 1441;

[0052] Figure 13 It is a schematic diagram of the geographical location of link 2305;

[0053] Figure 14 It is a graph of the hourly variation of the transmission loss of links 1441 and 2305;

[0054] Figures 15a to 15f They are the normalized delay power spectra of 6 scattering links calculated by the method of the present invention and the Bello method respectively;

[0055] Figure 16a and Figure 16b They are the hourly normalized delay power spectrum heat maps of link 1441 and link 2305 on August 25th, 26th and 27th, 2020 respectively.

[0056] The following further describes the specific embodiments of the present invention in detail with reference to the accompanying drawings. Specific Embodiments

[0057] Next, the technical solutions of the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0058] The existing calculation methods for tropospheric scatter transmission loss and propagation delay cannot effectively reflect the influence of meteorological environments in different regions and at different times on the propagation of tropospheric scatter electromagnetic waves, and cannot accurately calculate the path change, speed change, atmospheric absorption attenuation, and received power of electromagnetic waves when propagating in different atmospheric environments, thus affecting the accuracy of the calculation of tropospheric scatter transmission loss and propagation delay. To solve the above problems, the basic idea of the present invention is as follows: First, meteorological parameters of the scatter link calculation area are extracted based on a numerical meteorological model, and the electromagnetic wave propagation paths in the main axis directions of the transmitting antenna and the receiving antenna are calculated with the support of the meteorological parameters; the beam dissection and ray alignment method is used to simulate the scattering process, and the ray paths of the transmitting sub-beam and the receiving sub-beam are obtained by rotating the electromagnetic wave ray path of the main axis of the antenna, and then the path from the transmitting end to the scattering point and then to the receiving end is intercepted as the tropospheric scatter propagation path of each sub-beam; then, the volume of the common scatterer of the sub-beams is calculated, and the received power and propagation time of the sub-beams are calculated based on the volume of the common scatterer of the sub-beams. Finally, the tropospheric scatter transmission loss and propagation delay are calculated according to the received power and propagation time of each sub-beam.

[0059] Figure 1 For the flowchart of the method of the present invention, the following is combined with Figure 1 , to further illustrate the method of the present invention. In the calculation process, relevant parameters of the tropospheric scatter link need to be obtained. These relevant parameters of the tropospheric scatter link include the positions and altitudes of the transmitting antenna and the receiving antenna, the horizontal line-of-sight angles of the transmitting antenna and the receiving antenna, the elevation angles of the transmitting antenna and the receiving antenna, the electromagnetic wave frequency, the beam width, and the calculation time; after obtaining the parameters of the tropospheric scatter link, as Figure 1 shown, the method of the present invention includes the following steps:

[0060] S1. Obtain meteorological data on the scatter propagation path (scatter link calculation area) through a numerical meteorological model; the numerical meteorological model used in the present invention stores meteorological data in a grid form, the grid interval in the horizontal direction is the longitude and latitude resolution of the numerical meteorological model, and the vertical grid is an isobaric surface. Based on the numerical meteorological model, interpolation or extrapolation methods are used to obtain meteorological data at a specified position (scatter link calculation area). These meteorological data include temperature, specific humidity, air pressure, and wind speed (wind speed includes the eastward component of the horizontal wind and the northward component of the horizontal wind).

[0061] The method for obtaining meteorological data at a specified location based on a numerical meteorological model is as follows:

[0062] S101. Convert the geodetic height H at the specified location g to the geopotential height Numerical meteorological models generally use the geopotential height system, so it is necessary to convert the geodetic height and the geopotential height. The conversion formula is: where is the latitude at the specified location, is the normal gravity of the earth ellipsoid surface, γ 45 ° is the standard value of the gravitational acceleration at latitude 45°, is the local earth radius, γ 45o = 9.80665 m / s² -2 ,

[0063] S102. After obtaining the geopotential height at the specified location, determine the corresponding interpolation method and extrapolation method by comparing the relative position relationship between the interpolation point and the geopotential height of the isobaric surface of the numerical meteorological model, including the following three cases:

[0064] (1) When the geopotential height of the interpolation point is inside the isobaric surface, as shown in Figure 2 , first determine 8 grid points m1 - m8 adjacent to the interpolation point W in the numerical meteorological model through the longitude, latitude and geopotential height of the interpolation point W, and interpolate towards the geopotential height of the interpolation point for 4 groups of two vertically opposite grid points with the same longitude and latitude (m1 and m5, m2 and m6, m3 and m7, m4 and m8) respectively:

[0065] In the vertical direction, temperature, specific humidity and wind speed are interpolated linearly, and pressure is interpolated using an exponential model;

[0066] Linear interpolation: The m_value in the formula j represents the interpolation result of the meteorological data (temperature, specific humidity and wind speed) in the vertical direction between the upper and lower grid points. m_value g_j and m_value g_j+4 are the numerical values of the meteorological data of the upper and lower grid points respectively, h g_j is the distance from the geopotential height of the interpolation point to the upper grid point, h g_j+4 is the distance from the geopotential height of the interpolation point to the lower grid point. Here, j = 1, 2, 3, 4, as in Figure 2 the subscript number of the grid point;

[0067] Exponential model interpolation: The P in the formula j is the interpolation result of the pressure in the vertical direction, P c_j and Pc_j+4 are the barometric pressure values obtained by extrapolating the lower grid points and upper grid points through the exponential extrapolation formula, taking P c_j as an example, and its extrapolation formula is: T v_g_j = T g_j (1 + 0.6077q g_j ), where P g_j is the barometric pressure value at the lower grid point, T g_j is the temperature value at the lower grid point, q g_j is the specific humidity value at the lower grid point, g0 is the gravitational acceleration constant, and R d is the dry air gas constant, R d = 287.054 J / K·kg; the calculation method of the extrapolated value P c_j+4 of the upper grid point barometric pressure is the same;

[0068] Horizontally, all meteorological data (temperature, specific humidity, wind speed, barometric pressure) are interpolated bilinearly; after obtaining the meteorological data difference results m_value(w 1,5 ), m_value(w 2,6 ), m_value(w 3,7 ), m_value(w 4,8 ) of the four interpolation points w 1,5 , w 2,6 , w 3,7 , w 4,8 at the same geopotential height as the interpolation point W through linear interpolation and exponential model interpolation, the final difference result m_value(W) of the interpolation point W is obtained using bilinear interpolation:

[0069] where is the latitude of the interpolation point, λ c is the longitude of the interpolation point, λ 2,6 is the longitude corresponding to the grid points m2 and m6, and λ 4,8 is the longitude corresponding to the grid points m4 and m8, is the latitude corresponding to the grid points m1 and m5, is the latitude corresponding to the grid points m3 and m7;

[0070] (2) When the interpolated point geopotential height is below the lowest isobaric surface, the longitude, latitude, and geopotential height of the interpolated point W are also used to determine the eight grid points m1 to m8 adjacent to the interpolated point W. For the four grid points of the upper isobaric surface among the eight grid points, the bilinear interpolation method is used in the horizontal direction to obtain the meteorological parameters at the corresponding position of the interpolated point on the isobaric surface. In the vertical direction, the following methods are used to extrapolate and calculate the meteorological parameters at the corresponding position of the interpolated point on the isobaric surface:

[0071] The air pressure uses an exponential model;

[0072] The temperature varies linearly with a constant temperature gradient of -0.0065 °C / m: T j = T v_j -0.0065h g_j , T v_j is the temperature obtained by horizontal bilinear interpolation on the isobaric surface;

[0073] The specific humidity remains constant below the lowest isobaric surface: q j = q v_j , q v_j is the specific humidity obtained by horizontal bilinear interpolation on the isobaric surface;

[0074] The wind speed uses power-law interpolation: u v_j and v v_j are the eastward component of the horizontal wind and the northward component of the horizontal wind obtained by horizontal bilinear interpolation on the isobaric surface respectively, h j is the altitude of the interpolated point, h v_j is the altitude of the grid point, m is the wind shear coefficient, which can be obtained by looking up the table;

[0075] (3) When the interpolated point geopotential height is above the highest isobaric surface, the longitude, latitude, and geopotential height of the interpolated point W are used to determine the eight grid points m1 to m8 adjacent to the interpolated point W. For the four grid points of the lower isobaric surface among the eight grid points, the bilinear interpolation method is used in the horizontal direction to obtain the meteorological parameters at the corresponding position of the interpolated point on the isobaric surface. In the vertical direction, the following methods are used to extrapolate and calculate the meteorological parameters at the corresponding position of the interpolated point on the isobaric surface:

[0076] The specific humidity is set to zero;

[0077] The air pressure is extrapolated using an exponential model based on horizontal bilinear interpolation;

[0078] The temperature and wind speed components are obtained using the CIRA86 International Reference Atmosphere Model.

[0079] S2. Based on the meteorological data obtained through the numerical meteorological model, calculate the electromagnetic wave propagation path in the main axis direction of the transmitting antenna and the electromagnetic wave propagation path in the main axis direction of the receiving antenna;

[0080] In this embodiment, the ray tracing method based on the geometric optical ray theory is adopted, and the electromagnetic wave propagation path in the main axis direction of the transmitting antenna and the electromagnetic wave propagation path in the main axis direction of the receiving antenna are calculated based on the meteorological data obtained in step S1. The calculation processes of the electromagnetic wave propagation path in the main axis direction of the transmitting antenna and the electromagnetic wave propagation path in the main axis direction of the receiving antenna are the same. The following takes the electromagnetic wave propagation path in the main axis direction of the transmitting antenna as an example for description:

[0081] S201. Establish a local rectangular coordinate system for the ray tracing method. As Figure 3 shown, taking the position where the transmitting antenna is located as the coordinate origin, a local rectangular coordinate system xyz is established based on the right-hand screw rule. The x-axis of the coordinate system is the projection of the initial direction of the ray on the tangent plane of the earth's surface, the z-axis is the direction vector from the center of the earth to the position where the transmitter is located, and α0 and β0 are the azimuth angle and elevation angle respectively pointed by the main axis of the transmitting antenna;

[0082] In the local rectangular coordinate system xyz, the electromagnetic wave propagation path satisfies the linear differential equation system In the formula, n is the atmospheric refraction index, n = 10 -6 N + 1, N is the atmospheric refractive index, P, T, and q are the air pressure, temperature, and specific humidity obtained at the specified position from the numerical meteorological model through step S1. k1, k2, and k3 are the atmospheric refraction coefficients, k1 = 77.6890×10 -2 k / hPa, k2 = 71.2952×10 -2 k / hPa, k3 = 375.463×10 3 k 2 / hPa, R d is the dry air gas constant, R v is the water vapor gas constant, R v = 461.525 J / (kg·K), and e is the water vapor pressure;

[0083] Under the constraints of the above linear differential equation system under the initial boundary conditions and the tropopause height (86 km) limit, based on the meteorological data of the numerical meteorological model, the propagation path of the electromagnetic wave from the transmitting antenna to the tropopause can be obtained by iterative numerical solution, and thus the electromagnetic wave propagation path in the main axis direction of the transmitting antenna is obtained. For the receiving antenna, the same method can be used with the receiving antenna as the starting point for calculation.

[0084] S3, splitting the transmitting beam in the horizontal and vertical directions and rotating the path to obtain the ray path of each transmitting sub-beam, splitting the receiving beam in the horizontal direction, rotating the path and aligning the rays to obtain the ray path of the receiving sub-beam within the effective range, and intercepting the path from the transmitting antenna to the scattering point and then to the receiving antenna from the obtained ray path of the transmitting sub-beam and the ray path of the receiving sub-beam within the effective range, which is the tropospheric scattering propagation path of the scattering sub-beam;

[0085] The scattering beam includes a transmitting beam and a receiving beam. Tropospheric scattering is a typical multipath transmission channel. A single antenna main axis ray cannot reflect the real situation of the channel. The present invention adopts a beam splitting method to decompose the scattering process into a multi-sub-beam propagation process, and determines the tropospheric scattering transmission loss and propagation delay by calculating the receiving power and propagation time of each transmitting sub-beam.

[0086] Without considering atmospheric refraction, assuming that the transmitting and receiving antennas emit a conical beam (the vertical and horizontal widths of the beam are equal), the 3dB beam width of the antenna is the boundary of the conical beam; Figure 4 As shown, the rectangular cross-section partitioning method is used for the transmit beam. The transmit beam is partitioned into K equal parts along the horizontal axis and the vertical axis of the transmit beam cross section (beam propagation direction cross section) at equal angle intervals dw, and the transmit beam is decomposed into a plurality of closely arranged and non-overlapping regular tetrahedral transmit sub-beams. The center of each regular tetrahedral sub-beam corresponds to a ray. The elevation angle β corresponding to the ray of the transmit sub-beam is 0_ij and azimuth α 0_ij for α0 and β0 are the azimuth and elevation angles of the main axis of the transmitting antenna, respectively. i,j are the numbers of the horizontal and vertical divisions of the beam, i,j = -K / 2, -K / 2+1, ..., K / 2-1, K / 2;

[0087] Since the scattering beam width is usually narrow, the atmospheric refractive index changes of each transmitting / receiving sub-beam along the propagation direction within the beam width are basically the same, and it can be considered that the bending degree of each transmitting / receiving sub-beam is also the same. Therefore, the electromagnetic wave propagation path in the main axis direction of the transmitting / receiving antenna is rotated to obtain the ray path of each transmitting / receiving sub-beam; for the transmitting sub-beam, the electromagnetic wave propagation path in the main axis direction of the transmitting antenna obtained in step S2 is obtained. transmit_ray_main The local rectangular coordinate system xyz is converted to the spherical coordinate system BLH. The method of coordinate system conversion is a well-known technology in the art and is not the innovation of the present invention. It will not be described in detail here. Then, the initial direction of the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system BLH is rotated to the initial direction of each transmitting sub-beam ray to obtain the ray path of each transmitting sub-beam: the ray path l of the transmitting sub-beam with the serial number subscript ij transmit_ray_ij = l transmit_ray_main(B, L, H) + (dw·i, dw·j, 0), l transmit_ray_main (B, L, H) is the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system BLH;

[0088] According to the theory of turbulent incoherent scattering, the energy of the receiving beam comes from the secondary radiation of the dipoles excited by the transmitting beam in the common scatterers. Therefore, each transmitting sub - beam will excite secondary radiation in the common overlapping intersection area with the receiving antenna beam, that is, each transmitting sub - beam corresponds to multiple receiving sub - beams. For the receiving beam, its subdivision method is different from that of the transmitting beam. The elevation angles of the receiving sub - beams are evenly distributed at equal angular intervals within the width of the receiving beam, and the azimuth angles are adjusted using the ray alignment algorithm. The receiving sub - beams are specifically constructed through three processes: transverse subdivision, path rotation, and ray alignment:

[0089] Transverse subdivision: As Figure 5 shown, in the cross - section of the receiving beam, the receiving beam is transversely divided into K equal parts at an equal angular interval dw, obtaining several receiving sub - beams with the same azimuth angle but different elevation angles. The center of each receiving sub - beam also corresponds to a ray, and the azimuth angle corresponding to the ray of the receiving sub - beam is α0', and the elevation angle β 0_ij ' is β 0_ij ' = β0' + dw·i, where β0' is the elevation angle pointed by the main axis of the receiving antenna, and i is the serial number of the transverse subdivision of the beam, i = -K / 2, -K / 2 + 1, …, K / 2 - 1, K / 2;

[0090] Path rotation: The ray path of the receiving sub - beam is obtained by rotating the electromagnetic wave propagation path in the main axis direction of the receiving antenna, that is, the electromagnetic wave propagation path l r_ray_main in the local rectangular coordinate system xyz is converted to the spherical coordinate system BLH, and then the initial direction of the electromagnetic wave propagation path in the main axis direction of the receiving antenna in the spherical coordinate system BLH is rotated to the initial direction of each receiving sub - beam to obtain the ray paths of each receiving sub - beam: The ray path l r_ray_i of the receiving sub - beam with the serial number subscript ij is l r_ray_main =(B, L, H)+(dw·i, 0, 0), l r_ray_main (B, L, H) is the electromagnetic wave propagation path in the main axis direction of the receiving antenna in the spherical coordinate system BLH;

[0091] Ray alignment: Usually, the rays of the transmitting sub - beam and the corresponding receiving sub - beam obtained through path rotation are skew lines that do not intersect in three - dimensional space, and the azimuth angle of the receiving sub - beam needs to be adjusted to achieve the alignment of the transmitting and receiving rays. The azimuth angle of the receiving sub - beam is related to the transmitting sub - beam that excites the receiving sub - beam. As Figure 6As shown, for each transmitting sub-beam, the ray alignment method is used to determine the azimuth angles of the corresponding receiving sub-beams. Specifically, for each transmitting sub-beam, the minimum distance point p between the receiving sub-beam ray and the transmitting sub-beam ray is calculated. t_min p r_min and the minimum distance d between the minimum distance points min are obtained. The included angle Δα between p t_min and p r_min with respect to the transmitting point of the receiving antenna is calculated. For each receiving sub-beam ray, the azimuth angle of the receiving sub-beam ray is rotated according to the magnitude of the included angle Δα until the minimum distance d min is less than the set threshold (the threshold is set to the meter level, and the set threshold in this embodiment is 1 m). Calculate the included angle between the initial azimuth angle of the receiving sub-beam ray and the azimuth angle of the main axis of the receiving antenna at this time (when d min is less than the set threshold). Discard the receiving sub-beam rays whose azimuth angle included angle exceeds the beam width, and obtain the ray paths of the receiving sub-beams within the effective range. Intercept from the obtained ray paths of the transmitting sub-beams and the ray paths of the receiving sub-beams within the effective range the path from the transmitting point (transmitting antenna) to the scattering point (when d min is less than the set threshold, it is determined that p t_min and p r_min are the scattering points) and then to the receiving point (receiving antenna), which is the tropospheric scattering propagation path of the scattering sub-beam. The transmitting ray and the receiving ray are skew curves in space and do not intersect at a point. Therefore, there are two minimum distance points between them. One is the minimum distance point p t_min on the transmitting ray, and the other is the minimum distance point p r_min on the receiving ray. d min is the distance between these two points, and Δα is the included angle between the connecting lines of point p t_min and point p r_min and the transmitting point of the receiving antenna respectively. After the receiving ray is rotated by Δα, the value of d min is less than the set threshold, and its length can be ignored. p t_min and p r_min almost coincide in space position and can be determined as the scattering point.

[0092] S4. Calculate the volume of the common scatterer of each sub-beam; the common scatterer of the sub-beam is the intersection part of a transmitting sub-beam and a receiving sub-beam, that is, the common part of the transmitting sub-beam and the receiving sub-beam. Through step S3, the ray paths of the transmitting sub-beam and the receiving sub-beam are obtained. The transmitting sub-beam is a regular square pyramid, and the included angle between two opposite side faces of the regular square pyramid is known as dw. The boundary included angle of the elevation direction of the receiving sub-beam is known as dw. The common scatterer of the scattering sub-beam is approximately the hexahedron shown in Figure 7 , and its volume can be calculated by mathematical geometric methods.

[0093] Electromagnetic waves are affected by the change of refractive index in the atmosphere, and the propagation path bends towards the Earth, and its effect is equivalent to raising the transmitting and receiving antennas. To accurately calculate the volume of the common scatterer, in this embodiment, the volume of the common scatterer of the sub-beam is calculated by setting virtual transmitting points and virtual receiving points; as Figure 8 shown, based on the tropospheric scatter propagation path of the scattered sub-beam obtained through step S3, in the plane formed by the scatter point S (i.e., the intersection point of the transmitting sub-beam ray and the receiving sub-beam ray), the transmitting point T, and the receiving point R, starting from the scatter point S, sub-beam ray path tangent vectors are respectively made in the directions of the transmitting point T and the receiving point R, and the lengths of the tangent vectors are respectively set as the actual ray path lengths from the transmitting point T to the scatter point S and from the scatter point S to the receiving point R, and the end points T' and R' of the tangent vectors are the virtual transmitting point and the virtual receiving point;

[0094] In the plane formed by the virtual transmitting point T', the scatter point S, and the virtual receiving point R', with the virtual transmitting point T' as the origin o', the direction from the virtual transmitting point T' to the scatter point S as the x'-axis, and the direction from the center of the Earth to the virtual transmitting point T' as the y'-axis, a three-dimensional rectangular coordinate system x'y'z' is constructed based on the right-hand screw rule. In this three-dimensional rectangular coordinate system x'y'z', T'S and SR' are the virtual straight propagation paths of the electromagnetic waves. T'S and SR' are respectively rotated by 0.5dw and -0.5dw around the origin o' and the virtual receiving point R' in the x'o'y' plane, that is, the upper and lower boundaries of the transmitting sub-beam and the receiving sub-beam are obtained, and further the coordinates of the intersection points P1 - P4 of the upper and lower boundaries of the transmitting sub-beam and the receiving sub-beam in x'y'z' are obtained; the beam angle of the transmitting sub-beam is dw, and the coordinates of each vertex of the common scatterer (hexahedron) of the sub-beam can be obtained through the coordinates of P1 - P4 in x'y'z'; as Figure 9 shown, the common scatterer of the sub-beam is divided into trapezoidal bodies with a length of 1 km along the x'-axis direction (the propagation direction of the transmitting beam), the volume of each trapezoidal body is calculated, and the sum of the volumes of all trapezoidal bodies is the volume of the common scatterer of the sub-beam.

[0095] S5. Calculate the received power and propagation time of the scattered sub-beam;

[0096] In this embodiment, the bistatic radar equation is used to calculate the received power of the scattered sub-beam; first, the received power generated by the secondary excitation of each trapezoidal body is calculated respectively. The received power of the l-th trapezoidal body in the scattered sub-beam is: where P t_ij is the transmitting power of the scattered sub-beam with the serial number subscript ij, G t , G r are the transmitting antenna gain and the receiving antenna gain respectively, gt , g r are respectively the directivity function of the transmitting antenna and the directivity function of the receiving antenna, and σ ij_l is the scattering cross-section of the l-th frustum, ξ is the electromagnetic wave wavelength, and dV ij_l is the volume of the l-th frustum, R ij_l is the distance from the center point of the l-th frustum to the receiving antenna, and S ij_l is the distance from the center point of the l-th frustum to the transmitting antenna;

[0097] The gains and directivity functions of the transmitting and receiving antennas, the electromagnetic wave wavelength, etc. are independent of the meteorological environment and can be directly obtained through the basic parameters of the scattering link; the volume of the frustum and the distances from the center point of the frustum to the transmitting and receiving antennas are indirectly related to the meteorological environment and are obtained through the sub-beam ray paths obtained based on the numerical meteorological model; the scattering cross-section is closely related to the meteorological environment, and the variables M, u, v, and dT / dh related to the meteorological environment are all obtained through the numerical meteorological model; more specifically, the gains of the transmitting antenna and the receiving antenna can be calculated according to the antenna diameter D and the electromagnetic wave wavelength ξ, that is, the antenna gain G = 10lg(4.5(D / ξ) 2 );

[0098] The directivity function of the antenna is Gaussian: where ψ1 and ψ2 are respectively the 3dB widths of the transmitting sub-beam and the receiving sub-beam, β1 and β2 are respectively the elevation angles of the transmitting antenna and the receiving antenna, α1 and α2 are respectively the azimuth angles of the transmitting antenna and the receiving antenna, and β 10 , β 20 are respectively the elevation angles of the main axes of the transmitting antenna and the receiving antenna, and α 10 , α 20 are respectively the azimuth angles of the main axes of the transmitting antenna and the receiving antenna;

[0099] The scattering cross-section σ of the frustum ij_l is calculated according to the Kolmogorov theory and the Kolmogorov-Obukhov law, and σ ij_l = 2πk 4 cos(Θ ij ) 2 Φ(k), where k is the wave number, k = 2π / ξ, Θ ij is the scattering angle determined by the tropospheric scattering propagation path of the scattering sub-beam, and Φ(k) is the spatial spectrum function, is the refractive index structure constant, M is the refractive index gradient, L0 is the outer scale of turbulence, and can be calculated using the HMNSP99 outer scale model dT / dh is the temperature gradient, u and v are the eastward component and northward component of the horizontal wind respectively, h is the altitude, and M, u, v, and dT / dh are all obtained from the numerical weather model through step S1;

[0100] The received power P of the scattered sub-beam with the serial number subscript ij r_ij is the sum of the received powers of all trapezoids in the scattered sub-beam, P r_ij = ∑P r_ij_l ;

[0101] The propagation time τ of the scattered sub-beam with the serial number subscript ij ij is calculated by the refractive index integration method: S in the formula t is the ray path from the emission point to the scattering point, S r is the ray path from the scattering point to the receiving point, n is the atmospheric refractive index on the ray path, which can be obtained from the numerical weather model, c is the speed of light, and ds is the curve element along the ray path.

[0102] S6. Calculate the tropospheric scatter transmission loss and propagation delay;

[0103] The tropospheric scatter transmission loss is: where P r_ij is the received power of the scattered sub-beam with the serial number subscript ij, P t is the transmitted power of the transmitting antenna;

[0104] The propagation delay is represented by the delay power spectrum. According to step S5, the propagation time of each scattered sub-beam and the received power of each scattered sub-beam are obtained. The received powers of all scattered sub-beams are normalized, and the propagation times of all scattered sub-beams are arranged in ascending order to obtain the normalized delay power spectrum of the tropospheric scatter link, which is the tropospheric scatter propagation delay.

[0105] Atmospheric attenuation will affect the received power. Preferably, in this embodiment, based on the transmission loss L of the sub-beam ij the received power of the scattered sub-beam is corrected, and the received power of the corrected scattered sub-beam is: Then the tropospheric scatter transmission loss is: When calculating the tropospheric scatter propagation delay, the received power of each sub-beam is also preferably the received power of the corrected scattered sub-beam.

[0106] The transmission loss of the scattered sub-beam with the serial number subscript ij P r_ij is the received power of the scattered sub-beam with the serial number subscript ij, P t_ij is the transmitted power of the scattered sub-beam with the serial number subscript ij, is the atmospheric absorption loss of the scattered sub-beam with serial number subscript ij; considering that when electromagnetic waves propagate in the atmosphere, they are affected by dry air and water vapor in the atmosphere and generate atmospheric absorption loss, the atmospheric absorption loss of the scattered sub-beam is calculated by the line-by-line summation method of path integration in Recommendation ITU-R P.676-12: where ij_length is the segmented length of the tropospheric scatter propagation path of the scattered sub-beam (i.e., the length of the path array), a ij_b is the length of the b-th path segment, γ ij_b is the atmospheric attenuation ratio of the b-th path segment, γ ij_b is related to the temperature, atmospheric pressure, and specific humidity on the propagation path. The corresponding meteorological data is obtained through step S1. The specific calculation method refers to Recommendation ITU-R P.676-12 and will not be elaborated here.

[0107] In order to verify the calculation results of the transmission loss and transmission delay of the method of the present invention, the method of the present invention and the existing method are respectively used to calculate the transmission loss and transmission delay.

[0108] I. Calculation and analysis of transmission loss

[0109] The method of the present invention and two methods of ITU-R P.617-2 and ITU-R P.617-5 are used for comparative calculation of transmission loss. The numerical weather model used in the comparative calculation of the present invention is the ERA5 pressure stratification product provided by the European Centre for Medium Range Weather Forecasts (ECMWF). The ERA5 has a time resolution of 1 h, a spatial horizontal resolution of 0.25°, and is divided into 37 layers in the vertical direction according to isobaric surfaces from 1 hPa to 1000 hPa. The corresponding data is stored in the grid model.

[0110] The 10 tropospheric scatter link data used for comparison are selected from the CCIR Report OT / TRER 16 "Measurement and Prediction of Long-Term Distribution of Tropospheric Transmission Loss" published in 1986. The parameters of the 10 tropospheric scatter communication links are shown in Table 1 and Table 2 respectively.

[0111] Table 1

[0112]

[0113]

[0114] Table 2

[0115]

[0116] The ITU method uses the horizon angle instead of the actual antenna elevation angle and does not reflect the impact of beamwidth on transmission loss. The method of the present invention calculates the loss using the actual antenna elevation angle and beamwidth. In link construction, troposcatter generally uses a low elevation angle to reduce loss, and the antenna is adjusted to the optimal elevation angle with the lowest loss during communication. However, the optimal elevation angle varies with the meteorological environment. Therefore, the method of the present invention uses the ERA5 meteorological data on January 1, 2020 to iteratively calculate the antenna elevation angle with the minimum transmission loss, and the 3dB beamwidth of the antenna is obtained through the antenna gain. The antenna gains of the ten links are provided by the China Institute of Radio Wave Propagation, and the elevation angles and total gain magnitudes are shown in Table 1.

[0117] Table 3 gives the measured median of the transmission loss of the ten links relative to free space and the transmission loss calculation results of the three methods. The meteorological data of the present invention is the monthly average of ERA5 within the link observation time range.

[0118] Comparison of calculation results of three methods with observed values in Table 3

[0119]

[0120] It can be seen from the results in Table 3 that the calculation results of the method of the present invention are superior to the two ITU models in terms of mean error and root mean square error.

[0121] Transmission loss is closely related to the elevation angle. Figure 10 And Figure 11 respectively show the calculation results of the three methods for the two links LT2063 and LT2307 at different elevation angles. From Figure 10 And Figure 11 it can be seen that the transmission loss increases with the increase of the elevation angle. At low elevation angles, the differences between different methods are small; as the elevation angle increases, the differences tend to be more obvious. The results given by the ITU-R P.617-2 model at high elevation angles are much larger than those of the method of the present invention and ITU-R P.617-5. The difference between the method of the present invention and ITU-R P.617-5 also increases with the increase of the elevation angle, but the magnitude of the difference is smaller than the difference between ITU-R P.617-2 and the method of the present invention.

[0122] To demonstrate the hourly loss analysis ability of the method of the present invention based on meteorological data, the LT1441 and LT2305 links are selected to analyze the variation characteristics of troposcatter transmission loss based on the hourly ERA5 meteorological data. The positions of the two links are respectively as Figure 12 And Figure 13 shown, Figure 12 The LT1441 link located on the island of Newfoundland, Canada, is 277 kilometers long, and the scattered link part passes through the bay; Figure 13The LT2305 link is located between Tokyo and Fukushima in Japan, with a length of 226 kilometers. The scattering link passes through a piece of land. Figure 14 The hourly transmission losses calculated using the present invention on August 25th, 26th, and 27th, 2020 are given. As can be seen from the figure, the method of the present invention has successfully obtained the hourly change of transmission loss through the meteorological data of ERA5. Both links show a loss change trend of day and night alternation, with the loss increasing during the day and decreasing at night, but the rising and falling amplitudes during day and night are not the same. The LT1441 link is adjacent to the Atlantic Ocean and passes through a bay, with a complex and changeable meteorological environment. The transmission loss varies greatly at different times. The maximum hourly change can reach 6.3 dB, and the maximum day-night loss fluctuation is 20.4 dB; the LT2305 link passes through the land area, with relatively stable meteorological conditions, and the loss change is slight. The maximum hourly change is 2.4 dB, and the maximum day-night loss fluctuation is 5.3 dB. The loss change shows a periodic pattern.

[0123] II. Calculation and analysis of propagation delay

[0124] The propagation delay calculation ability of the method of the present invention is compared and analyzed with the Bello model based on the assumption of two-dimensional plane and straight-line light speed propagation of electromagnetic waves. To analyze and compare the performance of the two, six links with different distances and beam widths are set, and the parameters of the links are shown in Table 4. Figures 15a to 15f The normalized delay power spectra calculated by the two methods are given.

[0125] Table 4

[0126]

[0127]

[0128] Figure 5 Normalized delay power spectrum

[0129] From Figures 15a to 15f it can be seen that the delay power distributions of the method of the present invention and the Bello method are basically the same. When the beam width is the same, the farther the communication distance, the greater the time delay and the wider the delay spread; comparing Figure 15c and Figure 15d, with the same communication distance, the wider the beamwidth, the greater the delay spread. However, it should be noted that there is a certain offset in the delay spectrum distribution between the method of the present invention and the Bello method, that is, the calculated scattering propagation delays of the two are different. The propagation times calculated by the method of the present invention are all greater than those of the Bello method. This is because the Bello method derives the propagation delay based on the geometric configuration of the scattering link, assuming that the radio wave propagates in a straight line during the calculation process and does not consider the influence of the atmospheric environment on the radio wave propagation. While the method of the present invention calculates the bending and delay of the radio wave propagation path caused by the atmospheric environment in the generation of the scattering beam path and the delay statistics by introducing meteorological parameters. Therefore, the propagation delays are all greater than those of the Bello method. Obviously, the calculation results of the present invention are more accurate than those of the Bello method and are more suitable for situations that require accurate estimation of the transmission delay.

[0130] The change of the meteorological environment will affect the delay distribution of the scattering link. The ERA5 meteorological data of the 1441 link and the 2305 link on August 25th, 26th, and 27th, 2020 are selected to calculate the hourly change of the propagation delay. For the convenience of comparative analysis, except for the different positions and altitudes of the transceiver sites of the two links, the other link parameters are the same: frequency 4.09 GHz, transceiver antenna height 5.2 m, transceiver antenna horizontal angle 20 mrad, transceiver total gain 97 dB. The other link parameters refer to Table 1 and Table 2.

[0131] Figure 16a and Figure 16b The hourly normalized delay power spectrum heat maps of the two links calculated by the method of the present invention are given. It can be seen from the figure that the method of the present invention extracts the hourly delay change characteristics of the tropospheric scattering link. The changing meteorological environment at different times makes the propagation time and group delay of the scattering link different. The 1441 link passes through the Notre Dame Bay and is affected by the marine climate, with a complex and changeable meteorological environment. The scattering delay has weak regularity and strong randomness, and changes with the meteorological parameters. For example, Figure 16a from 19:00 to 23:00 on the 25th, the delay power spectrum shows a large mutation within a short time. This is mainly due to the sudden change of the temperature gradient in the vertical direction at the position of the scattering common body of some low-elevation sub-beams of the scattering link, which changes from about -5e-3 °C / m to about -3e-3 °C / m, thereby causing the received power of the sub-beams to decrease. And the propagation time of the low-elevation sub-beams is less than that of the high-elevation sub-beams, which is manifested as the peak of the maximum received power moving towards the region with a larger propagation delay in the normalized delay power spectrum. The 2305 link passes through land areas, and the meteorological environment is relatively stable. The scattering delay shows an obvious diurnal variation law. The propagation time during the day is significantly higher than that at night, and the variation range of the propagation time between day and night is about dozens of nanoseconds. Therefore, for time synchronization systems and passive detection systems that require accurate estimation of the propagation delay and group delay, the influence brought by the meteorological environment cannot be ignored.

[0132] The above are only the preferred embodiments of the present invention, and do not impose any form of limitation on the present invention. Although the present invention has been disclosed above with the preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some changes or modifications to equivalent embodiments with equivalent changes within the scope of the technical solution of the present invention by using the disclosed technical content above. However, as long as it does not depart from the content of the technical solution of the present invention, any simple modification, equivalent change and modification made to the above embodiments according to the technical essence of the present invention still fall within the scope of the technical solution of the present invention.

Claims

1. A method for calculating tropospheric scatter transmission loss and propagation delay, characterized in that, It includes the following steps: S1. Obtain meteorological data on the scattering propagation path through a numerical weather model; the numerical weather model stores meteorological data in a grid form, with a grid interval in the horizontal direction being the longitude and latitude resolution, and the vertical grid being an isobaric surface. Based on the numerical weather model, interpolation or extrapolation methods are used to obtain meteorological data at a specified location. The meteorological data includes temperature, specific humidity, air pressure, and wind speed. S2. Based on the meteorological data obtained through the numerical weather model, calculate the electromagnetic wave propagation path in the main axis direction of the transmitting antenna and the electromagnetic wave propagation path in the main axis direction of the receiving antenna. The ray tracing method is used to calculate the electromagnetic wave propagation path in the main axis direction of the transmitting antenna. The steps are as follows: S201. Take the location of the transmitting antenna as the coordinate origin, and establish a local rectangular coordinate system xyz based on the right-hand screw criterion. The x-axis of the coordinate system is the projection of the initial direction of the ray on the tangent plane of the earth's surface, the z-axis is the direction vector from the center of the earth to the location of the transmitter, and α0 and β0 are the azimuth angle and elevation angle respectively pointed by the main axis of the transmitting antenna. In this local rectangular coordinate system xyz, the electromagnetic wave propagation path satisfies a system of linear differential equations In the formula, n is the atmospheric refraction index and N is the atmospheric refractive index; The above linear differential equations are subject to the initial boundary conditions Under the constraints and the tropopause height limit, based on the meteorological data obtained from the numerical weather model, the electromagnetic wave propagation path in the main axis direction of the transmitting antenna can be obtained by iterative numerical solution. For the receiving antenna, the same steps are used with the receiving antenna as the starting point to obtain the electromagnetic wave propagation path in the main axis direction of the receiving antenna. S3. Perform horizontal and vertical dissections and path rotations on the transmitting beam to obtain the ray paths of each transmitting sub-beam, perform horizontal dissections, path rotations, and ray alignments on the receiving beam to obtain the ray paths of the receiving sub-beams within the effective range. Intercept the paths from the transmitting antenna to the scattering point and then to the receiving antenna from the obtained ray paths of the transmitting sub-beams and the ray paths of the receiving sub-beams within the effective range, which is the tropospheric scattering propagation path of the scattering sub-beam. S4. Calculate the volume of the common scatterer of each sub-beam. The common scatterer of the sub-beam is the intersection of a transmitting sub-beam and a receiving sub-beam. The volume of the common scatterer of the sub-beam is calculated through the following steps: In the plane formed by the scattering point S, the transmitting point T, and the receiving point R, starting from the scattering point S, make tangent vectors of the sub-beam ray paths in the directions of the transmitting point T and the receiving point R respectively. The lengths of the tangent vectors are set to the actual ray path lengths from the transmitting point T to the scattering point S and from the scattering point S to the receiving point R respectively. The end points T' and R' of the tangent vectors are the virtual transmitting point and the virtual receiving point. In the plane formed by the virtual emission point T', the scattering point S, and the virtual reception point R', with the virtual emission point T' as the origin o', the direction from the virtual emission point T' to the scattering point S as the x'-axis, and the direction from the center of the earth to the virtual emission point T' as the y'-axis, a three-dimensional rectangular coordinate system x'y'z' is constructed based on the right-hand screw rule. In this three-dimensional rectangular coordinate system x'y'z', T'S and SR' are the virtual straight-line propagation paths of the electromagnetic wave. T'S and SR' are respectively rotated by 0.5dw and -0.5dw around the origin o' and the virtual reception point R' in the x'o'y' plane to obtain the upper and lower boundaries of the transmitting sub-beam and the receiving sub-beam, and further obtain the coordinates of the intersection points P1 - P4 of the upper and lower boundaries of the transmitting sub-beam and the receiving sub-beam in x'y'z'; the coordinates of each vertex of the common scatterer of the sub-beam can be obtained through the coordinates of P1 - P4 in x'y'z'; the common scatterer of the sub-beam is subdivided into several trapezoidal bodies along the x'-axis direction, the volume of each trapezoidal body is calculated, and the sum of the volumes of all trapezoidal bodies is the volume of the common scatterer of the sub-beam; S5. Calculate the received power of the scattered sub-beam and the propagation time of the scattered sub-beam; S6. Calculate the tropospheric scattering transmission loss and propagation delay; The tropospheric scatter transmission loss is as follows: where P r_ij is the received power of the scatter sub-beam with the serial number subscript ij, and P t is the transmitted power of the transmitting antenna; Normalize the received power of all scattered sub-beams and arrange the propagation times of all scattered sub-beams in ascending order to obtain the normalized delay power spectrum of the tropospheric scattering link, which is the tropospheric scattering propagation delay.

2. The tropospheric scatter transmission loss and propagation delay calculation method according to claim 1, wherein: The steps of obtaining meteorological data at a specified location based on the numerical meteorological model in step S1 are specifically as follows: S101. Convert the geodetic height H at a specified position g to the geopotential height where is the latitude at the specified position, is the normal gravity on the earth ellipsoid surface, γ 45 is the standard value of the gravitational acceleration at latitude 45°, is the local earth radius; S102. Determine the corresponding interpolation method and extrapolation method by comparing the relative position relationship between the interpolation point and the geopotential height of the isobaric surface of the numerical meteorological model, including the following three cases: (1) When the geopotential height of the interpolation point is inside the isobaric surface, determine 8 grid points adjacent to the interpolation point in the numerical meteorological model through the longitude, latitude, and geopotential height of the interpolation point, and respectively interpolate the 4 groups of two relatively upper and lower grid points with the same longitude and latitude towards the geopotential height of the interpolation point: Vertically, linear interpolation is used for temperature, specific humidity, and wind speed, and exponential model interpolation is used for air pressure; Horizontally, bilinear interpolation is used for temperature, specific humidity, wind speed, and air pressure; (2) When the geopotential height of the interpolation point is below the bottom layer of the isobaric surface, determine 8 grid points adjacent to the interpolation point through the longitude, latitude, and geopotential height of the interpolation point. For the 4 grid points of the upper isobaric surface among the 8 grid points: Horizontally, use the bilinear interpolation method to obtain the meteorological parameters at the corresponding position of the interpolation point on the isobaric surface; Vertically, the following methods are respectively used to extrapolate and calculate the meteorological parameters at the corresponding positions of the interpolation points on the isobaric surface: For pressure, the exponential model is used, and for temperature T j = T v_j -0.0065h g_j , where T v_j is the temperature obtained by horizontal bilinear interpolation on the isobaric surface, and for specific humidity q j = q v_j , where q v_j is the specific humidity obtained by horizontal bilinear interpolation on the isobaric surface, and for wind speed, the power-law interpolation is used; (3) When the geopotential height of the interpolation point is above the top layer of the isobaric surface, determine 8 grid points adjacent to the interpolation point through the longitude, latitude, and geopotential height of the interpolation point. For the 4 grid points of the lower isobaric surface among the 8 grid points: Horizontally, use the bilinear interpolation method to obtain the meteorological parameters at the corresponding position of the interpolation point on the isobaric surface; Vertically, the following methods are respectively used to extrapolate and calculate the meteorological parameters at the corresponding position of the interpolation point on the isobaric surface: the specific humidity is zero; the air pressure is extrapolated using the exponential model based on horizontal bilinear interpolation; the temperature and wind speed components are obtained using the CIRA86 international reference atmosphere model.

3. The tropospheric scatter transmission loss and propagation delay calculation method according to claim 1, characterized in that: In the step S3, the emission beam is dissected by using a rectangular cross-section dissection method. The emission beam is dissected into K equal parts along the horizontal axis and the vertical axis of the cross-section of the emission beam at an equal angular interval dw, and the emission beam is decomposed into a plurality of non-overlapping emission sub-beams in the shape of regular quadrangular pyramids arranged closely; The electromagnetic wave propagation path \(l\) in the main axis direction of the transmitting antenna transmit_ray_main is converted from the local rectangular coordinate system \(xyz\) to the spherical coordinate system \(BLH\). Then, the initial direction of the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system \(BLH\) is rotated to the initial direction of each transmitting sub-beam ray. The ray path \(l\) of the transmitting sub-beam with the serial number subscript \(ij\) transmit_ray_ij = \(l\) transmit_ray_main (B, L, H)+(dw·i, dw·j, 0), where \(l\) transmit_ray_main (B, L, H) is the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system \(BLH\), and \(i,j = -K / 2, -K / 2 + 1, \cdots, K / 2 - 1, K / 2\).

4. The tropospheric scatter transmission loss and propagation delay calculation method according to claim 1, wherein: In the step S3, the received beam is dissected by a transverse dissection method. The received beam is transversely dissected into K equal parts at an equal angular interval dw in the cross-section of the received beam, and a number of received sub-beams with the same azimuth angle but different elevation angles are obtained; the electromagnetic wave propagation path l in the main axis direction of the transmitting antenna transmit_ray_main is converted from the local rectangular coordinate system xyz to the spherical coordinate system BLH, and then the initial direction of the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system BLH is rotated to the initial direction of each transmitting sub-beam ray. The ray path of the transmitting sub-beam with the serial number subscript ij l transmit_ray_ij = l transmit_ray_main (B, L, H) + (dw·i, dw·j, 0), l transmit_ray_main (B, L, H) is the electromagnetic wave propagation path in the main axis direction of the transmitting antenna in the spherical coordinate system BLH, i = -K / 2, -K / 2 + 1, …, K / 2 - 1, K / 2; for each transmitting sub-beam, calculate the minimum distance point between the receiving sub-beam ray and the transmitting sub-beam ray and the minimum distance between these two minimum distance points, obtain the included angle of these two minimum distance points with respect to the transmitting point of the receiving antenna, rotate the azimuth angle of the receiving sub-beam ray according to the size of the included angle until the minimum distance between the two minimum distance points is less than the set threshold, calculate the included angle between the initial azimuth angle of the receiving sub-beam ray and the azimuth angle of the main axis of the receiving antenna at this time, and discard the receiving sub-beam rays whose azimuth angle included angle exceeds the beam width to obtain the ray path of the receiving sub-beam within the effective range.

5. The tropospheric scatter transmission loss and propagation delay calculation method according to claim 1, wherein: In the step S5, the bistatic radar equation is used to calculate the received power of the scattered sub-beams. The steps are as follows: Calculate the received power of each frustum. The received power of the l-th frustum in the scattered sub-beam is: Among them, P t_ij is the transmission power of the scattered sub-beam with the sequence number subscript ij, G t , G r are the transmitting antenna gain and the receiving antenna gain respectively, g t , g r are the directivity functions of the transmitting antenna and the receiving antenna respectively, σ ij_l is the scattering cross section of the l-th frustum, ξ is the electromagnetic wave wavelength, dV ij_l is the volume of the l-th frustum, R ij_l is the distance from the center point of the l-th frustum to the receiving antenna, S ij_l is the distance from the center point of the l-th frustum to the transmitting antenna; The received power P of the scattered sub-beam with subscripts ij r_ij is the sum of the received powers of all the frustums in the scattered sub-beam, P r_ij = ∑P r_ij_l .

6. The tropospheric scatter transmission loss and propagation delay calculation method according to claim 1, characterized in that: In the step S5, the refractive index integration method is used to calculate the propagation time τ of the scattered sub-beam ij , where S in the formula t is the ray path from the emission point to the scattering point, S r is the ray path from the scattering point to the receiving point, n is the atmospheric refractive index on the ray path, c is the speed of light, and ds is the curve element.

7. The tropospheric scatter transmission loss and propagation delay calculation method according to claim 1, characterized in that: In step S6, first, based on the transmission loss L of the sub-beam ij correct the received power of the scattered sub-beam. The received power of the corrected scattered sub-beam is: In the formula, P r_ij is the received power of the scattered sub-beam with the serial number subscript ij, and P t_ij is the transmitted power of the scattered sub-beam with the serial number subscript ij. is the atmospheric absorption loss of the scattered sub-beam with the serial number subscript ij. Then, according to the received power P' of the corrected scattered sub-beam r_ij calculate the tropospheric scatter transmission loss. The tropospheric scatter transmission loss 8. The tropospheric scatter transmission loss and propagation delay calculation method according to claim 7, wherein: Calculate the tropospheric scatter transmission delay according to the received power P' of the corrected scatter sub-beams. r_ij ​

Images