Short-wave single-station direction-finding and elevation-finding combined positioning method and system

Abstract

The application discloses a short-wave single-station direction-finding and elevation-measuring combined positioning method and system, extracts a weighted azimuth angle sequence and an elevation angle confidence interval from array receiving data synchronously, compresses a reliable range of the azimuth angle to construct a coming wave angle parameter set by using a coupling relationship between the azimuth and the elevation; reversely calculates ionospheric height real-time deviation from a multipath time delay difference component of the receiving data and performs self-correction on the ionospheric virtual height, to correct a virtual height value and an elevation angle estimation value combined with the coming wave angle parameter set to jointly drive a reflection slant distance solution to generate a candidate positioning coordinate set; identifies a multi-hop propagation component, extracts each hop ground reflection constraint, implements geometric consistency verification on the candidate coordinates to form a reliability grading and optimizes a slant distance configuration according to the reliability grading; finally, the candidate coordinates are weighted and clustered in cooperation with the reliability grading and the slant distance optimization configuration to construct a multi-dimensional positioning parameter set, and the final positioning coordinates of a radiation source are determined through fusion solution, so that more reliable position solution results can be provided for short-wave single-station direction-finding and elevation-measuring combined positioning.

Description

Technical Field

[0001] This invention relates to the field of radio signal processing technology, and in particular to a shortwave single-station direction finding and elevation measurement combined positioning method and system. Background Technology

[0002] Shortwave electromagnetic waves can propagate beyond line of sight by means of ionospheric refraction, and have extensive engineering applications in long-range radio monitoring and radiation source location estimation. Single-station passive positioning systems rely solely on observation data from a single receiving station to estimate the azimuth and distance of radiation sources. This offers advantages such as flexible deployment and no need for time synchronization between multiple stations, making it suitable for applications with limited site selection.

[0003] In a single-station system, obtaining distance information depends on the accuracy of ionospheric reflection height. However, the actual ionospheric height fluctuates continuously due to factors such as day-night cycles and solar activity, easily leading to discrepancies between the reflection path length calculated using empirical models and the actual value. In the array direction finding process, there is a coupling effect between the estimation of horizontal azimuth and elevation angles. Under low elevation conditions, azimuth estimation is particularly sensitive to elevation deviations, further increasing the uncertainty in position calculation. Long-distance propagation of shortwave signals is often accompanied by multi-hop path effects, with different hop count components carrying different path geometric information. Existing schemes mostly drive calculations with fixed ionospheric parameters, lacking adaptive tracking of real-time ionospheric fluctuations and failing to fully utilize the geometric constraints inherent in multi-hop paths to cross-validate positioning results, thus limiting the applicability of single-station positioning in complex ionospheric environments. Summary of the Invention

[0004] This invention provides a shortwave single-station direction-finding and elevation-measuring joint positioning method and system. By synchronously analyzing the incoming wave angle information of the array received data and narrowing the angle search range using the azimuth-elevation coupling relationship, the real-time deviation of the ionospheric height is calculated by back-calculating the multipath delay difference component to correct the reflection path parameters online. Combined with the geometric constraints of the multi-hop propagation components, the candidate coordinates are subjected to credibility classification and weighted clustering fusion, thereby improving the accuracy and robustness of radiation source location calculation under single-station conditions.

[0005] The first aspect of this invention proposes a shortwave single-station direction finding and elevation finding combined positioning method, comprising the following steps: Acquire shortwave antenna array received data, and simultaneously extract weighted azimuth sequence and elevation confidence interval based on the shortwave antenna array received data. Compress the azimuth confidence range based on the weighted azimuth sequence and the elevation confidence interval to construct the arrival angle parameter set. Based on the arrival angle parameter set, the elevation angle estimate is extracted. The multipath delay difference component is separated from the data received by the shortwave antenna array to reversely calculate the real-time ionospheric height deviation. Based on the real-time ionospheric height deviation, the ionospheric virtual height is self-corrected to obtain the corrected virtual height value. The corrected virtual height value and the elevation angle estimate are combined with the arrival angle parameter set to jointly calculate the reflection slant range and determine the candidate positioning coordinate set. In the shortwave antenna array received data, multi-hop propagation components are identified and ground reflection constraints of each hop are extracted. Geometric consistency verification is performed between each hop ground reflection constraint and the candidate positioning coordinate set to form a coordinate confidence level table. The candidate positioning coordinate set is then weighted and optimized according to the coordinate confidence level table to establish a slant range optimization configuration. Based on the coordinate confidence level table and the slant distance optimization configuration, the candidate positioning coordinate set is collaboratively weighted and clustered to construct a multidimensional positioning parameter set. Based on the multidimensional positioning parameter set, the radiation source location fusion calculation is performed to determine the final positioning coordinates.

[0006] A second aspect of this invention provides a shortwave single-station direction finding and elevation measurement combined positioning system, comprising: An angle extraction module is used to acquire shortwave antenna array received data, and simultaneously extract weighted azimuth sequence and elevation confidence interval based on the shortwave antenna array received data. Based on the weighted azimuth sequence and the elevation confidence interval, the azimuth confidence range is compressed to construct the incoming wave angle parameter set. The virtual height correction module is used to extract the elevation angle estimate based on the arrival angle parameter set, separate the multipath delay difference component from the data received by the shortwave antenna array to back-calculate the real-time ionospheric height deviation, self-correct the ionospheric virtual height based on the real-time ionospheric height deviation to obtain the corrected virtual height value, and combine the corrected virtual height value with the elevation angle estimate and the arrival angle parameter set to jointly perform reflection slant range calculation to determine the candidate positioning coordinate set. The multi-hop verification module is used to identify multi-hop propagation components in the shortwave antenna array received data, extract ground reflection constraints for each hop, perform geometric consistency verification between the ground reflection constraints for each hop and the candidate positioning coordinate set to form a coordinate confidence level table, and establish a slant range optimization configuration by weighting the candidate positioning coordinate set according to the coordinate confidence level table. The fusion positioning module is used to construct a multi-dimensional positioning parameter set by collaborative weighted clustering of the candidate positioning coordinate set based on the coordinate confidence level table and the slant range preferred configuration, and to perform radiation source location fusion calculation based on the multi-dimensional positioning parameter set to determine the final positioning coordinates.

[0007] The beneficial effects of this invention are reflected in the following points: 1. For each elevation sampling point within the elevation confidence interval, the absolute value of the partial derivative of the azimuth estimation error with respect to the elevation angle is calculated one by one. A quantitative mapping between the elevation error amplification factor sequence and the upper bound of the azimuth coupling error is established. A compressed azimuth distribution is generated and jointly encapsulated with the elevation confidence interval to construct the arrival angle parameter set, realizing the transformation of the azimuth-elevation coupling effect from the source of positioning error to the basis for active compression of the angle confidence range; 2. The time-varying slope is extracted window by window from the multipath delay difference component of the received data. The real-time deviation of the ionospheric virtual height is calculated in reverse and corrected by window by window. The dependence of the slant range calculation on external ionospheric detection data is eliminated, and the self-consistent correction of the angle input and path parameters within the same data frame is realized; 3. A hop-by-hop delay difference ratio is constructed. The sequence is extracted and the linear trend of the proportional deviation is established. The trend slope is geometrically mapped to the magnitude and direction of the ionospheric horizontal gradient. Differential gradient compensation is applied to each hop reflection point hop by hop and then extended in the reverse direction to form the spatial constraint of the radiation source. The spatial overlap is defined as the credibility quantification index of the candidate positioning coordinate set. The transformation of the multi-hop propagation redundancy path from positioning interference to geometric cross-validation active resource is realized. 4. The clustering residual is defined as the evaluation deviation quantification basis of the initial credibility classification. The reverse mapping between the residual distribution and the credibility level adjustment direction is established. The second round of weighted clustering is driven according to the corrected credibility structure. The credibility classification is upgraded from a one-time static weight to an iterative self-correcting dynamic constraint mechanism, so that the converged positioning coordinates are established on the credibility structure verified by the residual. Attached Figure Description

[0008] The accompanying drawings illustrate specific examples of the technical solutions described in this invention and, together with the detailed embodiments, form part of the specification, serving to explain the technical solutions, principles, and effects of this invention.

[0009] Unless otherwise specified or defined, the same reference numerals in different figures represent the same or similar technical features, and different reference numerals may be used to represent the same or similar technical features.

[0010] Figure 1 This is a flowchart illustrating a shortwave single-station direction finding and elevation measurement combined positioning method according to the present invention.

[0011] Figure 2 This is a structural block diagram of a shortwave single-station direction finding and elevation measurement combined positioning system according to the present invention. Detailed Implementation

[0012] In the following description, specific details such as particular system architectures and techniques are set forth for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application may also be implemented in other embodiments without these specific details. In other instances, detailed descriptions of well-known systems, apparatuses, circuits, and methods have been omitted so as not to obscure the description of this application with unnecessary detail.

[0013] It should be understood that, when used in this application specification and the appended claims, the term "comprising" indicates the presence of the described features, integrals, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, integrals, steps, operations, elements, components and / or a collection thereof.

[0014] References to "one embodiment" or "some embodiments" as described in this specification mean that one or more embodiments of this application include a specific feature, structure, or characteristic described in connection with that embodiment. Therefore, the phrases "in one embodiment," "in some embodiments," "in other embodiments," "in still other embodiments," etc., appearing in different parts of this specification do not necessarily refer to the same embodiment, but rather mean "one or more, but not all, embodiments," unless otherwise specifically emphasized. The terms "comprising," "including," "having," and variations thereof mean "including but not limited to," unless otherwise specifically emphasized.

[0015] The technical solutions of the embodiments of this application will be described below.

[0016] like Figure 1 As shown, this embodiment of the invention provides a shortwave single-station direction finding and elevation measurement combined positioning method, including the following steps S110-S140: Step S110: Obtain the received data from the shortwave antenna array, and simultaneously extract the weighted azimuth sequence and elevation confidence interval based on the received data from the shortwave antenna array. Compress the azimuth confidence range based on the weighted azimuth sequence and elevation confidence interval to construct the arrival angle parameter set.

[0017] Specifically, the data received by the shortwave antenna array is acquired. The shortwave antenna array consists of multiple array elements arranged according to a preset geometric configuration. The spacing between the array elements is set according to the half-wavelength condition of the operating frequency band to meet spatial sampling requirements. The array configuration provides an azimuth aperture in the horizontal plane while having a height difference component in the vertical direction to support elevation angle extraction. The electromagnetic signals received by each array element are amplified by a low-noise amplifier and then digitally sampled by a synchronous analog-to-digital converter of each channel using a unified sampling clock. The sampling clock of each channel is driven by the same reference oscillator to ensure that the alignment error of the sampling time between array elements is much lower than the sampling period. The digital sampling sequence of each array element is stored in frames according to a preset frame length. The frame length is an integer fraction of the coherence time of the shortwave signal to ensure that the statistical characteristics of the signal within the frame are approximately stable. Frames are continuously acquired at fixed frame intervals to form a time-continuous shortwave antenna array received data stream. Each frame of data contains the synchronous sampling matrix of all array elements and is accompanied by a frame number and an absolute timestamp. Before data acquisition, each array element receiving channel undergoes channel consistency calibration. The calibration coefficient matrix performs array element-by-array amplitude and phase compensation on the sampling matrix of each frame of the shortwave antenna array received data to eliminate the influence of channel hardware differences on the extraction of incoming wave spatial features. The shortwave antenna array received data is in the form of the calibrated and compensated sampling matrix sequence of each frame, covering all monitoring frequency points within the working frequency band, and is organized by dual indexing of frequency point index and frame sequence number.

[0018] Weighted azimuth sequence and elevation confidence interval are extracted synchronously based on data received by shortwave antenna array. The extraction of the weighted azimuth sequence involves constructing spatial covariance matrices and performing spatial spectrum estimation on the shortwave antenna array received data at multiple monitoring frequencies. The spatial covariance matrices at each frequency are separated into signal and noise subspaces through eigenvalue decomposition. The peak position of the spatial spectrum function at each frequency corresponds to the estimated azimuth angle value at that frequency. The ratio of peak power to noise power is the signal-to-noise ratio (SNR) of the estimated azimuth angle at that frequency. The estimated azimuth angle values ​​at each frequency are normalized with SNR and then weighted and combined to form a single-frame weighted azimuth angle value. The weighted combination formula is φ_w=Σ_f(SNR_f×φ_f) / Σ_f(SNR_f), where φ_w is the single-frame weighted azimuth angle value (unit: °), SNR_f is the SNR of the f-th monitoring frequency (dimensionless), and φ_f is the estimated azimuth angle value of the f-th monitoring frequency (unit: °). The weighted azimuth angle values ​​of each frame are arranged according to the frame number to form a weighted azimuth angle sequence. The elevation confidence interval is extracted using a subset of array elements with vertical aperture components. Vertical aperture components refer to pairs or groups of elements in the array configuration that have a height difference. A vertical-dimensional spatial covariance matrix is ​​constructed for each frame of data from the vertical aperture element subset received by the shortwave antenna array, and spatial spectrum estimation is performed within the elevation search range. The position of the main peak of the elevation-dimensional spatial spectrum function corresponds to the center value of the elevation estimation, and the half-power width of the main peak determines the upper and lower boundaries of the elevation confidence interval. After calculating the elevation-dimensional spatial spectrum function at multiple frequency points, the intersection of the main peak intervals at each frequency point is taken as the final elevation confidence interval. This intersection operation narrows the elevation confidence interval compared to the single-frequency-point result. The elevation confidence interval and the weighted azimuth sequence are extracted synchronously on each frame of data and share the frame number index and timestamp. Synchronous extraction ensures that the azimuth and elevation angle information are strictly aligned in the time dimension.

[0019] In some embodiments, constructing the arrival angle parameter set by compressing the azimuth confidence range based on the weighted azimuth sequence and the elevation confidence interval includes: performing low elevation disturbance sensitivity analysis on the elevation confidence interval to generate an elevation error amplification factor sequence; performing error coupling quantization on the weighted azimuth sequence based on the elevation error amplification factor sequence to obtain an upper bound of azimuth coupling error; performing interval contraction on the weighted azimuth sequence according to the upper bound of azimuth coupling error to generate a compressed azimuth distribution; and jointly encapsulating the compressed azimuth distribution with the elevation confidence interval to construct the arrival angle parameter set.

[0020] A low-elevation-angle perturbation sensitivity analysis is performed on the elevation confidence interval to generate an elevation error amplification factor sequence. The amplification of the azimuth estimation error is uneven across different elevation angle values ​​within the elevation confidence interval; small perturbations at low elevation angles have a much stronger impact on azimuth estimation than those at high elevation angles. The low-elevation-angle perturbation sensitivity analysis calculates the absolute value of the partial derivative of the azimuth estimation error with respect to elevation angle for each elevation angle sampling point within the elevation confidence interval. This absolute value of the partial derivative is the elevation error amplification factor corresponding to that sampling point, defined as... Where k(θ) is the elevation error amplification factor (unit: ° / °) at elevation angle θ (unit: °), and φ is the azimuth estimate (unit: °). The partial derivatives are calculated numerically at each sampling point in the elevation confidence interval under the constraints of the geometric propagation model. The elevation error amplification factor sequence records the k(θ) results at each point with the elevation sampling point as the index. k(θ) near the lower boundary of the elevation confidence interval is usually significantly larger than that in the middle of the interval under the geometric constraints of radio wave grazing propagation. When shortwave signals are refracted by low elevation angles in the ionosphere, the propagation geometry is extremely sensitive to elevation angle disturbances. Small measurement deviations in elevation angle are significantly mapped to azimuth estimation errors. However, when the signal crosses the ionosphere at a higher elevation angle, almost vertically, the impact of the same elevation angle deviation on azimuth estimation is relatively limited. When the elevation confidence interval is wide, the elevation error amplification factor sequence covers a wider range of elevation angles, and the distribution span of factor values ​​at each sampling point is also larger. When the elevation confidence interval narrows to the high elevation angle segment, the peak factor of the elevation error amplification factor sequence decreases accordingly, and the overall distribution of the sequence tends to be flat. The lower boundary elevation angle provided by the elevation confidence interval determines the peak factor in the elevation error amplification factor sequence. The peak value of the elevation error amplification factor sequence corresponding to the elevation confidence interval containing low elevation angle values ​​can be several times or even more than ten times the factor value at the high elevation angle end. The elevation error amplification factor sequence completely describes the amplification law of the azimuth estimation error under each elevation angle assumption within the elevation confidence interval.

[0021] The upper bound of the azimuth coupling error is obtained by performing error coupling quantization on the weighted azimuth sequence based on the elevation error amplification factor sequence. The peak factor of the elevation error amplification factor sequence corresponds to the elevation angle condition that maximizes the azimuth estimation error within the elevation confidence interval. Error coupling quantization uses the elevation error amplification factor sequence and the elevation confidence interval as joint inputs. The maximum error offset of the weighted azimuth sequence under the most unfavorable elevation angle assumption is calculated by combining the peak factor of the elevation error amplification factor sequence with the width of the elevation confidence interval. The upper bound of the azimuth coupling error is defined as follows: Where Δφ_max is the upper bound of the azimuth coupling error (unit: °), k(θ) is the factor value of the elevation error amplification factor sequence at the elevation angle θ (unit: ° / °), θ_L and θ_U are the lower and upper boundaries of the elevation confidence interval, respectively (unit: °), and max{k(θ)} takes the maximum value of the factor value within the elevation confidence interval [θ_L, θ_U]. The error offset of each azimuth estimate in the weighted azimuth sequence increases monotonically with the factor value of the elevation error amplification factor sequence. When the factor value in the elevation error amplification factor sequence is high, the offset of the corresponding azimuth estimate increases significantly. The upper envelope of the offset of each azimuth estimate constitutes the upper bound of the azimuth coupling error. The coverage of the upper bound of the azimuth coupling error is determined by the concentration of entries in the weighted azimuth sequence and the peak value of the elevation error amplification factor, both of which are encapsulated in the calculation result of the upper bound of the azimuth coupling error. When there are multiple peak factors in the elevation error amplification factor sequence, the upper bound of the azimuth coupling error is taken as the maximum value of the envelope of the offsets corresponding to all peaks, ensuring that the upper bound of the azimuth coupling error effectively covers the error spread range of the weighted azimuth sequence under the assumption of all elevation angles in the elevation confidence interval.

[0022] The weighted azimuth sequence is compressed by applying interval shrinkage to generate a compressed azimuth distribution based on the upper bound of the azimuth coupling error. The upper bound of the azimuth coupling error defines the maximum offset boundary of each azimuth estimate in the weighted azimuth sequence due to elevation uncertainty. Interval shrinkage combines each azimuth estimate in the weighted azimuth sequence with the upper bound of the azimuth coupling error, retaining only entries where the azimuth estimate still falls within the central interval of the overall weighted azimuth sequence distribution under the offset boundary constraint. The main estimates of the weighted azimuth sequence are concentrated towards the true direction of arrival. Deviation estimates caused by multipath interference or ionospheric disturbances are identified and eliminated because they deviate from the central interval beyond the upper bound of the azimuth coupling error. Therefore, the compressed azimuth distribution is more reliably concentrated near the true direction of arrival. The shrinkage retention conditions for azimuth estimates with higher weights in the weighted azimuth sequence are relatively lenient, while marginal estimates with lower weights are more easily eliminated after shrinkage is applied to the upper bound of the azimuth coupling error. Thus, the angular coverage range of the compressed azimuth distribution is significantly narrower than the original range of the weighted azimuth sequence. When multiple azimuth estimates cluster in the same angular interval, the upper bound of the azimuth coupling error has a relatively small shrinkage effect on this clustered interval, and the compressed azimuth distribution retains a relatively complete set of azimuth entries within this interval. Low-weight estimates of the weighted azimuth sequence at isolated angles are no longer retained in the compressed azimuth distribution after interval shrinkage. The entry density distribution of the compressed azimuth distribution reflects the reliable concentrated region of the weighted azimuth sequence under the constraint of the upper bound of the azimuth coupling error. The interval shrinkage operation transforms the error amplification effect of elevation uncertainty on azimuth into a quantifiable basis for elimination. All entries retained in the compressed azimuth distribution have passed the boundary test of the upper bound of the coupling error.

[0023] The wave arrival angle parameter set is constructed by jointly encapsulating the compressed azimuth distribution and the elevation confidence interval. The compressed azimuth distribution provides the reliable angle range in the azimuth dimension and the weight distribution of each angle entry, while the elevation confidence interval provides the angle constraint boundary in the elevation dimension. The joint encapsulation of these two elements results in the wave arrival angle parameter set possessing angle constraint capabilities in both the azimuth and elevation dimensions. The wave arrival angle parameter set uses the angle entries from the compressed azimuth distribution as its core, associating each azimuth entry with the upper and lower boundaries of the elevation confidence interval, forming a two-dimensional angle coverage structure with joint azimuth and elevation constraints. The upper and lower boundaries of the elevation confidence interval uniformly apply elevation constraints to all azimuth entries in the wave arrival angle parameter set. The weight structure of the compressed azimuth distribution remains unchanged in the wave arrival angle parameter set. The angle coverage density of the wave arrival angle parameter set is determined by the concentration of entries in the compressed azimuth distribution; the sparser the entries in the compressed azimuth distribution, the smaller the effective angle coverage interval of the wave arrival angle parameter set. The incoming wave angle parameter set solidifies the joint expression of the compressed azimuth distribution and elevation confidence interval into a single data structure. The subsequent reflection slant range calculation and geometric consistency verification are both based on the incoming wave angle parameter set as a unified angle input source. The angle boundaries of the azimuth and elevation dimensions have been processed by error coupling compression.

[0024] Step S120: Extract elevation angle estimate based on arrival angle parameter set, separate multipath delay difference component from shortwave antenna array received data to reverse calculate real-time ionospheric height deviation, obtain corrected virtual height value by self-correcting ionospheric virtual height based on real-time ionospheric height deviation, and combine corrected virtual height value with elevation angle estimate to perform reflection slant range calculation to determine candidate positioning coordinate set.

[0025] Specifically, elevation angle estimates are extracted based on the arrival angle parameter set. Elevation angle sampling points are set at uniform intervals within the elevation angle confidence interval of the arrival angle parameter set. The elevation angle-dimensional spectral function value corresponding to each sampling point is used as the spectral energy weight. The spectral energy weight is multiplied together with the reciprocal of the corresponding factor value of the elevation angle error amplification factor sequence in the arrival angle parameter set to form a comprehensive weight, defined as w_j = E_j / k(θ_j), where w_j is the comprehensive weight of the j-th elevation angle sampling point (dimensionless), E_j is the elevation angle-dimensional spectral function value corresponding to that sampling point (dimensionless), and k(θ_j) is the elevation angle error amplification factor (unit: ° / °) corresponding to that sampling point in the arrival angle parameter set. When the elevation angle confidence interval simultaneously covers both low and medium-high elevation angle sampling points, the spectral energy of the low elevation angle sampling points may be stronger, but their corresponding error amplification factor is extremely large. The comprehensive weight mechanism automatically shifts the elevation angle estimate towards the medium-high elevation angle sampling points, where the error amplification effect is weaker, reducing the adverse influence of low elevation angle refraction geometry on the final positioning result. The elevation angle estimate is calculated by taking the weighted average of all sampling points using the comprehensive weighting formula: θ_e = Σ_j(w_j × θ_j) / Σ_j(w_j), where θ_e is the elevation angle estimate (unit: °), and θ_j is the elevation angle value of the j-th sampling point (unit: °). The more concentrated the comprehensive weight distribution, the closer the elevation angle estimate is to the sampling point corresponding to the weight peak; a flat distribution approaches the geometric center of the elevation angle confidence interval. The elevation angle estimate is accompanied by a weighted standard deviation as a measure of uncertainty.

[0026] In some embodiments, the step of separating multipath delay difference components from the data received by the shortwave antenna array to inversely calculate the real-time ionospheric height deviation includes: performing time-delay domain separation on the data received by the shortwave antenna array to extract the primary and secondary path components; generating a path delay difference sequence by differentiating the primary and secondary path components; extracting a time-varying slope from the path delay difference sequence to generate an ionospheric drift rate estimate; and performing dynamic weighted calculation on the path delay difference sequence based on the ionospheric drift rate estimate to generate the real-time ionospheric height deviation.

[0027] The received data from the shortwave antenna array is separated in the time-delay domain to extract the principal and secondary path components. The sampling sequences of each element of the received data are mapped to the time-frequency domain via windowed Fourier transform. The windowing process uses either a Hanning window or a Blackman window to achieve a balance between time-delay resolution and sidelobe suppression. The energy distribution in the time-frequency domain along the time-delay axis forms the time-delay power spectrum of each element. The principal path component corresponds to the energy peak with the shortest time delay in the time-delay power spectrum, while the secondary path component corresponds to the independent peak with the second shortest time delay among the energy peaks in the time-delay power spectrum, whose peak position distance from the principal path component exceeds the minimum resolvable time delay. The minimum resolvable time delay is determined by the signal bandwidth of the received data from the shortwave antenna array; the wider the signal bandwidth, the smaller the minimum resolvable time delay, and the more relaxed the resolvable conditions for the principal and secondary path components. The received data from the shortwave antenna array exhibits a slight time delay deviation along the element arrangement direction in the peak position of the main path component of the time delay power spectrum of each element. This deviation reflects the phase difference between elements caused by the direction of arrival. The main path component is extracted based on the median value of the peak position of the entire element time delay power spectrum. The main path peak position of the element time delay power spectrum with a deviation exceeding one time delay resolution unit is replaced by the median value. This replacement process ensures that the consistency of the main path component peak position across all elements is not affected by the abnormal deviation of individual elements. After removing the peak position of the primary path component from the time-delay power spectrum, the secondary path component is determined by selecting the first peak position that meets the resolvability condition in the remaining energy peaks in ascending time-delay order. The confirmation of the secondary path component peak position also requires that the peak power exceeds the noise floor to eliminate interference from false peak positions. When the secondary path component is missing in the shortwave antenna array received data, the array element will not generate a secondary path entry in the time-delay domain separation result. Array elements with missing secondary path entries will not participate in the subsequent cross-element time delay difference weighted average calculation. The extraction of the primary and secondary path components is performed independently on each frame of the shortwave antenna array received data. Time-delay domain separation automatically terminates at the component extraction result corresponding to the last energy peak that meets the resolvability condition when the resolvability condition is not met.

[0028] The path delay difference sequence is generated by differentiating the primary and secondary path components. The delay peak positions of the primary and secondary path components are calculated element-by-element for each array element. Before the difference operation, the peak positions of the primary and secondary path components are refined to sub-sampling precision using quadratic interpolation to reduce quantization errors introduced by discrete sampling. The path delay difference of all array elements is weighted and averaged using the signal-to-noise ratio (SNR) to generate a single-frame path delay difference value. Array elements with higher SNR contribute more weight to the average to improve the estimation stability of the single-frame path delay difference value. The path delay difference sequence continuously records the single-frame path delay difference values ​​using the observation frame number as an index. The values ​​of the path delay difference sequence fluctuate slightly around a certain central delay under normal ionospheric conditions. The propagation path height difference between the primary and secondary path components determines the baseline level of the path delay difference sequence. When the single-frame path delay difference value of the path delay difference sequence is larger than the baseline level, it indicates that the secondary path component has undergone refraction and return from a higher layer. Frames with extremely low signal-to-noise ratios (SNR) exhibit significant peak extraction errors in their principal path component, resulting in low reliability of the path delay difference sequence entries. The path delay difference sequence is categorized into three reliability levels: high, medium, and low, using an SNR threshold to assign a reliability level to each frame entry. High-reliability entries have an SNR above the upper quartile of the entire sequence, low-reliability entries have an SNR below a preset threshold, and medium-reliability entries have an SNR between the preset threshold and the upper quartile of the entire sequence. High-reliability entries participate with full weight in subsequent time-varying slope fitting, medium-reliability entries participate with a 0.5 weight, and low-reliability entries participate with reduced weight in subsequent time-varying slope extraction. The inter-frame differences of the path delay difference sequence are calculated between adjacent frames with higher reliability. Inter-frame differences spanning low-reliability frames are filled in using interpolation estimation to ensure continuous difference results in the path delay difference sequence over time.

[0029] A time-varying slope is extracted from the path delay difference sequence to generate an estimate of the ionospheric drift rate. The time-varying slope is calculated window by window from the linear fitting slope of the path delay difference sequence within a sliding time window. The sliding window length is a fixed time span corresponding to the frame rate of the path delay difference sequence. Adjacent windows slide with a half-window step size to balance temporal resolution and the number of fitting samples. Path delay difference sequence entries within a window participate in least squares slope fitting with reliability level as the weight. Entries with higher reliability levels contribute more residual constraint force in the fitting. The time-varying slope is measured in terms of the rate of change of time delay (in seconds). The ionospheric drift rate estimate is converted into the rate of change of ionospheric virtual height by multiplying the time-varying slope by the speed of light and then dividing by the two-way coefficient. The conversion formula is v_h = (c / 2) × s, where v_h is the estimated ionospheric drift rate (in meters per second), c is the speed of light (in meters per second), s is the linear fitting slope of the path delay difference sequence within the current sliding window (in seconds per second), and the coefficient 2 comes from the two-way conversion of the electromagnetic wave's round-trip propagation path. When the time-varying slope smoothly transitions between adjacent windows of the path delay difference sequence, the estimated ionospheric drift rate shows a gradual trend. When the slope of the path delay difference sequence abruptly changes within a certain window, the estimated ionospheric drift rate shows a corresponding rate jump. The magnitude of the estimated ionospheric drift rate in this window is higher than that in adjacent windows, and the fitting residual of the slope abrupt window also increases accordingly. The fitting quality of this window is marked as low confidence to remind subsequent calculation steps to reduce the contribution of this window. The ionospheric drift rate estimate is calculated for each sliding window covered by the path delay difference sequence. Windows without sufficient reliable entries in the path delay difference sequence are filled with the ionospheric drift rate estimate of the preceding and following windows by linear interpolation.

[0030] Real-time ionospheric height deviation is generated by dynamically weighting the path delay difference sequence based on the ionospheric drift rate estimate. Windows with larger ionospheric drift rate estimates correspond to periods of rapid ionospheric height change, during which the absolute deviation of the path delay difference sequence is difficult to directly and stably estimate. The dynamic weighting calculation uses the reciprocal of the ionospheric drift rate estimate as the calculation weight for each window, with windows having smaller drift rate estimates having higher weights and windows having larger drift rate estimates having lower weights. Before assignment, the weight values ​​are normalized using the median of the drift rate estimate sequence to avoid excessively large or small weight values ​​for extreme drift rate windows. The absolute delay difference of the path delay difference sequence in low drift rate windows is directly mapped to the ionospheric height deviation benchmark value at that moment through conversion. The conversion relationship uses the propagation speed of electromagnetic waves in free space and the two-way path conversion factor to convert the delay difference to the height difference. The ionospheric height deviation of high drift rate windows is estimated by integral extrapolation of the baseline values ​​of adjacent low drift rate windows and the estimated ionospheric drift rate. The extrapolation starting point is selected from the low drift rate baseline window closest to the high drift rate window. The extrapolation step size is measured by the time interval between adjacent windows (in seconds). The product of the extrapolation step size and the estimated ionospheric drift rate (in m / s) gives the cumulative height deviation of each high drift rate window relative to the baseline (in meters, converted to km and then superimposed on the baseline value). The real-time ionospheric height deviation is summarized by using the window timestamp as an index for all windows. Real-time ionospheric height deviation entries are generated for all windows covered by the path delay difference sequence. The real-time ionospheric height deviation entries corresponding to the continuous low drift rate segments in the path delay difference sequence have the highest reliability. The reliability of the extrapolation entries for high drift rate segments decreases as the extrapolation step size increases. When the number of consecutive extrapolation steps exceeds the set upper limit, the real-time ionospheric height deviation entries for that window are marked as pending verification.

[0031] The corrected virtual height value is obtained by self-correcting the virtual height of the ionosphere based on the real-time ionospheric height deviation. The virtual height of the ionosphere is the nominal reference value of the equivalent height of ionospheric reflection, provided by the empirical ionospheric model or historical data from the vertical sounder. The nominal virtual height value reflects the statistical average level of ionospheric reflection height at a specific time and geographical location. The actual ionospheric reflection height is affected by factors such as diurnal variation, solar activity, and geomagnetic disturbances, and continuously deviates from the nominal virtual height value. The self-correction of the virtual height of the ionosphere superimposes the corresponding real-time ionospheric height deviation on a window-by-window basis. The formula for calculating the corrected virtual height value is H_c=H_0+ΔH(t), where H_c is the corrected virtual height value (unit: km), H_0 is the nominal virtual height value (unit: km), and ΔH(t) is the real-time ionospheric height deviation (unit: km) corresponding to the current observation window time t. Deviation entries marked as unverified are replaced by linear interpolation of adjacent verified window entries before superposition to avoid uncertain deviation values ​​affecting the reliability of the correction. The update frequency of the corrected false height values ​​is consistent with the window sliding step size of the real-time ionospheric height deviation, and each sliding window outputs one corrected false height value entry. The corrected false height value comes with an uncertainty range, which is jointly determined by the nominal false height model accuracy and the reliability level of the deviation entry. Entries with high reliability correspond to smaller uncertainty ranges for their corrected false height values, while those with low reliability have correspondingly larger uncertainty ranges. The corrected false height values ​​and their uncertainty ranges together provide real-time corrected ionospheric path parameter inputs for subsequent reflection slant range calculations.

[0032] In some embodiments, the step of combining the corrected virtual height value and the elevation angle estimate with the wave arrival angle parameter set to jointly calculate the reflection slant range and determine the candidate positioning coordinate set includes: extracting an azimuth probability density sequence based on the wave arrival angle parameter set; generating a multi-hypothesis slant range set based on the azimuth probability density sequence using the corrected virtual height value and the elevation angle estimate; performing probability weighting filtering on the multi-hypothesis slant range set according to the azimuth probability density sequence to obtain a weighted slant range distribution; and expanding the polar coordinate projection according to the weighted slant range distribution to determine the candidate positioning coordinate set.

[0033] The azimuth probability density sequence is extracted based on the incoming wave angle parameter set. Each azimuth entry in the compressed azimuth distribution of the incoming wave angle parameter set carries a weight value. The azimuth probability density sequence is generated by normalizing the weights of the azimuth entries in the incoming wave angle parameter set and mapping them onto a uniform azimuth grid. The normalization operation scales the sum of the weights of all azimuth entries to a unit value to ensure the consistency of the integral of the probability density sequence. The probability density value of each grid cell is given by the sum of the weights of the azimuth entries in the incoming wave angle parameter set falling within that cell. The grid cell width is half the median value of the spacing between azimuth entries in the incoming wave angle parameter set to ensure that adjacent entries do not merge into the same cell due to an excessively wide grid. Angle intervals within the azimuth entry set of the incoming wave angle parameter set exhibit high-density peaks in the azimuth probability density sequence, while sparse edge intervals correspond to low-density segments. The peak width of the azimuth probability density sequence directly reflects the range of azimuth uncertainty in the incoming wave angle parameter set. The grid resolution of the azimuth probability density sequence matches the azimuth entry density of the incoming wave angle parameter set. The effective interval width and main peak shape of the azimuth probability density sequence directly inherit the azimuth entry distribution characteristics of the incoming wave angle parameter set.

[0034] A multi-hypothesis slant range set is generated based on the azimuth probability density sequence, driven by the combined correction of the virtual height and the elevation estimate. Each effective grid cell in the azimuth probability density sequence corresponds to an azimuth hypothesis. Each azimuth hypothesis, combined with the elevation estimate and the correction of the virtual height, calculates the reflection slant range from the receiving station to the radiation source using a spherical geometric reflection model. The formula for calculating the reflection slant range is R = 2H_c / sin(θ_e), where R is the reflection slant range (unit: km), H_c is the correction of the virtual height (unit: km), and θ_e is the elevation estimate (unit: °). The coefficient 2 is derived from the two-way propagation path of the electromagnetic wave reflected by the ionosphere. The projection distance of the reflection slant range onto the horizontal plane is given by D = R × cos(θ_e) = 2H_c / tan(θ_e). The polar coordinate projection uses D as the polar radius, and θ_e is accompanied by a weighted standard deviation σ_θ as an uncertainty measure. Each azimuth hypothesis sets multiple elevation angle sampling values ​​at uniform intervals within the elevation uncertainty range [θ_e-σ_θ, θ_e+σ_θ]. The combination of the azimuth hypothesis and each elevation angle sampling value for each effective grid cell in the azimuth probability density sequence is substituted into the reflection slant range formula to calculate the corresponding slant range hypothesis value. The distance difference between entries in the multi-hypothesis slant range set originates from both the azimuth distribution of the azimuth probability density sequence and the sampling difference within the elevation uncertainty range. The multi-hypothesis slant range set aggregates the slant range hypothesis values ​​corresponding to all combinations of effective grid cells and elevation angle sampling values. The number of entries in the multi-hypothesis slant range set is equal to the product of the number of effective grid cells and the number of elevation angle sampling values. Each entry carries an azimuth label and a comprehensive screening weight for subsequent weighted screening and polar coordinate projection. The comprehensive screening weight is determined by normalizing the product of the probability density value of the corresponding grid cell in the azimuth probability density sequence and the comprehensive elevation weight corresponding to that elevation angle sampling value.

[0035] The multi-hypothesis slant range set is probabilistically weighted and filtered based on the azimuth probability density sequence and the comprehensive weight of elevation angle to obtain a weighted slant range distribution. For each surface distance hypothesis value in the multi-hypothesis slant range set, the product of the probability density value of its corresponding azimuth grid cell and the comprehensive weight of elevation angle sampling is used as the filtering weight for that hypothesis entry. The probabilistic weighted filtering arranges each hypothesis entry in the multi-hypothesis slant range set from high to low according to the filtering weight. When the cumulative weight exceeds a set threshold of the total weight, the low-weight hypothesis entries at the bottom of the list are removed. The set threshold is 90% of the total weight, meaning the cumulative weight of the retained portion covers 90% of the total weight of all hypotheses. The retained portion constitutes the weighted slant range distribution. Hypotheses corresponding to the main peak interval of the azimuth probability density sequence and the high comprehensive weight interval of elevation angle have higher filtering weights and are retained in a high proportion in the weighted slant range distribution. Hypotheses corresponding to the low-density intervals at the edge of the azimuth probability density sequence or the edge of the uncertainty range of elevation angle are largely removed after probabilistic weighted filtering. The range of entries in the weighted slant range distribution is relatively large, and the hypothesis slant range set narrows simultaneously in both the azimuth and range dimensions. The number of entries in the weighted slant distance distribution increases monotonically with the width of the main peak and the range of uncertainty in the elevation angle of the azimuth probability density sequence. After probability weighting, the range of entries in the weighted slant distance distribution narrows synchronously in both the azimuth and range dimensions, and low-weighted hypothetical entries have been removed from the distribution.

[0036] Candidate positioning coordinate sets are determined by polar coordinate projection using a weighted slant range distribution. Each entry in the weighted slant range distribution is projected onto a ground coordinate system using its assumed ground distance as the polar radius and the associated azimuth grid cell azimuth as the polar angle, with the receiving station coordinates as the pole, according to the polar coordinate relationship between ground distance and azimuth. The projection process considers the influence of the Earth's curvature on long-distance projection results; the difference between spherical and planar projections increases with distance. Each polar coordinate projection result is converted into geographic coordinates to form a single coordinate entry in the candidate positioning coordinate set. The coordinate density distribution of the candidate positioning coordinate set is consistent with the weight distribution of the entries in the weighted slant range distribution. In the weighted slant range distribution, entries corresponding to geographic coordinate regions with higher weights in the slant range and azimuth are denser, while entries corresponding to lower weights in the slant range distribution are sparser. Correction errors in the estimation of artificial height values ​​cause an overall shift in the weighted slant range distribution, leading to a systematic shift in the coordinate system of the candidate positioning coordinate set. Uncertainty in the elevation angle estimation value expands the slant range range of the weighted slant range distribution, thus expanding the coordinate coverage area of ​​the candidate positioning coordinate set. The candidate positioning coordinate set completes polar coordinate projection on all slant range and azimuth combinations covered by the weighted slant range distribution. Slant ranges that are eliminated by probability weighting in the weighted slant range distribution are assumed not to participate in the generation of the candidate positioning coordinate set. The weight value and azimuth label of each entry in the candidate positioning coordinate set are written along with the coordinates during generation as input fields for geometric consistency verification and weighted clustering.

[0037] Step S130: Identify multi-hop propagation components in the shortwave antenna array received data, extract ground reflection constraints for each hop, perform geometric consistency verification between each hop ground reflection constraint and the candidate positioning coordinate set to form a coordinate reliability grading table, and establish a slant range optimization configuration by weighting the candidate positioning coordinate set according to the coordinate reliability grading table.

[0038] In some embodiments, the step of identifying multi-hop propagation components and extracting ground reflection constraints for each hop in the shortwave antenna array received data includes: separating the shortwave antenna array received data into hop number components to generate a set of arrival waveforms for each hop; extracting the time delay difference between adjacent hops from the set of arrival waveforms to generate a sequence of time delay differences between hops; identifying the proportional deviation of the time delay difference sequence to generate a horizontal gradient correction amount; and performing gradient compensation on the position of each hop reflection point based on the horizontal gradient correction amount to generate ground reflection constraints for each hop.

[0039] Hop-count component separation is performed on the received data from the shortwave antenna array to generate waveform sets for each hop. After broadband correlation processing, the signals from each element of the received data form a time-delay power spectrum on the time-delay axis. The energy peaks on the time-delay power spectrum correspond to propagation components with one, two, and higher hop counts, in ascending order of time delay. Each hop component reaches the receiving station after undergoing a corresponding number of alternating reflections from the ionosphere and the ground. Effective energy peak identification uses a threshold of 6 dB exceeding the noise floor. The noise floor is estimated from the mean energy value in the tail region of the time-delay power spectrum received by the shortwave antenna array. Adjacent energy peaks are considered adjacent hop component pairs if the deviation between their time-delay interval and the theoretical single-hop time-delay increment is within 20%. The theoretical single-hop time-delay increment is calculated jointly from the great circle distance from the receiving station to the radiation source and the current elevation angle estimate. Isolated energy peaks with deviations exceeding the threshold are marked as suspected spurious peaks and do not participate in hop count allocation. After assigning each energy peak to a hop sequence number, a time-domain waveform is extracted from the signal of each element in the shortwave antenna array received data using a corresponding time delay window. The width of the extraction window is three times the time delay resolution unit of that hop to preserve the complete pulse shape of the hop component. The window boundary is weighted with a cosine conical window to suppress truncated sidelobe leakage. The extraction results of each element are stacked in the order of element arrangement to form the element waveform matrix of that hop. The element waveform matrices of all hop sequences are arranged in hop order to form the arrival waveform set of each hop. When a single hop entry is missing in the arrival waveform set of each hop, it is filled with a zero matrix to ensure continuous alignment of the hop sequence index. When the resolvability condition is not met, the effective number of hops in the arrival waveform set of each hop is automatically capped at the hop sequence corresponding to the last resolvable peak. The effective number of hops in the arrival waveform set of each hop determines the upper limit of the number of entries in the subsequent inter-hop time delay difference sequence.

[0040] For each hop arriving at the waveform set, the inter-hop delay difference is extracted to generate an inter-hop delay difference sequence. The peak delay positions of the array element waveform matrices of adjacent hop sequences in each hop arriving at the waveform set are extracted. The peak delay positions are refined to sub-sampling precision using quadratic interpolation of the peak values ​​of discrete sampling points. The difference in peak delay between adjacent hop pairs is the single-pair inter-hop delay difference. The single-pair inter-hop delay differences of all array elements are weighted and averaged using the signal-to-noise ratio as the weight to obtain the representative delay difference value for that hop pair. The extraction of representative delay differences starts from the first and second hop pairs in each hop arriving at the waveform set and proceeds sequentially along the hop sequence to the last complete hop pair. The representative delay difference values ​​for each hop pair are arranged according to the hop pair number to form the inter-hop delay difference sequence. The number of entries in the inter-hop delay difference sequence is equal to the effective hop number of each hop arriving at the waveform set minus one. For example, a three-hop waveform set generates two inter-hop delay difference entries. When the actual ionospheric horizontal gradient exists, the values ​​of the inter-hop delay difference sequence entries tend to deviate from the proportional benchmark along the hop order. The direction of deviation corresponds to the gradient direction. For example, as the virtual ionospheric height gradually increases along the signal propagation direction, each entry in the inter-hop delay difference sequence monotonically increases with the hop order. Frames with extremely low signal-to-noise ratios have large delay peak extraction errors, resulting in low reliability of the inter-hop delay difference sequence entries. The inter-hop delay difference sequence uses a signal-to-noise ratio threshold to label the reliability level of each frame entry, with low-reliability entries participating in subsequent proportional deviation identification in a weighted manner. When a hop arrives at a waveform set and a hop is filled with a zero matrix, the hop pair delay difference adjacent to that hop is marked as invalid. Invalid entries in the inter-hop delay difference sequence are supplemented by linear interpolation of the preceding and following valid entries to ensure the continuity and completeness of the hop pair index in the inter-hop delay difference sequence.

[0041] For example, the step of identifying the proportional deviation of the inter-hop delay difference sequence to generate a horizontal gradient correction includes: constructing a hop-by-hop delay difference ratio sequence based on adjacent hop pairs of the inter-hop delay difference sequence; extracting a median reference value from the hop-by-hop delay difference ratio sequence to generate a proportional reference sequence; fitting a linear trend between the hop-by-hop delay difference ratio sequence and the proportional reference sequence along the hop order direction to extract the trend slope and trend sign; and mapping the ionospheric horizontal gradient amplitude and gradient direction based on the trend slope and trend sign to generate a horizontal gradient correction.

[0042] A hop-by-hop delay difference ratio sequence is constructed from adjacent hop pairs for the inter-hop delay difference sequence. The ratio of the (i+1)th entry to the ith entry in the inter-hop delay difference sequence is defined as the ith hop ratio. Each hop ratio is derived by iterating through all adjacent entry pairs in the inter-hop delay difference sequence, and these ratios are arranged according to the hop pair number to form the hop-by-hop delay difference ratio sequence. Under a uniform horizontal ionosphere, the hop delay differences are equal, and the ratios of all entries in the inter-hop delay difference sequence approach 1. Entries in the hop-by-hop delay difference ratio sequence that deviate from 1 reflect unequal changes in the path lengths of adjacent hop pairs. The source of these unequal changes is the horizontal gradient of the ionospheric virtual height within the propagation path coverage area. The length of the hop-by-hop delay difference ratio sequence is equal to the number of entries in the inter-hop delay difference sequence minus one. For example, a three-hop waveform set generates a two-entry inter-hop delay difference sequence, corresponding to a one-entry hop-by-hop delay difference ratio sequence; a four-hop waveform set generates a three-entry inter-hop delay difference sequence, corresponding to a two-entry hop-by-hop delay difference ratio sequence. When the inter-hop delay difference sequence contains invalid entries, the hop ratio adjacent to the invalid entry is marked as unreliable. Unreliable entries in the hop-by-hop delay difference ratio sequence are excluded in the subsequent median benchmark extraction. An overall hop-by-hop delay difference ratio sequence greater than 1 indicates that the ionospheric virtual height gradually increases along the propagation direction, while an overall value less than 1 indicates that the ionospheric virtual height gradually decreases along the propagation direction. The absolute amount of the ratio deviating from 1 is monotonically correlated with the amplitude of the virtual height gradient within the corresponding segment. When the shortwave propagation path crosses a region with a significant ionospheric gradient, the equivalent reflection height experienced by each hop changes sequentially, and the hop delay differences of adjacent hops therefore lose their proportional relationship. The systematic deviation of the hop-by-hop delay difference ratio sequence directly quantifies the ionospheric tilt amplitude of the area covered by the propagation path. The horizontal gradient correction amount is used accordingly to apply differential compensation to the position of each hop reflection point to restore the true reflection point distribution.

[0043] A proportional reference sequence is generated by extracting the median reference value from the hop-by-hop delay difference ratio sequence. The median of each valid entry in the hop-by-hop delay difference ratio sequence is used as the median reference value. The median is insensitive to extreme deviations in the hop-by-hop delay difference ratio sequence and can robustly reflect the reference ratio level of the ionospheric gradient within the propagation path coverage area. The median reference value is expanded and filled along the hop pair numbers in the form of a constant sequence to form a sequence of the same length as the hop-by-hop delay difference ratio sequence. This constant sequence is defined as the proportional reference sequence, which represents the reference level that the hop-by-hop delay difference ratio sequence should present when there is no trend gradient drift. When the hop-by-hop delay difference ratio sequence has only one valid entry, the median reference value is the entry itself, and the proportional reference sequence degenerates into a single-valued constant. Subsequent linear trend fitting uses this constant as the reference for zero slope estimation. When all entries in the hop-by-hop delay difference ratio sequence are unreliable, the proportional reference sequence is assigned a sequence of all 1s, corresponding to the homogeneous ionospheric assumption. Subsequent trend slope and trend sign are also assigned zero values. The larger the deviation of the constant value of the equidistant reference sequence from 1, the more systematic the non-uniformity of the ionosphere in the area covered by the propagation path of this batch is. The direction of deviation indicates the overall trend of the virtual height changing with the propagation distance.

[0044] The trend slope and trend sign are extracted by fitting a linear trend between the hop-by-hop delay difference ratio sequence and the equidistant reference sequence along the hop order. The item-by-item differences between the hop-by-hop delay difference ratio sequence and the equidistant reference sequence form a residual sequence. The residual sequence is then fitted with a least-squares linear relationship using the hop pair number as the independent variable. The slope of the fitted line is the trend slope, and the sign of the slope indicates the trend sign. The dimension of the trend slope is the change in ratio per hop pair. A large absolute value of the trend slope indicates that the deviation of the hop-by-hop delay difference ratio sequence from the equidistant reference sequence increases rapidly with the hop order, corresponding to a significant horizontal gradient of the ionosphere along the propagation path. A trend slope close to zero indicates that the equidistant relationship of the delay differences between each hop is stable, corresponding to an insignificant horizontal gradient of the ionosphere along the propagation path. A positive trend sign indicates that the hop-by-hop delay difference ratio sequence increases relative to the equidistant reference sequence with the hop order, indicating that the artificial height of the ionosphere at the far end of the propagation path is higher than that at the near end; a negative trend sign indicates that the artificial height at the far end is lower than that at the near end. When the residual sequence has only one valid entry, the trend slope is approximated by the residual value at that point as the change within the step size of a single jump pair, and the trend sign is taken as the positive or negative sign of the residual at that point. When the root mean square residual of the linear fit is greater than 0.1, the fitting quality is labeled as low, and the reliability of the trend slope and trend sign is reduced. The subsequent horizontal gradient correction is accompanied by uncertainty labeling. The uncertainty value is characterized by the ratio of the root mean square residual to the product of the absolute value of the trend slope and the number of valid entries in the residual sequence. This ratio is a dimensionless quantity, reflecting the relative uncertainty of the residual relative to the trend amplitude.

[0045] The horizontal gradient correction is generated by mapping the ionospheric horizontal gradient amplitude and direction based on the trend slope and trend sign. The trend slope is mapped to the ionospheric virtual height horizontal gradient amplitude through the ionospheric horizontal gradient geometric propagation model. The mapping formula is g = s × τ0 × c / (2 × D_hop), where g is the ionospheric horizontal gradient amplitude (unit: km / km), s is the trend slope (unit: ratio / hop pair), τ0 is the median time delay difference of the inter-hop time delay difference sequence (unit: s), c is the speed of light (unit: km / s), D_hop is the single-hop ground propagation distance (unit: km), and the coefficient 2 comes from the electromagnetic wave two-way propagation conversion. D_hop is jointly calculated from the receiver coordinates, the azimuth angle of the arrival angle parameter set, and the corrected virtual height value. The ground propagation distance of each hop is approximately equal under the assumption of a uniform ionosphere. A positive trend sign indicates that the gradient direction points to the far end along the propagation path, corresponding to a greater virtual height of the ionosphere at the far end than at the near end; a negative trend sign indicates that the gradient direction points to the near end, corresponding to a greater virtual height of the ionosphere at the near end than at the far end. The magnitude and direction of the ionospheric horizontal gradient are jointly encapsulated into a horizontal gradient correction. The magnitude component of the horizontal gradient correction gives the virtual height change per unit horizontal distance for each bounce point, while the direction component gives the gradient ascent direction in azimuth angles. When the absolute value of the trend slope is less than 0.02, the magnitude of the ionospheric horizontal gradient is determined to be negligible, the magnitude of the horizontal gradient correction is assigned zero, the gradient direction is marked as uncertain, and the positions of each bounce point remain unchanged from the initial estimate under the assumption of a homogeneous ionosphere.

[0046] Gradient compensation is applied to the positions of each hop reflection point based on the horizontal gradient correction to generate ground reflection constraints for each hop. Under the assumption of a uniform ionosphere, the initial position of each reflection point within each hop ground reflection constraint is determined by spherical geometric transformation using the receiver coordinates, hop delay values, and elevation angle estimates. The initial positions are uniformly distributed along the propagation azimuth direction according to the hop order, with the distance between adjacent initial reflection points equal to D_hop. The horizontal gradient correction provides two components: gradient magnitude and gradient direction. The gradient compensation amount for each hop reflection point within each hop ground reflection constraint is determined by the product of the horizontal distance from the reflection point to the receiver and the gradient magnitude of the horizontal gradient correction. The farther the hop reflection point is from the receiver, the greater the cumulative compensation amount. The component of the gradient compensation amount perpendicular to the propagation azimuth causes a lateral offset of each reflection point within each hop ground reflection constraint, while the component along the propagation azimuth direction causes a longitudinal offset. The longitudinal offset causes the reflection point to move forward or backward along the propagation path, while the lateral offset causes the reflection point to deviate from the original propagation great circle. The superposition of the two offset vector components gives the compensated position of each hop reflection point. The set of compensated reflection point locations for all jump sequences is encapsulated as ground reflection constraints for each jump. Each entry in each ground reflection constraint is associated with the corresponding jump sequence number, the compensated coordinates, and the uncertainty range of the compensation amount. The uncertainty range is jointly determined by the uncertainty associated with the horizontal gradient correction amount and the distance from the jump to the receiving station. When the horizontal gradient correction amount is zero, each entry in each ground reflection constraint for each jump maintains its initial position under the assumption of a uniform ionosphere.

[0047] In some embodiments, the step of performing geometric consistency verification between the ground reflection constraints of each jump and the candidate positioning coordinate set to form a coordinate reliability grading table includes: arranging the ground reflection constraints of each jump in order of jump number to establish a multi-hop geometric sequence; calculating the radiation source position constraints of each jump by extending backward in the direction of the radiation source based on the multi-hop geometric sequence to generate a radiation source position constraint set; performing spatial overlap quantification on the radiation source position constraint set and the candidate positioning coordinate set to generate an overlap sequence; and classifying the coordinate reliability grading table according to the overlap sequence.

[0048] A multi-hop geometric sequence is established by arranging the ground reflection constraints in hop order. Each hop reflection point in the ground reflection constraint carries a hop sequence number, compensated coordinates, and uncertainty range. Arranged in ascending order of hop sequence number, these form an ordered chain of reflection points extending from the receiving station towards the radiation source; this ordered chain constitutes the multi-hop geometric sequence. The multi-hop geometric sequence starts at the receiving station coordinates, with the first hop reflection point as the first node. Subsequent hop reflection points are added sequentially, ending in the direction of the calculated radiation source. The distance between adjacent nodes corresponds to the horizontal propagation distance of a single hop. Under uniform ionospheric conditions, the distance between adjacent nodes in the multi-hop geometric sequence is approximately equal. After horizontal gradient compensation, the non-uniformity of the node distances reflects the actual tilt distribution of the ionosphere within the propagation path coverage area. The spatial error ellipse of the nodes in the multi-hop geometric sequence expands cumulatively with the hop sequence, with the largest uncertainty range at the terminal node. The accuracy of the radiation source location constraint calculation increases with the number of nodes in the multi-hop geometric sequence because more nodes provide more independent geometric constraints. When the multi-hop geometric sequence contains only one valid node, it degenerates into a single-hop geometric sequence, and the coverage area of ​​the subsequent radiation source location constraint set expands accordingly.

[0049] For example, the step of generating a radiation source position constraint set by extending the multi-hop geometric sequence backwards towards the radiation source direction to calculate the position constraints of each hop radiation source includes: extracting the incident elevation angle estimate of each hop reflection point in the multi-hop geometric sequence to generate a set of hop incident elevation angles; calculating the upper and lower bounds of the radiation source distances by extending each hop reflection point in the multi-hop geometric sequence along the backward propagation path based on the set of hop incident elevation angles to generate a set of hop radiation source distance intervals; weighting and overlapping the set of hop radiation source distance intervals according to the number of hops to generate a fused radiation source distance interval; and expanding the spatial projection based on the fused radiation source distance interval and the azimuth information in the multi-hop geometric sequence to generate a set of radiation source position constraints.

[0050] An estimated incident elevation angle is extracted from each hop reflection point in the multi-hop geometric sequence to generate a set of hop incident elevation angles. The incident elevation angle of each hop reflection point in the multi-hop geometric sequence is determined by the horizontal distance from the reflection point to the previous node and the corresponding corrected dummy height value. The node above the first hop reflection point is the receiving station, and the node above subsequent hop reflection points is the previous hop reflection point. The formula for calculating the incident elevation angle estimate is θ_i = arctan(H_i / L_i), where θ_i is the estimated incident elevation angle of the i-th hop reflection point (unit: °), H_i is the corrected dummy height value corresponding to the i-th hop (unit: km), and L_i is the horizontal great circle distance between the i-th hop reflection point and its previous node (unit: km). The corrected dummy height value H_i for each hop is obtained by hop-by-hop superimposing an offset on the baseline corrected dummy height value based on the gradient magnitude when the horizontal gradient correction is not zero. The offset is equal to the product of the horizontal distance of the hop reflection point and the gradient magnitude. Each jump incidence elevation angle estimate is arranged in jump order and set as a jump incidence elevation angle set. The number of entries in each jump incidence elevation angle set is equal to the number of effective nodes in the multi-hop geometric sequence. Jump reflection points with larger horizontal gradient compensation in the multi-hop geometric sequence correspond to larger uncertainties in the corrected virtual height value, and the uncertainty range of the jump incidence elevation angle estimate is correspondingly expanded. Entries with larger uncertainty ranges in each jump incidence elevation angle set contribute a wider interval in the subsequent radiation source distance range estimation. The trend of elevation angle values ​​along the jump order reflects the systematic change of the virtual ionospheric height within the propagation path coverage area. An increase in elevation angle value with jump order indicates that the virtual ionospheric height increases along the propagation direction.

[0051] Based on the set of incident elevation angles for each jump, the distance interval set of each jump reflection point in the multi-hop geometric sequence is generated by extending along the reverse propagation path. Each jump incident elevation angle estimate in the set of incident elevation angles for each jump is accompanied by an uncertainty range. The upper and lower bounds of the elevation angle within the uncertainty range correspond to the shortest and longest ground propagation distances from the radiation source to the jump reflection point, respectively. The calculation formula is d_i=H_i / tan(θ_i), where d_i is the horizontal distance (in km) from the radiation source to the reflection point calculated for the i-th jump. The lower bound of the distance is determined by substituting the upper bound of the elevation angle, and the upper bound of the distance is determined by substituting the lower bound of the elevation angle. Starting from each hop reflection point, the distance is extended along the reverse propagation azimuth. The reverse azimuth is determined by the opposite direction of the line connecting the hop node and the previous node in the multi-hop geometric sequence. Each hop landing interval is marked as a line segment with a distance range width in the ground coordinate system. The radial width of the line segment is equal to the difference between the upper and lower distance bounds. The set of line segments for all hop sequences is encapsulated as the distance interval set for each hop radiation source. Each entry carries the hop sequence number, upper distance bound, lower distance bound, and reverse azimuth. The interval width of high-hop sequence entries in the distance interval set for each hop radiation source is usually larger than that of low-hop sequence entries. This rule stems from the fact that high-hop sequence reflection points are farther from the receiving station, corresponding to a wider range of uncertainty in elevation angle estimation, and the difference between the upper and lower distance bounds increases with the hop sequence. Low-hop sequence reflection points are closer to the receiving station and have better elevation angle measurement conditions, corresponding to a narrower radiation source distance interval. The monotonic relationship between interval width and hop sequence makes the reliability of each hop radiation source distance interval set entry decrease with the hop sequence. The constraint ability of low-hop sequence entries on the radiation source location is significantly stronger than that of high-hop sequence entries.

[0052] The fused radiation source distance interval is generated by weighting and overlapping the hop distance intervals of each hop set according to the number of hops. The reliability of each hop distance interval in the hop set decreases with the hop order. The weighted overlap uses the reciprocal of the hop order as the weight coefficient for each hop distance interval, with lower hop order intervals having higher weights. The weighted overlap takes the weighted intersection of the radiation source distance intervals of all valid hops in each hop distance interval set. The lower bound of the weighted intersection is determined by taking the high-weight quantile after weighting and sorting the lower bound values ​​of each hop according to the weight coefficients of each hop. The upper bound is determined by taking the low-weight quantile after weighting and sorting the upper bound values ​​of each hop according to the weight coefficients of each hop. The position of the quantile is determined by the cumulative weight after normalization of the weight coefficients of each hop. When there is a non-empty common overlap among the radiation source distance intervals of all hops in the hop distance interval set, the overlap and narrowing result is the fused radiation source distance interval. The width of the fused radiation source distance interval narrows as the number of valid hops participating in the overlap increases. The three-hop overlap result is usually more than 50% narrower than the width of a single-hop interval. When there is no common overlap among the range intervals of each hop radiation source, the range interval of the hop with the highest weight is merged and a warning mark for inconsistency among multiple hops is attached. This status indicates that there is abnormal scattering or measurement error in the propagation path, and the uncertainty settings of each item in the range of incident elevation angles of each hop need to be reviewed again.

[0053] A spatial projection is used to generate a radiation source position constraint set based on the fused radiation source distance interval and the azimuth information of the multi-hop geometric sequence. The coordinates of each node in the multi-hop geometric sequence determine the azimuth reference direction of the propagation path. The azimuth extending backward from each hop reflection point is obtained by reversing the vectors between nodes in the multi-hop geometric sequence. The dispersion of the reverse azimuth of each hop reflects the curvature of the propagation link; the greater the curvature, the more significant the difference between the reverse azimuths of each hop. The fused radiation source distance interval is extended backward from each hop reflection point, forming an arc segment with a distance range width in the ground coordinate system. The radial width of the arc segment is determined by the upper and lower bounds of the fused radiation source distance interval, and the angular width is determined by the dispersion of the reverse azimuth of each hop. The arc segments of all hop sequences are superimposed in the ground coordinate system, and the intersection area is the spatial coverage of the radiation source position constraint set. The area of ​​the radiation source position constraint set decreases significantly with the increase of the effective hop count. The smaller the uncertainty of the azimuth angle of each node in the multi-hop geometric sequence, the narrower the angle width of each hop arc segment, and the smaller the horizontal width of the radiation source position constraint set; the narrower the distance interval of the fused radiation source, the smaller the vertical width of the radiation source position constraint set. The radiation source position constraint set is encapsulated by a set of polygon vertex coordinates, and the three statistical measures of region center coordinates, vertical range, and horizontal range are output along with the encapsulation structure as input parameters for spatial overlap quantification.

[0054] The spatial overlap metric is used to generate an overlap sequence by quantifying the spatial overlap between the radiation source location constraint set and the candidate positioning coordinate set. The radiation source location constraint set is represented as a polygonal region with joint azimuth and distance constraints on the ground coordinate system. Each coordinate entry in the candidate positioning coordinate set is distributed as a point, falling at different locations inside and outside this polygonal region. The spatial overlap metric is based on each coordinate entry in the candidate positioning coordinate set, calculating the nearest distance between each coordinate entry and the boundary of the radiation source location constraint set. A negative distance value indicates an inclusion relationship when the coordinate entry falls inside the radiation source location constraint set, and a positive value indicates an inclusion relationship when it falls outside. The overlap value of each coordinate entry is mapped to the interval 0 to 1 from the nearest distance using an exponential decay function. The mapping formula is O_p = 1 / (1 + exp(r_p / λ_d)), where O_p is the overlap value of coordinate entry p (dimensionless, ranging from 0 to 1), r_p is the nearest distance between coordinate entry p and the boundary of the radiation source location constraint set (unit: km, negative when it falls inside the constraint set, positive when it falls outside), and λ_d is the decay constant (unit: km), λ_d = β × R_eq, where R_eq is the equivalent half of the radiation source location constraint set. The radius (unit: km) and β are proportionality coefficients. A larger equivalent radius results in a larger attenuation constant and a wider boundary transition zone. When multi-hop constraints are insufficient, the radiation source location constraint set covers a larger area. The wider boundary transition zone allows candidate positioning coordinate entries near the constraint set boundary to still obtain continuously varying moderate overlap values ​​instead of being rigidly truncated, ensuring the continuity of the reliability distribution even with limited positioning constraint accuracy. When multi-hop constraints are sufficient, the radiation source location constraint set narrows significantly, the attenuation constant decreases accordingly, and the overlap drops sharply at the boundary, efficiently distinguishing entries inside and outside the constraint set. Coordinate entries falling within the core region of the radiation source location constraint set have an overlap value close to 1, entries near the boundary have an overlap value of approximately 0.5, and entries outside the boundary rapidly approach 0 with increasing deviation distance. The overlap values ​​of each coordinate entry in the candidate positioning coordinate set are arranged by entry number to form an overlap sequence. This overlap sequence corresponds one-to-one with the candidate positioning coordinate set. The numerical distribution of the overlap sequence directly reflects the screening effectiveness of multi-hop geometric constraints on the candidate positioning coordinate set. When the constraint set narrows, high-value entries concentrate in the core region.

[0055] A coordinate reliability grading table is formed by classifying coordinates according to their overlap sequence. The overlap sequence is divided into three levels based on numerical range: entries with an overlap value greater than 0.7 are classified as high reliability, those between 0.3 and 0.7 as medium reliability, and those less than 0.3 as low reliability. The grading threshold is adaptively adjusted based on the number of effective hops. When the number of effective hops is one, the lower limit for high reliability is lowered to 0.6, and the upper limit for low reliability is raised to 0.4 to avoid too few reliable entries. When the number of effective hops reaches three or more, the grading threshold returns to its default value. Each coordinate entry in the candidate positioning coordinate set is labeled with a reliability level according to the corresponding grading result in the overlap sequence: high reliability with a weight coefficient of 1.0, medium reliability with a weight coefficient of 0.5, and low reliability with a weight coefficient of 0.1. The coordinate value, reliability level, and weight coefficient of each coordinate entry are combined to form a single record in the coordinate reliability grading table. The set of records for all coordinate entries constitutes the complete coordinate reliability grading table. When the number of high-confidence entries in the overlap sequence is less than 3, the coordinate confidence level table is marked with an insufficient overlap warning. Subsequent slope distance optimization configurations are then supplemented with medium-confidence entries for weighting. The hierarchical structure of the coordinate confidence level table determines the narrowing degree of the slope distance optimization configuration. When the overlap sequence is concentrated at high values, the proportion of high-confidence entries in the coordinate confidence level table is large, correspondingly narrowing the slope distance optimization range.

[0056] The candidate positioning coordinate set is weighted and optimized using a coordinate reliability grading table to establish the optimal slant range configuration. The coordinate reliability grading table assigns a weight to each coordinate entry in the candidate positioning coordinate set: a weight coefficient of 1.0 for high reliability, 0.5 for medium reliability, and 0.1 for low reliability. This differentiated weighting causes the candidate positioning coordinate set to concentrate in the high reliability region during the weighted optimization. The reflection slant range corresponding to each coordinate entry is jointly determined by the azimuth angle and the corrected dummy height value from the incoming wave angle parameter set. The weights of the coordinate reliability grading table are mapped to the slant range values ​​of each coordinate entry to form a weighted slant range candidate set. High-weight slant range entries have a higher retention priority than low-weight entries. After the slant range candidate set is sorted in descending order of weight, low-weight tail entries whose cumulative weight exceeds a set retention ratio threshold are removed. The retention ratio threshold is 90% of the total weight; that is, retention stops when the cumulative weight reaches 90% of the total weight. The retained portion is grouped by azimuth angle intervals. The slant range range within each group and the weight distribution within each group together constitute the grouped entries for the optimal slant range configuration. The more concentrated the coordinate confidence level table is in the high confidence level, the narrower the range of the optimal slope distance configuration. When the high confidence level entries are concentrated in the azimuth angle range of 180° to 185°, the slope distance range of the corresponding group of the optimal slope distance configuration narrows to the concentrated distribution band of slope distance in this range. The optimal slope distance configuration encapsulates the slope distance range and weight distribution of each group using azimuth angle grouping as an index. The optimal slope distance configuration is encapsulated using azimuth angle grouping as an index, and the grouped slope distance range and weight distribution have established a traceable association with the coordinate confidence level table.

[0057] Step S140: Based on the coordinate confidence level table and the slant range optimization configuration, the candidate positioning coordinate set is collaboratively weighted and clustered to construct a multidimensional positioning parameter set. Based on the multidimensional positioning parameter set, the radiation source location fusion calculation is performed to determine the final positioning coordinates.

[0058] In some embodiments, the step of constructing a multidimensional positioning parameter set by collaborative weighted clustering of the candidate positioning coordinate set based on the coordinate confidence level table and the preferred slope distance configuration includes: weighted filtering and clustering of the candidate positioning coordinate set based on the coordinate confidence level table and the preferred slope distance configuration to generate a clustering residual sequence; performing reverse correction on the coordinate confidence level table based on the clustering residual sequence to generate an updated confidence level table; re-weighting and clustering the candidate positioning coordinate set according to the updated confidence level table to generate convergent positioning coordinates; and jointly encapsulating the convergent positioning coordinates with the updated confidence level table and the preferred slope distance configuration to construct a multidimensional positioning parameter set.

[0059] Based on the coordinate confidence level table and the optimal slant range configuration, the candidate positioning coordinate set is weighted and clustered to generate a clustering residual sequence. In the coordinate confidence level table, high-confidence coordinate entries participate in weighted clustering with a weight coefficient of 1.0, medium-confidence entries with 0.5, and low-confidence entries with 0.1. This differentiated weighting causes the cluster centroids to preferentially gravitate towards high-confidence coordinate regions. The optimal slant range configuration defines permissible slant ranges for each group of coordinate entries based on azimuth angle. Coordinate entries in the candidate positioning coordinate set whose slant range values ​​exceed the permissible range of their corresponding azimuth angle group are removed before weighted clustering. After screening, the number of participating entries in the candidate positioning coordinate set is reduced compared to before screening, thus decreasing the clustering computation scale. Weighted clustering performs weighted k-means iteration on the selected candidate location coordinate set. The number of clusters, k, is automatically determined by density gap scanning of the candidate location coordinate set. The density gap scan uses the number of significant breakpoints in the spacing distribution of coordinate entries plus one as the value of k. When the candidate location coordinate set is highly concentrated, k is 1; when there are clearly spatially separated clusters of entries, k increases accordingly. The initial cluster centers are the weighted centroids of each region in the k-equal division of the high-confidence entries in the coordinate confidence grading table. During the iteration process, each coordinate entry participates in the centroid update with its grading weight. The iteration converges when the centroid displacement of adjacent steps is less than the resolution threshold. The Euclidean distance between each coordinate entry and its corresponding cluster center constitutes the residual value of that entry. The residual values ​​of all participating entries are arranged by entry number to form a cluster residual sequence. The distribution pattern of the cluster residual sequence is jointly determined by the narrowing degree of the slant distance optimization configuration and the proportion of high-confidence entries in the coordinate confidence grading table. Both factors work together to contribute to the stability of the cluster centroids.

[0060] The coordinate confidence level table is updated by performing a reverse correction based on the clustering residual sequence. The residual value of each entry in the clustering residual sequence reflects the geometric deviation of the corresponding coordinate entry from the cluster center. Entries with large residuals indicate that the initial confidence level table overestimated geometric consistency, while entries with small residuals indicate that the initial level was either accurate or underestimated. The reverse correction downgrades entries in the coordinate confidence level table with residuals exceeding a set upper limit by one level, and upgrades entries with residuals below a set lower limit by one level. The upgrade magnitude does not exceed the boundary between adjacent levels. Entries already at a low confidence level in the coordinate confidence level table are not downgraded further. The upper and lower limits of the residuals are adaptively determined based on the median residual of the clustered residual sequence. 1.5 times the median residual is used as the down-adjustment trigger threshold, and 0.5 times the median residual is used as the up-adjustment trigger threshold. The thresholds are dynamically adjusted according to the actual dispersion of the clustered residual sequence entries. When ionospheric disturbances are strong, the overall dispersion of the clustered residual sequence is large, and the trigger thresholds are proportionally expanded to avoid misjudgments. When the ionosphere is stable, the trigger thresholds are narrowed, enabling precise differentiation of slightly deviated entries and timely correction of their confidence levels. After reverse correction, the coordinate confidence level table updates the level and weight coefficients of each entry synchronously, solidifying into an updated confidence level table. The updated confidence level table maintains the same three-level structure as the coordinate confidence level table, but the level assignment of each entry has been redistributed based on the residual evidence of the clustered residual sequence. The number of originally high-confidence entries in the updated confidence level table that were downgraded due to excessive residuals reflects the degree of deviation in the initial classification. The level distribution of the updated confidence level table is closer to the actual geometric consistency structure of the candidate positioning coordinate set than that of the coordinate confidence level table.

[0061] The candidate location coordinate set is re-weighted and clustered according to the updated confidence level table to generate convergent location coordinates. The updated confidence level table assigns updated weight coefficients to each entry in the candidate location coordinate set. Re-weighted clustering uses the updated weights from the updated confidence level table to replace the initial weights, driving the k-means iteration of the candidate location coordinate set. The initial cluster center is taken as the convergence centroid of the first round of clustering to accelerate convergence. During re-weighted clustering, the high-weight entries in the updated confidence level table exert a stronger pull on the cluster centroids than in the first round. Low-consistency entries that were overestimated in the initial grading contribute less to the centroids due to weight reduction, causing a shift in the cluster centroids compared to the first round. This shift reflects the correction magnitude of the confidence level table's reverse adjustment on the location results. After re-weighted clustering convergence, the centroid coordinates of the master cluster are defined as the convergent location coordinates. The master cluster is determined by the cluster with the largest number of weighted entries in the candidate location coordinate set. When the candidate location coordinate set is highly concentrated, the master cluster covers the vast majority of valid entries. The representativeness of the convergent location coordinates increases with the weight ratio of the master cluster entries. The convergent positioning coordinates, along with the weighted covariance matrix of the main cluster, serve as a description of positioning uncertainty. The major axis of the covariance matrix corresponds to the direction of maximum expansion of positioning error in the convergent positioning coordinates, which is usually consistent with the direction of the azimuth of the incoming wave. The minor axis corresponds to the direction of minimum error. The major and minor axis directions of the covariance matrix and the area of ​​the ellipse are determined by the weighted coordinate deviation of each entry in the main cluster. The major axis of the ellipse is usually consistent with the direction of the azimuth of the incoming wave.

[0062] A multidimensional positioning parameter set is constructed by jointly encapsulating converged positioning coordinates, an updated confidence level table, and a slant range optimization configuration. The converged positioning coordinates provide the optimal estimate of the radiation source location. The updated confidence level table records the confidence structure of each entry in the candidate positioning coordinate set after two rounds of weighted clustering iterations. The slant range optimization configuration stores the permissible slant range range and corresponding weight distribution for azimuth groupings. The multidimensional positioning parameter set, jointly encapsulated with the converged positioning coordinates, the updated confidence level table, and the slant range optimization configuration, simultaneously carries three types of positioning information: point estimation, uncertainty description, and slant range constraints. The multidimensional positioning parameter set uses the converged positioning coordinates as the main field, with associated fields storing the level distribution statistics of the updated confidence level table, the azimuth grouping boundaries and slant range ranges of the slant range optimization configuration, and the weighted covariance matrix of the converged positioning coordinates. The proportion of high-confidence entries in the updated confidence level table is recorded as a location confidence index in the multidimensional positioning parameter set. A high confidence index indicates strong reliability in the subsequent fusion calculation driven by the multidimensional positioning parameter set, while a low confidence index assigns a lower fusion weight to the multidimensional positioning parameter set. The encapsulation structure of the multidimensional positioning parameter set supports parallel input of positioning results from multiple time periods. Each multidimensional positioning parameter set entry carries a timestamp and a receiving station identifier. Subsequent radiation source location fusion calculations can align the multidimensional positioning parameter sets from different time periods based on the timestamps for joint processing.

[0063] The final positioning coordinates are determined by performing radiation source location fusion calculation based on a multi-dimensional positioning parameter set. The converged positioning coordinates from the multi-dimensional positioning parameter set serve as the initial input for the fusion calculation. The confidence index of the updated confidence level table is used as the overall fusion weight, and the slant range of the optimal slant range configuration is used as the slant range allowable boundary constraint. The fusion calculation performs weighted centroid merging on the converged positioning coordinates generated from observations in each time period, using the confidence index as the weight. For single-time period inputs, the converged positioning coordinates are directly output; for multi-time period inputs, the coordinates are merged into a unified geographic coordinate system using the confidence index of the multi-dimensional positioning parameter set for each time period as the weight. The initial coordinates after weighted centroid merging are constrained by the slant range allowable range of each azimuth group in the optimal slant range configuration as the boundary for projection. If the coordinates fall outside the allowable range, they are pulled back to the nearest allowable boundary. The pull-back distance is recorded as the slant range constraint correction. A large correction triggers confidence reduction processing, decreasing the overall contribution of this multi-dimensional positioning parameter set in the fusion. The coordinates after slant-range constraint projection are corrected using a maximum a posteriori (MAS) estimation, which is performed by combining the weighted covariance matrix. The correction amount is weighted by the inverse of the covariance matrix as the accuracy matrix. The MAS-corrected coordinates are the final positioning coordinates, with their source indicated by the receiving station identifier and observation period. The final positioning coordinates are accompanied by a comprehensive positioning uncertainty ellipse. The uncertainty ellipse is derived by merging the weighted covariance matrices of each participating multidimensional positioning parameter set. Its area shrinks as the sum of the confidence indices increases. It is output in a structured format along with the final positioning coordinates. The structure contains three types of fields: geographic coordinates, the major and minor axes of the uncertainty ellipse, and the azimuth of the major axis.

[0064] To implement the shortwave single-station direction finding and elevation measurement combined positioning method corresponding to the above method embodiments, in order to achieve the corresponding functions and technical effects. See also Figure 2 , Figure 2 This diagram illustrates a structural block diagram of a shortwave monostation direction-finding and elevation-measuring combined positioning system 200 provided in an embodiment of this application. For ease of explanation, only the parts relevant to this embodiment are shown. The shortwave monostation direction-finding and elevation-measuring combined positioning system 200 provided in this embodiment includes: Angle extraction module 201 is used to acquire shortwave antenna array received data, and simultaneously extract weighted azimuth sequence and elevation confidence interval based on the shortwave antenna array received data, and compress azimuth confidence range based on the weighted azimuth sequence and elevation confidence interval to construct wave-arrival angle parameter set; The virtual height correction module 202 is used to extract the elevation angle estimate based on the arrival angle parameter set, separate the multipath delay difference component from the data received by the shortwave antenna array to reverse calculate the real-time ionospheric height deviation, perform self-correction of the ionospheric virtual height based on the real-time ionospheric height deviation to obtain the corrected virtual height value, and combine the corrected virtual height value with the elevation angle estimate and the arrival angle parameter set to jointly perform reflection slant range calculation to determine the candidate positioning coordinate set. The multi-hop verification module 203 is used to identify multi-hop propagation components in the shortwave antenna array received data, extract ground reflection constraints for each hop, perform geometric consistency verification between the ground reflection constraints for each hop and the candidate positioning coordinate set to form a coordinate confidence level table, and establish a slant range optimization configuration by weighting the candidate positioning coordinate set according to the coordinate confidence level table. The fusion positioning module 204 is used to construct a multi-dimensional positioning parameter set by collaborative weighted clustering of the candidate positioning coordinate set based on the coordinate confidence level table and the slant distance preferred configuration, and to perform radiation source location fusion calculation based on the multi-dimensional positioning parameter set to determine the final positioning coordinates.

[0065] The aforementioned shortwave single-station direction-finding and elevation-measuring joint positioning system 200 can implement one of the shortwave single-station direction-finding and elevation-measuring joint positioning methods described in the above-described method embodiments. The options in the above method embodiments are also applicable to this embodiment and will not be detailed here. The remaining content of this application's embodiments can be referred to the content of the above method embodiments, and will not be repeated in this embodiment.

[0066] The above embodiments are not an exhaustive list based on the present invention, and there may be many other embodiments not listed. Any substitutions and improvements made without departing from the concept of the present invention are within the protection scope of the present invention.

Claims

1. A shortwave single-station direction finding and elevation measurement combined positioning method, characterized in that, include: Acquire shortwave antenna array received data, and simultaneously extract weighted azimuth sequence and elevation confidence interval based on the shortwave antenna array received data. Compress the azimuth confidence range based on the weighted azimuth sequence and the elevation confidence interval to construct the arrival angle parameter set. Based on the arrival angle parameter set, the elevation angle estimate is extracted. The multipath delay difference component is separated from the data received by the shortwave antenna array to reversely calculate the real-time ionospheric height deviation. Based on the real-time ionospheric height deviation, the ionospheric virtual height is self-corrected to obtain the corrected virtual height value. The corrected virtual height value and the elevation angle estimate are combined with the arrival angle parameter set to jointly calculate the reflection slant range and determine the candidate positioning coordinate set. In the shortwave antenna array received data, multi-hop propagation components are identified and ground reflection constraints of each hop are extracted. Geometric consistency verification is performed between each hop ground reflection constraint and the candidate positioning coordinate set to form a coordinate confidence level table. The candidate positioning coordinate set is then weighted and optimized according to the coordinate confidence level table to establish a slant range optimization configuration. Based on the coordinate confidence level table and the slant distance optimization configuration, the candidate positioning coordinate set is collaboratively weighted and clustered to construct a multidimensional positioning parameter set. Based on the multidimensional positioning parameter set, the radiation source location fusion calculation is performed to determine the final positioning coordinates.

2. The method according to claim 1, characterized in that, The step of constructing the arrival angle parameter set by compressing the azimuth confidence range based on the weighted azimuth sequence and the elevation confidence interval includes: A low elevation angle perturbation sensitivity analysis is performed on the elevation angle confidence interval to generate an elevation angle error amplification factor sequence; Based on the elevation error amplification factor sequence, the weighted azimuth sequence is subjected to error coupling quantization to obtain the upper bound of the azimuth coupling error; The weighted azimuth sequence is subjected to interval shrinkage based on the upper bound of the azimuth coupling error to generate a compressed azimuth distribution. The compressed azimuth distribution and the elevation confidence interval are jointly encapsulated to construct the incoming wave angle parameter set.

3. The method according to claim 1, characterized in that, The step of separating multipath delay difference components from the data received by the shortwave antenna array and then calculating the real-time ionospheric height deviation includes: The received data from the shortwave antenna array is subjected to time-delay domain separation to extract the main path component and the secondary path component; The path delay difference sequence is generated by differentiating the primary path component and the secondary path component. The time-varying slope is extracted from the path delay difference sequence to generate an estimate of the ionospheric drift rate; The real-time ionospheric height deviation is generated by dynamically weighting the path delay difference sequence based on the estimated ionospheric drift rate.

4. The method according to claim 1, characterized in that, The step of combining the corrected false height value with the estimated elevation angle value and the wave arrival angle parameter set to jointly calculate the reflection slant range and determine the candidate positioning coordinate set includes: Based on the set of incoming wave angle parameters, an azimuth probability density sequence is extracted; Based on the azimuth probability density sequence, the corrected false altitude value and the elevation angle estimate are used to jointly drive the generation of a multi-hypothesis slant range set; The weighted slope range distribution is obtained by performing probability weighting filtering on the azimuth probability density sequence of the multi-hypothesis slope range set. The candidate positioning coordinate set is determined by expanding the polar coordinate projection according to the weighted slant distance distribution.

5. The method according to claim 1, characterized in that, The step of identifying multi-hop propagation components and extracting ground reflection constraints for each hop in the received data from the shortwave antenna array includes: The shortwave antenna array receives data and performs hop count component separation to generate waveform sets for each hop. For each hop arrival waveform set, extract the inter-hop delay difference to generate an inter-hop delay difference sequence; The horizontal gradient correction amount is generated by identifying the proportional deviation of the jump delay difference sequence. Gradient compensation is applied to the position of each hop reflection point based on the horizontal gradient correction amount to generate ground reflection constraints for each hop.

6. The method according to claim 1, characterized in that, The step of performing geometric consistency verification between the ground reflection constraints of each jump and the candidate positioning coordinate set to form a coordinate reliability grading table includes: A multi-hop geometric sequence is established by arranging the ground reflection constraints of each hop in order of hop number; Based on the multi-hop geometric sequence, the radiation source position constraints of each hop are calculated by extending backwards towards the radiation source direction to generate a radiation source position constraint set. The spatial overlap quantification of the radiation source location constraint set and the candidate positioning coordinate set generates an overlap sequence. The coordinate reliability grading table is formed by classifying the overlap sequence.

7. The method according to claim 1, characterized in that, The step of constructing a multidimensional positioning parameter set by collaborative weighted clustering of the candidate positioning coordinate set based on the coordinate confidence level table and the slant range optimization configuration includes: Based on the coordinate confidence level table and the preferred slant distance configuration, the candidate positioning coordinate set is weighted, filtered, and clustered to generate a clustering residual sequence; The coordinate confidence level table is reverse-corrected based on the clustering residual sequence to generate an updated confidence level table; The candidate location coordinate set is re-weighted and clustered according to the updated confidence level table to generate convergent location coordinates; The converged positioning coordinates, the updated confidence level table, and the slant distance optimization configuration are jointly encapsulated to construct a multi-dimensional positioning parameter set.

8. The method according to claim 5, characterized in that, The step of generating a horizontal gradient correction amount by identifying the proportional deviation of the jump delay difference sequence includes: Construct a hop-by-hop delay difference ratio sequence for the hop-to-hop delay difference sequence based on adjacent hop pairs; Based on the hop-by-hop delay difference ratio sequence, the median reference value is extracted to generate an equal-ratio reference sequence; The trend slope and trend sign are extracted by fitting a linear trend between the hop-by-hop delay difference ratio sequence and the equal ratio benchmark sequence along the hop order direction. The horizontal gradient correction amount is generated by mapping the ionospheric horizontal gradient magnitude and gradient direction based on the trend slope and trend sign.

9. The method according to claim 6, characterized in that, The step of generating a radiation source position constraint set by extending the multi-hop geometric sequence backwards towards the radiation source direction to calculate the position constraints of each hop of the radiation source includes: The incident elevation angle estimates are extracted from each hop reflection point in the multi-hop geometric sequence to generate a set of hop incident elevation angles; Based on the set of incident elevation angles for each jump, the upper and lower bounds of the radiation source distances are calculated by extending each jump reflection point in the multi-hop geometric sequence along the reverse propagation path to generate a set of radiation source distance intervals for each jump. The distance intervals of each hop radiation source are weighted and overlapped according to the number of hops to generate a fused radiation source distance interval; Based on the fused radiation source distance interval and the azimuth information of the multi-hop geometric sequence, a spatial projection is performed to generate a radiation source position constraint set.

10. A shortwave single-station direction finding and elevation measurement combined positioning system, characterized in that, include: An angle extraction module is used to acquire shortwave antenna array received data, and simultaneously extract weighted azimuth sequence and elevation confidence interval based on the shortwave antenna array received data. Based on the weighted azimuth sequence and the elevation confidence interval, the azimuth confidence range is compressed to construct the incoming wave angle parameter set. The virtual height correction module is used to extract the elevation angle estimate based on the arrival angle parameter set, separate the multipath delay difference component from the data received by the shortwave antenna array to back-calculate the real-time ionospheric height deviation, self-correct the ionospheric virtual height based on the real-time ionospheric height deviation to obtain the corrected virtual height value, and combine the corrected virtual height value with the elevation angle estimate and the arrival angle parameter set to jointly perform reflection slant range calculation to determine the candidate positioning coordinate set. The multi-hop verification module is used to identify multi-hop propagation components in the shortwave antenna array received data, extract ground reflection constraints for each hop, perform geometric consistency verification between the ground reflection constraints for each hop and the candidate positioning coordinate set to form a coordinate confidence level table, and establish a slant range optimization configuration by weighting the candidate positioning coordinate set according to the coordinate confidence level table. The fusion positioning module is used to construct a multi-dimensional positioning parameter set by collaborative weighted clustering of the candidate positioning coordinate set based on the coordinate confidence level table and the slant range preferred configuration, and to perform radiation source location fusion calculation based on the multi-dimensional positioning parameter set to determine the final positioning coordinates.

Images