Transmission matrix calculation method and system based on three-dimensional path discretization ionospheric model

Abstract

This invention belongs to the field of electromagnetic technology and discloses a method and system for calculating the transmission matrix based on a three-dimensional path discretization ionospheric model. The algorithm models the electromagnetic wave propagation process as a three-dimensional spatial curve from the transmission point to the receiving point, and constructs local grid cells along this curve, achieving a refined description of the non-uniform structure of the ionosphere along the propagation path. Based on this, by establishing a coordinate system aligned with the local propagation direction and introducing a coordinate rotation mechanism between the medium and field quantities, the traditional vertically layered transmission matrix is ​​extended to any three-dimensional path. This significantly improves the prediction accuracy of spaceborne VLF signals in long-distance, curved path propagation scenarios while maintaining the computational efficiency of the transmission matrix method. This invention systematically elaborates on the three-dimensional path modeling and grid discretization method, the establishment of the local coordinate system and the derivation of the rotation matrix, and finally presents a complete transmission matrix algorithm applicable to curved path propagation in a three-dimensional non-uniform ionosphere.

Description

Technical Field

[0001] This invention belongs to the field of electromagnetic technology, and in particular relates to a method and system for calculating the transmission matrix based on a three-dimensional path discretization ionospheric model. Background Technology

[0002] Very low frequency (VLF) electromagnetic waves refer to radio waves with operating frequencies between 3 kHz and 30 kHz. Due to their extremely long wavelengths and ability to penetrate media such as seawater, they have significant application value in long-distance communication, navigation, and military communications. VLF waves propagate in Earth-ionospheric waveguides, and their propagation characteristics are influenced by various factors such as changes in the medium parameters of the low and mid-elevation ionosphere, the curvature of the Earth, and waveguide boundary conditions. For a long time, researchers have proposed various theoretical and numerical models for VLF wave propagation, including methods based on full-wave theory, ray theory, mode theory, and layered transmission matrices. These methods provide a theoretical foundation and numerical tools for understanding and calculating the propagation behavior of VLF waves in different scenarios.

[0003] In existing technologies, a major propagation calculation method is the layered model based on the transmission matrix method. This model divides the ionospheric medium into a series of horizontal layers, assuming isotropic medium and uniform horizontal parameters within each layer. The transmission matrices of each layer are derived and cascaded to provide a holistic description of VLF electromagnetic wave propagation in a layered structure. Related literature has used this method to numerically simulate the field strength and phase of VLF waves within a propagation range of 500 km to 1000 km, and compared the results with experimental data. The simulations show that the traditional layered transmission matrix model achieves good consistency even when the ionospheric medium is not uniform in the vertical direction. This type of model is widely used in ionospheric propagation simulation and communication prediction due to its low computational cost and high algorithm stability.

[0004] However, the aforementioned existing technologies have significant flaws in their model assumptions. On the one hand, these layered models assume that the medium parameters within each altitude layer are uniform in the horizontal direction, neglecting the horizontal non-uniformity caused by changes in longitude and latitude, ionospheric disturbances, and the influence of the geomagnetic field. On the other hand, path curvature under the influence of the Earth's curvature, variations in propagation altitude, and three-dimensional geometric effects are not adequately considered. This leads to significant deviations in the numerical results of such models when the propagation distance exceeds a certain threshold (e.g., over 1000 km), particularly the core issue of the calculated propagation field strength and phase values ​​exceeding the preset tolerance (>8%). This deviation mainly stems from the neglect of horizontal non-uniformity and three-dimensional path effects. Traditional models do not incorporate any spatial variation descriptions or local coordinate adjustment mechanisms in their mathematical processing, resulting in insufficient response to factors such as changes in electron density and the influence of magnetic field direction during long-distance propagation, thus affecting prediction accuracy.

[0005] To address the aforementioned technical problems, while existing technologies such as the full-wave solution and ray theory can more accurately describe complex propagation scenarios, these methods either have extremely high computational complexity, making them difficult to apply in real time, or are still limited by discontinuous medium parameters and numerical instability in practical engineering. Therefore, the inadequacy of existing layered transmission matrix models in handling real-world ionospheric spatial inhomogeneities and ground curvature propagation problems has become a core technical pain point that urgently needs to be addressed. This deficiency directly leads to a significant reduction in the prediction accuracy of traditional methods in long-distance, non-uniform propagation scenarios. This technical problem is precisely the objective that the subsequent technical solutions of this invention must solve. Summary of the Invention

[0006] To address the problems existing in the prior art, this invention provides a method and system for calculating the transfer matrix based on a three-dimensional path discretization ionospheric model.

[0007] This invention is implemented as follows: a method for calculating the transfer matrix based on a three-dimensional path discretization ionospheric model includes: Step 1: Obtain the spatial coordinates of the transmitting and receiving points, establish a spatial rectangular coordinate system with the Earth's center as the origin, and construct a three-dimensional propagation path connecting the transmitting and receiving points; Step 2: Establish a height distribution function on the three-dimensional propagation path, assign a propagation height that changes continuously with the path parameters to each location point on the path, and form a three-dimensional propagation trajectory in the form of a spatial curve; Step 3: Discretize the three-dimensional propagation trajectory according to the axial scale adaptive principle to generate a spatial node sequence arranged in ascending order of arc length; Step four: Using the line connecting adjacent spatial nodes as the axis, establish local grid cells in each discrete segment, so that the medium parameters in each grid cell satisfy the condition of approximate uniformity, and use them as the basic unit for subsequent transmission matrix calculation.

[0008] Furthermore, the spatial rectangular coordinate system is a geocentric coordinate system with the Earth's center of mass as the origin. The transmitting and receiving points are obtained by converting geographical longitude, geographical latitude, and vertical altitude, and the conversion relationship includes: Calculate the radius of curvature of the ellipsoid, which is the Earth's semi-major axis divided by the square root of 1 minus the square of the first eccentricity multiplied by the square of the latitude sine. The X component of spatial coordinates is the sum of the radius of curvature and the height, multiplied by the cosine of latitude, and then multiplied by the cosine of longitude. The Y component of the spatial coordinates is the sum of the radius of curvature and the height, multiplied by the cosine of the latitude, and then multiplied by the sine of the longitude. The Z component of the spatial coordinates is the sum of the radius of curvature multiplied by 1, the square of the first eccentricity, and the altitude, multiplied by the latitude sine.

[0009] Furthermore, the planar reference path of the three-dimensional propagation path is represented by a geocentric great circle path, and the position vector of any point on the path is represented by spherical linear interpolation, which is expressed as follows: The path point vector equals the first endpoint vector multiplied by the sine function 1 minus the path parameter multiplied by the geocentric angle divided by the geocentric angle sine, plus the second endpoint vector multiplied by the sine function path parameter multiplied by the geocentric angle divided by the geocentric angle sine, where the path parameter takes values ​​from 0 to 1.

[0010] Furthermore, the height distribution function adopts a unimodal continuous function form, the value of which is equal to the minimum height plus the maximum height minus the minimum height multiplied by the shape function. The shape function takes a value of 0 when the path parameter is 0 and 1, and a value of 1 when the path parameter is 0.5.

[0011] Furthermore, the shape function is a Gaussian function, which is expressed as the negative path parameter of the exponential function minus 0.5 squared divided by 2 multiplied by the square of the standard deviation, where the standard deviation is 0.2.

[0012] Furthermore, the axial dimension adaptive principle includes wavelength constraints and medium gradient constraints, with the axial length taken as the smaller of the two calculated results, wherein: The wavelength constraint length is equal to the guided wave wavelength divided by the number of sampling points per wavelength; The length of the medium gradient constraint is equal to the preset relative change threshold divided by the electron density gradient factor.

[0013] Another objective of this invention is to provide a method for calculating the local coordinate rotation transfer matrix for a three-dimensional path discretization ionospheric model, comprising the following steps: Step 1: For each discretized grid cell, construct a local rectangular coordinate system with the grid axis as the propagation direction. The local coordinate system includes a first axis along the propagation direction, a second axis located in the plane formed by the propagation direction and the geographic perpendicular direction and perpendicular to the first axis, and a third axis determined according to the right-hand rule. Step 2: Calculate the azimuth and tilt corrections for the propagation direction based on the electron density at the grid midpoint and the geomagnetic field parameters, and update the local propagation direction. Step 3: Construct a rotation matrix from the geographic coordinate system to the local coordinate system. The rotation matrix consists of two basic rotations: the first rotation is around the geographic vertical axis to rotate the azimuth angle, and the second rotation is around the updated horizontal axis to rotate the tilt angle. Step 4: Transform the relative permittivity tensor in the geographic coordinate system to the local coordinate system using the rotation matrix to obtain the local equivalent relative permittivity tensor. Step 5: Establish a hierarchical transmission matrix in the local coordinate system and perform matrix concatenation calculations in grid order.

[0014] Furthermore, the azimuth correction is related to the electron density horizontal gradient, refractive index, electron cyclotron frequency, and grid axial length. The correction is equal to the proportionality constant multiplied by the electron density gradient, divided by the refractive index, multiplied by the geomagnetic correction factor, and then multiplied by the axial length.

[0015] Furthermore, the relative permittivity tensor in the local coordinate system is represented by a rotation matrix multiplied by the geographic coordinate system relative permittivity tensor, and then multiplied by the transpose of the rotation matrix.

[0016] Another objective of this invention is to provide a transfer matrix calculation system based on a three-dimensional path discretization ionospheric model, comprising: The path modeling module is used to construct a three-dimensional propagation path based on the coordinates of the transmitter and receiver points, and generate a spatial propagation trajectory that includes height distribution. The path discretization module is used to adaptively discretize the spatial propagation trajectory according to wavelength constraints and medium gradient constraints to generate a grid cell sequence. The local coordinate construction module is used to establish a local Cartesian coordinate system consistent with the propagation direction within each grid cell; The medium transformation module is used to construct the rotation matrix from the geographic coordinate system to the local coordinate system and to complete the coordinate transformation of the relative permittivity tensor. The transmission matrix calculation module is used to establish the unit transmission matrix in the local coordinate system and perform sequential cascade operations to output the overall transmission characteristics from the transmitting point to the receiving point.

[0017] Based on the above technical solutions and the technical problems solved, please analyze the advantages and positive effects of the technical solution to be protected by this invention from the following aspects: First, addressing the technical problems existing in the prior art and the difficulty in solving them, this paper closely analyzes, in conjunction with the technical solution to be protected by this invention and the results and data obtained during the research and development process, how the technical solution of this invention solves the technical problems, and the inventive technical effects brought about by solving these problems. The specific description is as follows: This invention proposes a three-dimensional path discretization transmission matrix algorithm. The algorithm models the electromagnetic wave propagation process as a three-dimensional spatial curve from the transmission point to the receiving point, and constructs local grid cells along this curve, achieving a refined description of the non-uniform structure of the ionosphere along the propagation path. Based on this, by establishing a coordinate system aligned with the local propagation direction and introducing a coordinate rotation mechanism between the medium and field quantities, the traditional vertically layered transmission matrix is ​​extended to arbitrary three-dimensional paths. This significantly improves the prediction accuracy of spaceborne VLF signals in long-distance, curved path propagation scenarios while maintaining the computational efficiency of the transmission matrix method. This invention systematically elaborates on the three-dimensional path modeling and grid discretization methods, the establishment of the local coordinate system and the derivation of the rotation matrix, and finally presents a complete transmission matrix algorithm applicable to curved path propagation in a three-dimensional non-uniform ionosphere.

[0018] Secondly, as supporting evidence of the inventiveness of this invention, it is also reflected in the following important aspects: (1) The expected benefits and commercial value of the technical solution of this invention after transformation are as follows: Prediction accuracy is significantly improved: In the long-distance (3200 km) VLF propagation scenario from ground transmission to ground reception, the field strength error rate between the three-dimensional algorithm and the measured value is 7.4%, which is significantly better than the 25% of the two-dimensional algorithm, and the phase error rate is reduced by 17.6%. In the spaceborne reception scenario, the three-dimensional corrected model is closest to the measured value, while the uncorrected model has a maximum deviation of 9.4 dB in the main ionospheric region.

[0019] Eliminating numerical divergence: In high-latitude, strong geomagnetic regions, the traditional two-dimensional algorithm without rotation mechanism drops to 0 dB of cross-polarization isolation at 1300 km and exhibits numerical divergence, while the three-dimensional algorithm with local rotation mechanism consistently maintains above 20 dB within a range of 2000 km.

[0020] Wide applicability: The algorithm is adaptive to frequency changes (5 kHz, 12 kHz, 25 kHz), and the computation time increases linearly with the number of grids without exponential explosion or numerical divergence; in ultra-long-distance cross-hemispheric propagation of 7500 km, the average relative deviation from the IRI model is 7.6%.

[0021] The above technical advantages demonstrate that this algorithm has clear engineering transformation value and commercial prospects in VLF application fields such as long-distance communication, navigation, military communication, and spaceborne reception.

[0022] (2) The technical solution of this invention fills a technical gap in the industry both domestically and internationally: The three-dimensional path discretization transmission matrix algorithm models the electromagnetic wave propagation process as a three-dimensional spatial curve from the transmission point to the receiving point, and constructs local grid cells along the curve, thereby realizing a fine description of the non-uniform structure of the ionosphere along the propagation path.

[0023] Local coordinate rotation mechanism: By establishing a coordinate system aligned with the local propagation direction and introducing a coordinate rotation mechanism for the medium and field quantities, the traditional vertical layered transmission matrix is ​​extended to any three-dimensional path.

[0024] Adaptive Discretization Strategy: Path discretization is performed based on both wavelength constraints and medium gradient constraints. The grid is automatically densified in regions with drastic changes in medium parameters, eliminating numerical oscillations in fixed-step-size models under strong gradient conditions. (3) Whether the technical solution of the present invention solves the technical problem that people have long wanted to solve but have never been able to solve successfully: Long-distance ground waveguide propagation (3200 km): The three-dimensional algorithm successfully reproduced the field strength fluctuation characteristics (fluctuation amplitude ±9~13 dB) caused by multimode interferometry, with an average deviation of less than 1.2 dB from the measured values; the two-dimensional algorithm and the ITU-R model had some point errors of 3~6 dB. At the receiving point, the error rate of the three-dimensional algorithm was 7.4%, which was significantly better than the 25% of the two-dimensional algorithm.

[0025] Strong electron density gradient environment: The fixed step size model exhibits significant numerical oscillations in the region of intense gradient (150~500 km), with a maximum deviation of 4.1 dB; the adaptive discrete model yields smooth and stable results, with a deviation consistently less than 0.6 dB.

[0026] In high-latitude, strong geomagnetic regions: the traditional two-dimensional algorithm exhibits numerical divergence at 1300 km, while the three-dimensional algorithm achieves an average cross-polarization isolation of 29 dB across the entire field, a 23.2 dB improvement over the traditional algorithm. (4) Does the technical solution of the present invention overcome technical bias? Bias 1: It is widely believed in the industry that the transfer matrix method is only applicable to ionospheric models with horizontal stratification and uniform horizontal parameters. This invention extends the transfer matrix method to arbitrary three-dimensional curved paths through three-dimensional path discretization and local coordinate rotation mechanisms, proving that the transfer matrix method is also applicable to non-uniform, curved path scenarios.

[0027] Bias 2: Traditional view holds that high-frequency (e.g., 25 kHz) VLF signals are prone to numerical instability when computed on large-scale paths. This invention, through adaptive grid discretization (automatically shortening grid length at high frequencies) and a local coordinate rotation mechanism, maintains numerical stability at three frequencies: 5 kHz, 12 kHz, and 25 kHz. The computation time increases linearly with the number of grids, verifying the algorithm's universality in adapting to frequency variations.

[0028] Bias 3: The industry tends to believe that VLF propagation prediction must be uniformly processed in a global coordinate system. This invention demonstrates the effectiveness and numerical stability of the local coordinate system method in 3D curved path propagation by establishing a local coordinate system aligned with the propagation direction within each grid and introducing a field rotation matrix between adjacent grids to achieve state vector transformation. Attached Figure Description

[0029] Figure 1 This is a flowchart of the transfer matrix calculation method based on a three-dimensional path discretization ionospheric model provided in an embodiment of the present invention.

[0030] Figure 2 This is a block diagram of the transmission matrix system structure based on a three-dimensional path discretization ionospheric model provided in an embodiment of the present invention.

[0031] Figure 3 This is a schematic diagram of the ECEF coordinate system provided in an embodiment of the present invention.

[0032] Figure 4 This is a schematic diagram of three-dimensional path parameterization and discretization provided in an embodiment of the present invention.

[0033] Figure 5 This is a schematic diagram of a local coordinate system provided in an embodiment of the present invention. Figure 6 This is a comparison chart of the field strength calculation results and measured data provided by the three-dimensional algorithm, the traditional two-dimensional layered model, and the ITU-R P.372-11 standard model in the embodiments of the present invention.

[0034] Figure 7 This is a comparison curve of the three-dimensional algorithm, the two-dimensional algorithm, and the measured phase data in this scenario provided by the embodiments of the present invention.

[0035] Figure 8 This is a comparison curve of the amplitude of the three-dimensional corrected model and the uncorrected model provided in the embodiments of the present invention with the measured satellite data.

[0036] Figure 9 This is a comparison chart of the cross-polarization isolation of the two-dimensional and three-dimensional algorithms provided in the embodiments of the present invention.

[0037] Figure 10 This is a comparison curve between the calculation results and the measured values ​​of the adaptive discrete model and the fixed step size model provided in the embodiments of the present invention.

[0038] Figure 11 This is a comparison curve of the field strength variation with distance between the three-dimensional algorithm and the IRI model provided in this embodiment of the invention.

[0039] Figure 12 This is a comparison curve of the field strength changing with distance at different frequencies, provided by an embodiment of the present invention.

[0040] Figure 13 This is a schematic diagram of the relative calculation time of field distribution at different frequencies provided in the embodiments of the present invention. Detailed Implementation

[0041] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0042] like Figure 1 As shown, the method for calculating the transfer matrix based on a three-dimensional path discretization ionospheric model provided by this embodiment of the invention includes the following steps: S101, obtain the spatial coordinates of the transmitting point and the receiving point, establish a spatial rectangular coordinate system with the Earth's center as the origin, and construct a three-dimensional propagation path connecting the transmitting point and the receiving point; S102, establish a height distribution function on the three-dimensional propagation path, assign a propagation height that changes continuously with the path parameters to each position point on the path, and form a three-dimensional propagation trajectory in the form of a spatial curve; S103, along the three-dimensional propagation trajectory, the space nodes are discretized according to the axial scale adaptive principle to generate a spatial node sequence arranged in ascending arc length. S104 establishes local grid cells in each discrete segment using the line connecting adjacent spatial nodes as the axis, so that the medium parameters in each grid cell satisfy the condition of approximate uniformity, and serves as the basic unit for subsequent transmission matrix calculation.

[0043] The spatial rectangular coordinate system provided in this embodiment of the invention is a geocentric coordinate system with the Earth's center of mass as the origin. The transmitting point and receiving point are obtained by converting geographical longitude, geographical latitude, and vertical altitude. The conversion relationship includes: Calculate the radius of curvature of the ellipsoid, which is the Earth's semi-major axis divided by the square root of 1 minus the square of the first eccentricity multiplied by the square of the latitude sine. The X component of spatial coordinates is the sum of the radius of curvature and the height, multiplied by the cosine of latitude, and then multiplied by the cosine of longitude. The Y component of the spatial coordinates is the sum of the radius of curvature and the height, multiplied by the cosine of the latitude, and then multiplied by the sine of the longitude. The Z component of the spatial coordinates is the sum of the radius of curvature multiplied by 1, the square of the first eccentricity, and the altitude, multiplied by the latitude sine.

[0044] The planar reference path of the three-dimensional propagation path provided in this embodiment of the invention is represented by a geocentric great circle path, and the position vector of any point on the path is represented by spherical linear interpolation, which is expressed as follows: The path point vector equals the first endpoint vector multiplied by the sine function 1 minus the path parameter multiplied by the geocentric angle divided by the geocentric angle sine, plus the second endpoint vector multiplied by the sine function path parameter multiplied by the geocentric angle divided by the geocentric angle sine, where the path parameter takes values ​​from 0 to 1.

[0045] The height distribution function provided in this embodiment of the invention adopts a unimodal continuous function form, the value of which is equal to the minimum height plus the maximum height minus the minimum height multiplied by the shape function. The shape function takes a value of 0 when the path parameter is 0 and 1, and takes a value of 1 when the path parameter is 0.5.

[0046] The shape function provided in this embodiment of the invention is a Gaussian function, which is expressed as the negative path parameter of the exponential function minus 0.5 squared divided by 2 multiplied by the square of the standard deviation, where the standard deviation is 0.2.

[0047] The axial dimension adaptive principle provided in this embodiment of the invention includes wavelength constraint and medium gradient constraint, and the axial length is taken as the smaller value of the two calculated results, wherein: The wavelength constraint length is equal to the guided wave wavelength divided by the number of sampling points per wavelength; The length of the medium gradient constraint is equal to the preset relative change threshold divided by the electron density gradient factor.

[0048] This invention provides a method for calculating the local coordinate rotation transfer matrix of a three-dimensional path discretized ionospheric model, comprising the following steps: Step 1: For each discretized grid cell, construct a local rectangular coordinate system with the grid axis as the propagation direction. The local coordinate system includes a first axis along the propagation direction, a second axis located in the plane formed by the propagation direction and the geographic perpendicular direction and perpendicular to the first axis, and a third axis determined according to the right-hand rule. Step 2: Calculate the azimuth and tilt corrections for the propagation direction based on the electron density at the grid midpoint and the geomagnetic field parameters, and update the local propagation direction. Step 3: Construct a rotation matrix from the geographic coordinate system to the local coordinate system. The rotation matrix consists of two basic rotations: the first rotation is around the geographic vertical axis to rotate the azimuth angle, and the second rotation is around the updated horizontal axis to rotate the tilt angle. Step 4: Transform the relative permittivity tensor in the geographic coordinate system to the local coordinate system using the rotation matrix to obtain the local equivalent relative permittivity tensor. Step 5: Establish a hierarchical transmission matrix in the local coordinate system and perform matrix concatenation calculations in grid order.

[0049] The azimuth correction provided in this embodiment of the invention is related to the electron density horizontal gradient, refractive index, electron cyclotron frequency, and grid axial length. The correction is equal to the proportionality constant multiplied by the electron density gradient, divided by the refractive index, multiplied by the geomagnetic correction factor, and then multiplied by the axial length.

[0050] The representation of the relative permittivity tensor in the local coordinate system provided in this embodiment of the invention is a rotation matrix multiplied by the geographic coordinate system relative permittivity tensor, and then multiplied by the transpose of the rotation matrix.

[0051] like Figure 2 As shown in the figure, an embodiment of the present invention provides a transfer matrix calculation system based on a three-dimensional path discretization ionospheric model, comprising: The path modeling module is used to construct a three-dimensional propagation path based on the coordinates of the transmitter and receiver points, and generate a spatial propagation trajectory that includes height distribution. The path discretization module is used to adaptively discretize the spatial propagation trajectory according to wavelength constraints and medium gradient constraints to generate a grid cell sequence. The local coordinate construction module is used to establish a local Cartesian coordinate system consistent with the propagation direction within each grid cell; The medium transformation module is used to construct the rotation matrix from the geographic coordinate system to the local coordinate system and to complete the coordinate transformation of the relative permittivity tensor. The transmission matrix calculation module is used to establish the unit transmission matrix in the local coordinate system and perform sequential cascade operations to output the overall transmission characteristics from the transmitting point to the receiving point.

[0052] Specific implementation of the present invention: 1. This invention proposes a three-dimensional path discretization transmission matrix algorithm. This algorithm models the electromagnetic wave propagation process as a three-dimensional spatial curve from the transmission point to the receiving point, and constructs local grid cells along this curve, achieving a refined description of the non-uniform structure of the ionosphere along the propagation path. Based on this, by establishing a coordinate system aligned with the local propagation direction and introducing a coordinate rotation mechanism between the medium and field quantities, the traditional vertically layered transmission matrix is ​​extended to any three-dimensional path. This significantly improves the prediction accuracy of spaceborne VLF signals in long-distance, curved path propagation scenarios while maintaining the computational efficiency of the transmission matrix method. This invention systematically elaborates on the three-dimensional path modeling and grid discretization methods, the establishment of the local coordinate system and the derivation of the rotation matrix, and finally presents a complete transmission matrix algorithm applicable to curved path propagation in a three-dimensional non-uniform ionosphere.

[0053] 2. 3D path modeling and raster discretization To achieve accurate propagation modeling of very low frequency (VLF) waves in a three-dimensional non-uniform ionosphere, the continuous physical propagation path needs to be transformed into a three-dimensional grid system suitable for numerical calculation. This invention first establishes a transformation relationship from the geographic coordinate system to the Earth-Centered Earth-Fixed (ECEF) coordinate system, providing a unified framework for path geometry calculation and parameter lookup. Based on this, a three-dimensional spatial curve describing the true wave trajectory is constructed using a great circle path and height profile function, and then adaptively discretized into a series of grid cells. Subsequently, the geometric dimensions and equivalent electromagnetic parameters of the grid are determined, and finally, a local Cartesian coordinate system aligned with the propagation direction is established for each grid, thus laying the geometric and physical foundation for the implementation of the three-dimensional transmission matrix algorithm.

[0054] 2.1 Coordinate system establishment and great circle path description First, a geographic coordinate system is established. In ionospheric very low frequency propagation modeling, the geographic coordinate system serves as the fundamental framework for data input and physical parameter lookup, and its definition directly affects the accuracy of subsequent calculations. This invention employs a geodetic coordinate system with the geocenter as the reference. For any point P in space, its geographic center coordinates are... Unique identifier, among which Indicates geographical longitude. Indicates geographical latitude, It indicates vertical height, usually measured from sea level.

[0055] The actual propagation path of electromagnetic waves in a three-dimensional non-uniform ionosphere follows Fermat's principle and is a complex curve determined by the spatial distribution of the medium's refractive index. To perform three-dimensional transmission matrix calculations in a numerically tractable and physically interpretable manner, this invention selects a great circle path as the baseline path for spatial discretization.

[0056] Although geographic coordinate systems can directly describe the observation location, their spherical geometry necessitates complex spherical trigonometric formulas for calculating great circle paths. To address the Earth's curvature effect and achieve efficient geometric calculations in three-dimensional space, the Earth-centered Earth-fixed (ECEF) Cartesian coordinate system is introduced. This system uses the Earth's center of mass as its origin, mapping the Earth's surface into three-dimensional Euclidean space. This simplifies great circle paths to linear interpolation within a plane, reducing geometric calculations such as propagation direction and local coordinate system establishment to direct vector operations.

[0057] like Figure 3 The ECEF coordinate system (X,Y,Z) is defined as follows: with the Earth's center of mass as the origin, the Z-axis points towards the North Pole, the X-axis points towards the intersection of the Prime Meridian and the equatorial plane, and the Y-axis forms a right-handed orthogonal system with the X and Z axes.

[0058] ECEF coordinates can be obtained from geographic coordinates through the following transformation relationship: (1) in Let be the radius of curvature of the Earth's ellipsoid in that direction. The semi-major axis of the Earth's ellipsoid. It is the first eccentricity of the Earth's ellipsoid.

[0059] The ECEF coordinate system simplifies the calculation of straight-line distances, assuming the launch point... Receiving point The straight-line distance between the two points is... .

[0060] 2.2 3D Path Parameterization and Discretization Based on the established coordinate system, to achieve numerical processing of the three-dimensional propagation path, the continuous wave propagation trajectory needs to be parameterized into a spatial curve suitable for discrete calculation. This invention first describes the horizontal propagation trajectory based on a great circle path in the ECEF coordinate system, then introduces a height profile function to characterize the refraction and reflection behavior of the wave in the vertical direction, finally synthesizing a complete three-dimensional spatial path. This path will be further discretized into a series of nodes, serving as the geometric basis for subsequent grid cell partitioning, thus achieving a crucial transition from a continuous propagation model to a discrete grid system.

[0061] First, the three-dimensional path is discretized. The Earth's curvature and the launch point are considered. With receiving point The propagation path between them is a great circle path in three-dimensional space. In the ECEF coordinate system, this path can be accurately described by a parametric method. Let... It is the position vector of the launch point in the ECEF coordinate system. This is the position vector of the receiving point in the ECEF coordinate system. Let the geocentric angle be... Then the position vector of any point on the path This can be expressed using the Slerp formula: (2) parameter Corresponding launch point For the corresponding receiving point, this equation gives the planar great circle path determined by the Earth's center and two points. However, the propagation height of VLF waves in the ionosphere is not fixed but exhibits continuous variation, thus requiring the planar great circle path to be extended into a three-dimensional spatial curve. To modulate the path in the vertical direction and obtain the true three-dimensional electromagnetic wave propagation trajectory, a height profile function is introduced. : (3) in This is the lowest elevation of the path (usually the ground level). This represents the highest altitude along the path (maximum penetration height of the ionosphere). Let be a shape function, satisfying This invention employs a Gaussian contour function. : (4) in This function avoids numerical abrupt changes and smoothly describes the refraction-reflection propagation characteristics of VLF waves in ionospheric waveguides.

[0062] Finally, a 3D path is synthesized using a position vector function and a height contour function. : (5) This yields a three-dimensional spatial curve describing the VLF wave propagation trajectory, facilitating subsequent calculations. To apply the transmission matrix method, this curve needs to be... The path is divided into a series of sufficiently short segments such that the medium parameters are approximately uniform within each segment, and the wave propagation behavior can be described by a linear system. This division process is called path discretization, which sets a series of spatial location nodes along the path. and the corresponding arc length sequence ,satisfy: (6) in Indicates from the starting point to the... The cumulative arc length of each node, Corresponding launch point For the corresponding receiving point, the straight-line distance between adjacent nodes is the axial length of the subsequent grid cell, denoted as . Discretization effect as Figure 4 As shown.

[0063] This discretization process transforms a continuous path into a series of straight line segments connected end to end, providing a geometric basis for the construction of subsequent grid cells.

[0064] 2.3 Determination of Grid Geometric Parameters Based on the established node sequence, the geometric parameters of each grid cell need to be determined. In this invention, the grid cell is conceptually defined as a three-dimensional volume element, whose geometric parameters include axial length and cross-sectional dimensions. The axial length directly participates in the subsequent transmission matrix calculation, while the cross-sectional dimensions are mainly used to fully describe the spatial morphology of the grid, and its range is determined by the physical mechanism of VLF wave propagation.

[0065] (1) Axial length Sure First, determine the axial length of each grid cell. This length can be considered as an arc length increment within the grid. For the first... One grid, The wavelength constraint and the gradient constraint of the medium parameter together ensure a balance between computational efficiency and numerical accuracy.

[0066] Based on the spatial sampling theorem, each guided wave wavelength must contain enough sampling points to accurately describe the phase change. Guided wave wavelength The wavelength can be approximated by the vacuum wavelength. For the VLF frequency band (3-30 kHz), the guided wave wavelength is... To meet the requirements of numerical stability and phase accuracy, the minimum number of sampling points within each wavelength is set to be... Then the corresponding wavelength constraint length for: (7) Besides wavelength constraints, local variations in dielectric parameters also affect the selection of grating size. To satisfy the assumption of dielectric homogeneity within the grating, the relative rate of change of electron density must be controlled within a threshold value. Within. For adjacent path points and The electron density gradient factor is defined as : (8) Then the gradient constraint length for: (9) Combining the above two constraints, the grid axial length Take the minimum value: (10) Determined based on the above axial dimensions The sequence, along the path starting from the launch point, is sequentially truncated to a length of... Line segments, generating spatial location node sequences and the corresponding arc length parameter sequence .

[0067] These spatial locations These are the boundary nodes of the grid cells. The regions between adjacent nodes directly constitute the grid cells, ensuring consistency between the discretization density and the grid size, and simultaneously realizing path discretization and grid length construction. Furthermore, it achieves adaptive variable step-size discretization, especially in regions with drastic changes in medium parameters or significant path curvature. Automatically reduce resolution to improve local detail; in areas with uniform medium and straight paths, It can be appropriately increased to improve computational efficiency.

[0068] (2) Conceptual definition of cross-sectional range To conceptually fully describe the three-dimensional spatial extent of a grid cell, it is necessary to define the extent of its cross-section perpendicular to the propagation direction. This extent is determined by the physical characteristics of VLF wave propagation in the earth-ionospheric waveguide.

[0069] The physical scale of the cross-sectional area is based on the lateral confinement characteristics of VLF waveguide modes. Waveguide mode theory indicates that when a VLF wave propagates in an ionospheric waveguide, its energy decays exponentially in the plane perpendicular to the propagation direction, typically on the order of tens of kilometers. Measured data also show that the lateral coherence distance of VLF signals is usually within the tens of kilometers range. Therefore, the grid's extension in the plane perpendicular to the propagation direction should be comparable to this physical scale to conceptually encompass the main distribution of wave energy.

[0070] It is important to emphasize that, since the propagation path is a three-dimensional spatial curve, the local propagation direction of each grid cell is different, and therefore the orientation of the cross-sectional plane of each grid cell in the global geographic coordinate system is also different. This model assumes that within each grid cell, the medium parameters are uniformly distributed within the cross-sectional plane of that grid cell, and the lateral distribution pattern of the VLF wave is relatively stable. This assumption simplifies the propagation process of electromagnetic waves into a one-dimensional transmission problem along the grid axis.

[0071] Thus, by determining the axial length and cross-sectional range of the grid cell, the grid cell is conceptually defined as a spatial region with a finite extension in the vertical plane, centered on the propagation direction. Specifically, for the first... One grid cell (corresponding to a path segment) Its central axis is a line segment. ,length .

[0072] 2.4 Equivalent Electromagnetic Parameters and Local Coordinate System The previous section determined the geometric parameters of the grid (including the axial length and cross-sectional area definition). The next step is to determine the parameters of each... Equivalent electromagnetic parameters within the first grid cell. This invention employs an axis midpoint equivalence strategy for the first... Each grid is defined by its central axis midpoint. The electromagnetic parameters at that location are equivalent to the electromagnetic parameters of that grid cell. Key physical parameters of the ionosphere are obtained using the international standard model. Electron density. and collision frequency Geomagnetic field strength can be obtained through IRI model query. Geomagnetic tilt Magnetic declination It can be obtained through the World Geomagnetic Model (WMM).

[0073] After completing the construction of the three-dimensional grid cells, in preparation for the subsequent application of the hierarchical transmission matrix algorithm to the electromagnetic wave propagation problem, the next step is to establish a local Cartesian coordinate system within each grid cell that matches the propagation direction. The establishment of this coordinate system is independent of the specific shape and size of the grid cross-section. The origin is set at the grid starting point. , The axis along the grid axis, which is the actual three-dimensional propagation direction, represents the local propagation direction of the wave within this grid. The axis lies in the plane formed by the direction of propagation and the geographic perpendicular direction (zenith direction), and is perpendicular to... axis; Depend on and Determined according to the right-hand rule, such as Figure 5 As shown.

[0074] In the above local coordinate system The axial direction is initially determined by the geometric relationship of the great circle path, without considering the refraction effect caused by the horizontal inhomogeneity of the ionosphere. To more accurately describe the wave propagation path, this invention introduces a first-order direction correction based on ray optics theory. The core idea of ​​the horizontal refraction correction is that the horizontal gradient of electron density changes the direction of the local wave normal. Let the azimuth angle be the angle between the propagation direction and geographic north. The angle between the plane and the horizontal plane is the angle of inclination. For the first For each grid cell, the azimuth correction is calculated based on the ray equation and the refractive index expression under the quasi-longitudinal approximation. It can be represented as: (10) in, It is the basic proportionality constant. The refractive index; This is the geomagnetic field correction factor. It is the electron cyclotron frequency; The electron density at the midpoint of the lattice; This represents the axial length of the grid.

[0075] The corrected propagation direction azimuth angle for: (10) Let the azimuth be the angle between the (corrected) direction of propagation and geographic north. The angle between the plane and the horizontal plane is the angle of inclination. Because in a three-dimensional, non-uniform ionosphere, the propagation direction of electromagnetic waves, the direction of the geomagnetic field (principal axis of the dielectric tensor), and the geographic reference direction do not coincide, in order to correctly describe the wave propagation process and in the... Constructing the transmission matrix in the local coordinate system of each grid requires establishing a transformation chain from the geographic coordinate system to the local propagation coordinate system of each grid. Although the aforementioned ECEF coordinate system simplifies geometric calculations, solving for the electromagnetic field requires establishing a local Cartesian coordinate system aligned with the propagation direction within each grid, and transforming the anisotropic medium parameters defined in the geographic coordinate system to this local coordinate system.

[0076] Let the unit basis vector of the geographic coordinate system be... The unit basis vector of the local rectangular coordinate system is Rotation operations follow the rules of rigid body rotation in three-dimensional Euclidean space, and these two sets of basis rules... arrive The transformation can be achieved through two elementary rotations: (11) The following are two basic rotations. First, around the geography... Axis rotation azimuth angle , to make new The axis is aligned with the horizontal propagation direction. This rotation corresponds to a rotation matrix written in three-dimensional Cartesian coordinates as... ,in This is equivalent to the geography department. axis: (12) Subsequently, around Xin Axis rotation tilt angle , to make new The axis is aligned with the actual three-dimensional propagation direction; this rotation corresponds to the matrix. : (13) The required medium rotation matrix It is the product of two fundamental rotation matrices. (14) The relative permittivity of a point in the geographic coordinate system Relative permittivity obtained by transformation to local coordinate system for: (15) 3. Three-dimensional transfer matrix algorithm based on local coordinate system 3.1 Calculation of a single grid transmission matrix In Chapter 2, we established the Berreman-form equation of state for VLF wave propagation in an anisotropic horizontally layered ionosphere. Consider the... Layered medium, its system matrix If the (Belemann matrix) is approximately constant within the layer, then the electromagnetic field state vector within that layer is... Satisfying a system of linear differential equations with constant coefficients: (16) Where the state vector It includes the transverse component of the electromagnetic field. This is achieved by solving the matrix. The problem of eigenvalues, to obtain their eigenvalues. With the eigenvector matrix (4×4 matrix), the transfer matrix within this layer can be represented as: (17) in It is a diagonal matrix with internal elements of . , For the first The thickness of the layer.

[0077] In actual three-dimensional propagation, electromagnetic waves travel along spatial curves. Propagation. To extend the above-described horizontal layered transfer matrix algorithm to three-dimensional curved paths, based on the grid discretization method in Section 3.2, the complete three-dimensional path of the ionosphere is discretized into... There are 1 grid cells, and the arc length of each grid cell is 1. And satisfy the following assumptions: the electromagnetic parameters (electron density, collision frequency, geomagnetic field direction) of the medium within the grid are approximately uniform; the propagation direction within the grid is approximately a straight line; the grid size is much smaller than the characteristic length of the medium parameter variation and the waveguide wavelength.

[0078] Under these conditions, the propagation process of each grid cell in its local coordinate system can still be considered as one-dimensional propagation along the axis. Therefore, the coordinates in the original horizontally layered vertical propagation model can be... Replace with arc length coordinates along the grid axis The state equation can be naturally generalized to: (18) At this time, the system matrix For the first The system matrix of a grid in its local coordinate system is determined by the equivalent electromagnetic parameters within the grid and the local propagation direction. Since the medium within the grid is homogeneous and the propagation direction is fixed, its transmission matrix can be directly adopted from the eigenvalue decomposition form in Chapter 2: (19) At this point, we have obtained the transfer matrix of each grid in its local coordinate system. However, each grid cell exists in its local coordinate system. The z-axis defined in the code is not uniform along the axis of each individual grid. The local coordinate systems of adjacent grids have different orientations, therefore the transfer matrix between grids varies. Direct cascading is not possible; a coordinate rotation matrix must be introduced to transform the state vectors of different grids to the same coordinate system. Therefore, the field rotation matrix between the local coordinate systems of adjacent grids will be derived next to achieve the correct cascading of the coordinate transformation of the state vectors and the global transfer matrix.

[0079] 3.2 Derivation of the Field Rotation Matrix In a horizontally layered model, all layers share the same global coordinate system. However, in a three-dimensional curved path, adjacent grids have different local coordinate systems. When a wave moves from one grid to the next grid with a different orientation, the state vector needs to be rotated from the local coordinate system of the previous grid to the local coordinate system of the next grid. This satisfies the field rotation matrix between the local coordinate systems of adjacent grids. Should meet: (20) in For the first Raster output (in its own coordinate system). For the first Raster input (in) (In the coordinate system). In Section 3.2, the medium rotation matrix. Used to rotate the dielectric tensor from a geographic coordinate system to a local propagation coordinate system. In this invention, the field rotation matrix... It is used to rotate the state vector from the previous grid local coordinate system to the next grid local coordinate system. The two have different physical meanings but similar mathematical forms, both originating from three-dimensional rotation transformation.

[0080] From electromagnetic field vector Derivation: Through extraction using a 4×6 projection matrix and algebraic operations, the two rotation matrices of the four-dimensional state vector can be obtained. Let the... With the The azimuth and inclination angles of the grid-based regional coordinate system relative to the geographic coordinate system are respectively... and .

[0081] Define the relative rotation angle between two coordinate systems: (twenty one) From the From the grid coordinate system to the first The transformation of the grid-based regional coordinate system can be decomposed into two basic rotations: First, rotate the angle around the z-axis. The corresponding four-dimensional matrix is: (twenty two) Then, rotate around the x-axis by an angle. The resulting matrix is: (twenty three) Therefore, the field rotation matrix is ​​also composed of rotations about the z-axis and rotations about the x-axis: (twenty four) This field matrix ensures the physical consistency between the transformation of the state vector and the complete three-dimensional rotational transformation of the electromagnetic field, and is a key link in the transmission calculation of connecting adjacent grids.

[0082] 3.3 Three-dimensional transfer matrix concatenation method Based on the above derivation, the first The transmission matrix of each grid in its local coordinate system is: The field rotation matrix between adjacent grid cells is For a complete propagation path consisting of N grids, the total transmission matrix of the wave from the emission point to the receiving point can be obtained by alternately multiplying each stage of the cascaded transmission matrix with the field rotation matrix: (25) This formula describes the field evolution as a wave sequentially traverses each grid in a three-dimensional curved path. It shows that when an electromagnetic wave traverses each grid, it first undergoes a propagation process within that grid, and then adapts to the local coordinate system of the next grid through coordinate rotation. This cascaded form retains the numerical accuracy of the transfer matrix method while effectively handling the directional changes in the three-dimensional path through the local coordinate system and rotation matrix mechanism, providing a unified and flexible computational framework for VLF wave propagation modeling in complex ionospheric environments.

[0083] Example 1: Prediction of long-distance VLF propagation from ground transmission to ground reception. A frequency of 15 kHz was selected, with a transmission point 3200 km away from the receiver. A three-dimensional great circle path was established using geocentric and geofixed coordinates, and a Gaussian height modulation function was introduced, with a minimum height of 0 km, a maximum height of 85 km, and a standard deviation of 0.2. Adaptive discretization along the spatial curve was performed under dual constraints of wavelength and electron density gradient, with the average grid length taken as 1 / 12 of the guided wave wavelength. A local propagation coordinate system was established within each grid, and dielectric tensor rotation transformation was performed before matrix cascade calculation. The results show that compared to the two-dimensional layered model, the field strength error rate (compared to the measured value) was reduced by 14.5%, and the phase error rate (compared to the measured value) was reduced by 17.6%, verifying the improvement effect of the three-dimensional path discretization and local coordinate rotation mechanism on curved propagation.

[0084] To verify the effectiveness of the three-dimensional path discretization algorithm of this invention in mid-range ground waveguide propagation, a 15kHz very low frequency signal was selected, and propagation calculations and experimental measurements were compared along the path from the transmitter (N 45.8°, E 126.5°) to the receiver (N 23.1°, E 113.3°). This path has a great circle distance of approximately 3200 km, covering the complete propagation process from the near-field multimode interferometry region and transition region to the far-field single-mode dominance region, providing an ideal scenario for testing the algorithm's predictive capabilities under different propagation mechanisms.

[0085] Figure 6 A comparison of the field strength calculation results from the three-dimensional algorithm of this invention, the traditional two-dimensional layered model, the ITU-R P.372-11 standard model, and measured data is presented. The typical characteristics of VLF waveguide propagation can be clearly identified from the figures. The full-field curve exhibits a three-segment structure: a rapid initial decrease, followed by fluctuating changes, and then smooth attenuation, reflecting the complex physical processes of multimode interference and single-mode dominance.

[0086] In the near-field region of 0–515 km, all curves show a rapid attenuation trend. The measured value drops rapidly from 108 dBμV / m at the emission point to 65 dBμV / m at 515 km, an attenuation of 43 dB. This stage is mainly affected by the combined influence of ground waves and low-order modes, resulting in rapid energy diffusion.

[0087] In the 515–1800 km range, the curve exhibits a distinct non-monotonic undulation characteristic, which is the most characteristic physical phenomenon in VLF waveguide propagation. The measured value is 68 dBμV / m at 515 km, then fluctuates upwards, reaching a peak of 81 dBμV / m at 1000 km, before slowly decreasing to 72 dBμV / m at 1800 km, with a fluctuation amplitude of approximately ±9–13 dB. This phenomenon originates from coherent interference caused by the coexistence of multiple waveguide modes—when the relative phase of different modes changes with distance, constructive interference (increased field strength) occurs at some distances, while destructive interference (decreased field strength) occurs at others. The three-dimensional algorithm successfully reproduced this undulation characteristic in this section, with an average deviation from the measured value of less than 1.2 dB, accurately characterizing the position and amplitude of the interference peak. While the two-dimensional algorithm and ITU-R model can reflect the overall undulation trend, there are deviations in the simulation of the peak position, with errors reaching 3–6 dB at some points.

[0088] Beyond 1800 km, higher-order modes gradually decay to nothing, and the fundamental mode begins to dominate propagation, resulting in a significant reduction in the amplitude of the propagation curve and a smoother decline. The measured value decays from 72 dBμV / m at 1800 km to 20.7 dBμV / m at 3200 km, with a total attenuation of 51 dB and an average attenuation rate of approximately 3.64 dB / 100 km. The 3D algorithm maintains a high degree of consistency with the measured values ​​throughout this region, with the field strength at the 3200 km endpoint at 24 dBμV / m, which is closest to the measured value. Compared to the 2D algorithm, the 3D algorithm significantly reduces error accumulation in far-field propagation because it takes into account the horizontal inhomogeneities of the ionosphere.

[0089] Figure 7 Comparison curves of the 3D algorithm, 2D algorithm, and measured phase data in this scenario are presented.

[0090] The phase across the entire field exhibits a combination of overall linear growth and local fluctuations. The measured phase linearly increases from 0° at 0 km to 168° at 3200 km, with an average phase gradient of 5.25° / 100 km, reflecting the fundamental phase delay of electromagnetic wave propagation. In the 100–1500 km range, the phase exhibits periodic fluctuations of approximately ±5–8°, a typical phenomenon caused by multimode interference.

[0091] The three-dimensional algorithm shows a high degree of agreement with the measured phase: the average deviation in the near-field region is less than 1°, the maximum deviation in the multimode interference region does not exceed 3°, the deviation in the far-field region is less than 1°, and the endpoint phase of 155.6° is closest to the measured value. Although the two-dimensional algorithm can reflect the linear growth trend, the maximum deviation in the multimode interference region reaches 8–10°, and the deviation gradually increases with the distance.

[0092] Statistics show that the 3D algorithm has an error rate of 7.4% compared to the measured value at the receiving point, which is significantly better than the 25% of the 2D algorithm, and the phase error rate is reduced by 17.6%. The 3D algorithm successfully reproduces the phase fluctuation characteristics caused by multimode interferometry, verifies its ability to accurately characterize complex interferometric processes, and provides a high-precision calculation tool for VLF phase prediction.

[0093] Example 2: Spaceborne Reception Scenario. The transmitter is located on the ground, and the receiver is at an altitude of 600 km with a frequency of 20 kHz. The three-dimensional path rises from the ground to the ionosphere and then extends to the satellite altitude. The altitude modulation function uses a single-peak continuous function with a peak altitude of 90 km. The local orientation correction considers the influence of the horizontal gradient of electron density and the geomagnetic cyclotron frequency to update the azimuth and inclination angles. The calculated results are compared with measured satellite data, and the amplitude prediction error is less than 2 dB, which is better than the model without orientation correction, demonstrating the technical effectiveness of the orientation correction mechanism.

[0094] To verify the effectiveness of the three-dimensional path discretization algorithm in spaceborne receiving scenarios, a comparative study was conducted on the 20 kHz VLF propagation calculation and actual measurement of ground-launched satellites at an altitude of 600 km. Figure 8 Comparison curves of the amplitudes of the three-dimensional corrected model and the uncorrected model with measured satellite data are presented.

[0095] The measured amplitude monotonically decayed from 100 dBV / m at ground level to 22.9 dBV / m at 600 km, with a total attenuation of 77.1 dB and an average attenuation rate of approximately 12.85 dB / 100 km. The attenuation was most significant in the main ionospheric region (80–120 km), reflecting enhanced absorption due to increased electron density. The three-dimensional corrected model best approximated the measured values. The uncorrected model, neglecting the horizontal gradient of electron density and the influence of geomagnetic cyclotron frequency, exhibited a maximum deviation of 9.4 dB in the main ionospheric region, with the deviation gradually increasing with altitude. These results demonstrate that the local orientation correction mechanism effectively improves the amplitude prediction accuracy in spaceborne scenarios, validating the algorithm's applicability in vertical penetration propagation.

[0096] Example 3: Propagation in a high-latitude, strong geomagnetic region. Transmission and reception are located at 65°N latitude, with a frequency of 10 kHz. Due to the large tilt angle of the geomagnetic field, the electron cyclotron frequency significantly affects the refractive index. The geomagnetic field direction is included in the correction factor calculation when constructing the local coordinate system. The rotation matrix consists of two fundamental rotations, completing the tensor transformation. Numerical results show more stable polarization mode separation and the elimination of numerical divergence, indicating that the local rotation mechanism ensures the applicability of the transfer matrix method in high-magnetic regions.

[0097] from Figure 9 As can be seen, the traditional two-dimensional algorithm without rotation mechanism drops the cross-polarization isolation to 0 dB at around 1300 km, and then becomes negative, indicating that the cross-polarization component exceeds the common polarization component, resulting in numerical divergence. In contrast, the three-dimensional algorithm using a local rotation mechanism maintains a value above 20 dB within a range of 2000 km, with stable polarization mode separation and effective suppression of the cross-polarization component.

[0098] The results show that the three-dimensional algorithm achieves an average cross-polarization isolation of 29 dB across the entire field, which is 23.2 dB higher than that of the traditional algorithm. This verifies the key role of the local rotation mechanism in polarization mode separation and numerical stability in high-latitude, strong geomagnetic regions.

[0099] Example 4: Strong Electron Density Gradient Environment. A region of abrupt change in the E layer of the solar ionosphere was selected, where the electron density gradient varied by 8% per kilometer. Path discretization adopted a gradient constraint-first principle, with the axial length taken as a threshold divided by the gradient factor. The grid length was shortened to 0.6 km. The medium within each grid was kept approximately homogeneous during the calculation. Compared with the fixed step size model, the numerical stability was significantly improved, and the oscillation error disappeared, indicating that the adaptive discretization strategy ensured physical consistency.

[0100] Figure 10 Comparison curves between the calculation results and the measured values ​​of the adaptive discrete model and the fixed step size model are presented.

[0101] As can be seen from the figure, the fixed step size model (5 km step size) exhibits significant numerical oscillations in the region of intense gradients (150–500 km), with a maximum deviation of 4.1 dB, and the oscillations persist with distance. The adaptive discrete model automatically densifies the grid based on gradient constraints (minimum step size 0.6 km), and the calculation results are smooth and stable, with a deviation from the reference value consistently less than 0.6 dB.

[0102] The results show that the adaptive discretization strategy effectively guarantees the assumption of homogeneity of the medium within the grid by prioritizing gradient constraints, eliminates numerical oscillations, and significantly improves computational stability and physical consistency under strong gradient conditions.

[0103] Example 5: Long-distance transhemispheric propagation. Propagation distance: 7500 km, frequency: 8 kHz. Peak path height: 95 km. Total number of discrete nodes: 420. Local coordinate system continuously updated along the path, rotation matrix recursively derived segment by segment. Cascaded matrix calculation outputs field strength distribution. Comparison with the international ionospheric reference model shows consistent trends with a deviation of less than 10%, indicating that the three-dimensional curved path and tensor rotation mechanism are suitable for ultra-long-distance propagation.

[0104] To verify the applicability of the three-dimensional algorithm in ultra-long-distance propagation, propagation calculations were carried out along an east-west transoceanic path (7500 km, 8 kHz) from Shanghai (31.2°N, 121.5°E) to Los Angeles (34.1°N, 118.2°W). Figure 11 Comparison curves of field strength variation with distance between the 3D algorithm and the IRI model are presented.

[0105] The graph shows that the attenuation curve roughly follows a three-segment structure: the attenuation is relatively rapid in the nearshore region (0–1200 km) (approximately 13 dB / 1000 km), the attenuation is gradual in the open ocean region (1200–6000 km) (approximately 9–10 dB / 1000 km), and the attenuation rate further decreases to approximately 5 dB / 1000 km in the nearshore region (6000–7500 km). This characteristic is consistent with the theoretical laws governing the east-west propagation of VLF.

[0106] The three-dimensional algorithm and the IRI model curves showed a relatively consistent trend with an average relative deviation of 7.6%. The field strengths at the 7500 km endpoint were -70.5 dB and -64.2 dB, respectively, with a relative deviation of 9.8%, which verified the effectiveness and accuracy of the algorithm in ultra-long-distance cross-hemispheric propagation.

[0107] Example 6: Multi-band Comparison and Verification. The frequencies were 5 kHz, 12 kHz, and 25 kHz, with all other conditions remaining the same. Different frequency field distributions were calculated using the same three-dimensional path discretization framework. Wavelength variations automatically affected the selection of the axial grid size, with larger grid lengths for low frequencies and automatically shorter grid lengths for high frequencies. The computation time maintained a linear increase, and no numerical instability was observed. The results demonstrate that the modeling structure has adaptive capabilities to frequency variations, verifying its engineering feasibility and universality.

[0108] To verify the adaptive capability of the three-dimensional algorithm to frequency changes, path propagation calculations over 2000 km were performed at three frequencies: 5 kHz, 12 kHz, and 25 kHz. Figure 12 Comparative curves of field strength versus distance at different frequencies are presented. The images show that the higher the frequency, the greater the attenuation rate and the steeper the curve, consistent with the attenuation rate ∝ f in the VLF band. 2 The physical laws governing this process are as follows. Different frequency wavelengths influence the selection of the axial grid size. For a 5 kHz wavelength of 60 km, the grid size is 30 km (67 grids, N=2); for a 12 kHz wavelength of 25 km, the grid size is 12.5 km (160 grids, N=2); and for a 25 kHz wavelength of 12 km, the grid size is 12 km (167 grids, N=1). The grid size is automatically adjusted according to the wavelength, ensuring the calculation accuracy for each frequency band.

[0109] Figure 13 The relative computation time for the field distribution at different frequencies is presented. The results show that the computation time increases linearly with the number of grid cells (1.0 at 5 kHz, 2.4 times at 12 kHz, and 2.5 times at 25 kHz), without exponential explosion or numerical divergence. The results demonstrate that the three-dimensional algorithm has good adaptability to frequency changes and engineering versatility.

[0110] It should be noted that embodiments of the present invention can be implemented in hardware, software, or a combination of both. The hardware portion can be implemented using dedicated logic; the software portion can be stored in memory and executed by a suitable instruction execution system, such as a microprocessor or dedicated-design hardware. Those skilled in the art will understand that the above-described devices and methods can be implemented using computer-executable instructions and / or included in processor control code, for example, such code provided on a carrier medium such as a disk, CD, or DVD-ROM, a programmable memory such as read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The devices and modules of the present invention can be implemented by hardware circuitry such as very large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, or programmable hardware devices such as field-programmable gate arrays, programmable logic devices, etc., or by software executed by various types of processors, or by a combination of the above-described hardware circuitry and software, such as firmware.

[0111] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions, and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention, and within the spirit and principles of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for calculating the transfer matrix based on a three-dimensional path discretization ionospheric model, characterized in that, Includes the following steps: Step 1: Obtain the spatial coordinates of the transmitting and receiving points, establish a spatial rectangular coordinate system with the Earth's center as the origin, and construct a three-dimensional propagation path connecting the transmitting and receiving points; Step 2: Establish a height distribution function on the three-dimensional propagation path, assign a propagation height that changes continuously with the path parameters to each location point on the path, and form a three-dimensional propagation trajectory in the form of a spatial curve; Step 3: Discretize the three-dimensional propagation trajectory according to the axial scale adaptive principle to generate a spatial node sequence arranged in ascending order of arc length; Step four: Using the line connecting adjacent spatial nodes as the axis, establish local grid cells in each discrete segment, so that the medium parameters in each grid cell satisfy the condition of approximate uniformity, and use them as the basic unit for subsequent transmission matrix calculation.

2. The method as described in claim 1, characterized in that, The spatial rectangular coordinate system is a geocentric coordinate system with the Earth's center of mass as the origin. The transmitting and receiving points are obtained by converting geographical longitude, geographical latitude, and vertical altitude. The conversion relationship includes: Calculate the radius of curvature of the ellipsoid, which is the Earth's semi-major axis divided by the square root of 1 minus the square of the first eccentricity multiplied by the square of the latitude sine. The X component of spatial coordinates is the sum of the radius of curvature and the height, multiplied by the cosine of latitude, and then multiplied by the cosine of longitude. The Y component of the spatial coordinates is the sum of the radius of curvature and the height, multiplied by the cosine of the latitude, and then multiplied by the sine of the longitude. The Z component of the spatial coordinates is the sum of the radius of curvature multiplied by 1, the square of the first eccentricity, and the altitude, multiplied by the latitude sine.

3. The method as described in claim 1, characterized in that, The planar reference path of the three-dimensional propagation path is represented by a geocentric great circle path, and the position vector of any point on the path is represented by spherical linear interpolation, which is expressed as follows: The path point vector equals the first endpoint vector multiplied by the sine function 1 minus the path parameter multiplied by the geocentric angle divided by the geocentric angle sine, plus the second endpoint vector multiplied by the sine function path parameter multiplied by the geocentric angle divided by the geocentric angle sine, where the path parameter takes values ​​from 0 to 1.

4. The method as described in claim 1, characterized in that, The height distribution function adopts a unimodal continuous function form, and its value is equal to the minimum height plus the maximum height minus the minimum height multiplied by the shape function. The shape function takes a value of 0 when the path parameter is 0 and 1, and a value of 1 when the path parameter is 0.

5.

5. The method as described in claim 4, characterized in that, The shape function is a Gaussian function, which is expressed as the negative path parameter of the exponential function minus 0.5 squared, divided by 2, and multiplied by the square of the standard deviation, where the standard deviation is 0.

2.

6. The method as described in claim 1, characterized in that, The axial dimension adaptive principle includes wavelength constraints and medium gradient constraints, with the axial length taken as the smaller of the two calculated values, wherein: The wavelength constraint length is equal to the guided wave wavelength divided by the number of sampling points per wavelength; The length of the medium gradient constraint is equal to the preset relative change threshold divided by the electron density gradient factor.

7. A method for calculating the local coordinate rotation transfer matrix of a three-dimensional path-discrete ionospheric model, implementing the transfer matrix calculation method based on any one of claims 1-6, characterized in that, Includes the following steps: Step 1: For each discretized grid cell, construct a local rectangular coordinate system with the grid axis as the propagation direction. The local coordinate system includes a first axis along the propagation direction, a second axis located in the plane formed by the propagation direction and the geographic perpendicular direction and perpendicular to the first axis, and a third axis determined according to the right-hand rule. Step 2: Calculate the azimuth and tilt corrections for the propagation direction based on the electron density at the grid midpoint and the geomagnetic field parameters, and update the local propagation direction. Step 3: Construct a rotation matrix from the geographic coordinate system to the local coordinate system. The rotation matrix consists of two basic rotations: the first rotation is around the geographic vertical axis to rotate the azimuth angle, and the second rotation is around the updated horizontal axis to rotate the tilt angle. Step 4: Transform the relative permittivity tensor in the geographic coordinate system to the local coordinate system using the rotation matrix to obtain the local equivalent relative permittivity tensor. Step 5: Establish a hierarchical transmission matrix in the local coordinate system and perform matrix concatenation calculations in grid order.

8. The method as described in claim 7, characterized in that, The azimuth correction is related to the electron density horizontal gradient, refractive index, electron cyclotron frequency, and grid axial length. The correction is equal to the proportionality constant multiplied by the electron density gradient, divided by the refractive index, multiplied by the geomagnetic correction factor, and then multiplied by the axial length.

9. The method as described in claim 7, characterized in that, The relative permittivity tensor in the local coordinate system is represented by a rotation matrix multiplied by the geographic coordinate system relative permittivity tensor, and then multiplied by the transpose of the rotation matrix.

10. A transfer matrix calculation system based on a three-dimensional path discretized ionospheric model, implementing the transfer matrix calculation method based on any one of claims 1-6, characterized in that, include: The path modeling module is used to construct a three-dimensional propagation path based on the coordinates of the transmitter and receiver points, and generate a spatial propagation trajectory that includes height distribution. The path discretization module is used to adaptively discretize the spatial propagation trajectory according to wavelength constraints and medium gradient constraints to generate a grid cell sequence. The local coordinate construction module is used to establish a local Cartesian coordinate system consistent with the propagation direction within each grid cell; The medium transformation module is used to construct the rotation matrix from the geographic coordinate system to the local coordinate system and to complete the coordinate transformation of the relative permittivity tensor. The transmission matrix calculation module is used to establish the unit transmission matrix in the local coordinate system and perform sequential cascade operations to output the overall transmission characteristics from the transmitting point to the receiving point.

Images