Method for calculating azimuth angle of weather modification operation in real time based on Beidou satellite

Abstract

The invention provides a weather modification operation azimuth angle real-time calculation method based on Beidou satellites, which comprises the following steps: firstly, acquiring a Beidou satellite positioning data set of a target area, including a three-dimensional coordinate sequence of a plurality of positioning nodes in a continuous time period, and then, carrying out dynamic track fitting processing on the Beidou satellite positioning data set, performing environment interference feature analysis based on the dynamic positioning feature set, generating an environment interference feature set, determining a real-time azimuth angle compensation parameter set according to a coupling association relationship between the dynamic positioning feature and the environment interference feature, and finally, performing real-time azimuth angle compensation on the real-time azimuth angle compensation parameter set. And azimuth angle iteration resolving processing is carried out on the positioning data set based on the real-time azimuth angle compensation parameter set to generate a weather modification operation real-time azimuth angle parameter of the target area, so that real-time and accurate calculation of the weather modification operation azimuth angle is realized, and the operation effect and safety are improved.

Description

Technical Field

[0001] The present invention relates to the fields of meteorological technology and satellite navigation and positioning technology, and in particular to a real-time calculation method for the azimuth angle of a human shadow operation based on Beidou satellites. Background Art

[0002] In the field of meteorology, weather modification (WMO) operations, such as cloud seeding and hail suppression, are of great significance to agricultural production, water resource management, and disaster prevention. The precise implementation of WMO operations relies on accurate calculation of the operation azimuth to maximize the effectiveness of the operation. Traditional methods for calculating WMO azimuths rely primarily on ground-based observation equipment and empirical formulas. These methods can provide relatively accurate results in static or low-dynamic environments. However, in complex and changing natural environments, such as mountainous areas, oceans, or areas with drastic changes in meteorological conditions, ground-based observation equipment may be affected by factors such as terrain obstruction and meteorological interference, resulting in significant errors in azimuth calculations. Furthermore, with technological advancements, the requirements for real-time and accurate WMO azimuths are becoming increasingly stringent, and traditional methods are no longer able to meet these requirements. Therefore, developing a method that can accurately calculate WMO azimuths in real time has become an urgent task. Summary of the Invention

[0003] In view of the above-mentioned problems, in combination with the first aspect of the present invention, an embodiment of the present invention provides a method for real-time calculation of the azimuth angle of human shadow operation based on Beidou satellites, the method comprising: Acquire a Beidou satellite positioning data set for a target area, the Beidou satellite positioning data set comprising a three-dimensional coordinate sequence of multiple positioning nodes within a continuous time period; Performing dynamic trajectory fitting processing on the Beidou satellite positioning data set to obtain a dynamic positioning feature set including an instantaneous displacement deviation and a continuous trajectory curvature change; Performing environmental interference feature analysis processing on the target area based on the dynamic positioning feature set to generate an environmental interference feature set; Determining a real-time azimuth compensation parameter set according to a coupling association relationship between the dynamic positioning feature set and the environmental interference feature set; An azimuth iterative solution process is performed on the Beidou satellite positioning data set based on the real-time azimuth compensation parameter set to generate a real-time azimuth parameter of the human shadow operation in the target area.

[0004] In a possible implementation of the first aspect, performing dynamic trajectory fitting processing on the Beidou satellite positioning data set to obtain a dynamic positioning feature set including an instantaneous displacement deviation and a continuous trajectory curvature change includes: Extracting spatiotemporal distribution information of positioning nodes in the three-dimensional coordinate sequence to generate spatiotemporal correlation features; Adaptively segmenting the spatiotemporal correlation features using a sliding window to obtain multiple windowed trajectory segments; dynamically adjusting the time span of the sliding window based on the spatial distribution density of the positioning nodes; the lower the spatial density, the smaller the time span of the sliding window; and employing a Kalman filter to perform noise smoothing on discrete coordinate points within the sliding window. Adopting an adaptive piecewise interpolation method to perform nonlinear reconstruction processing on the discrete coordinate points in the windowed trajectory segments to generate a smooth trajectory curve; Calculating the displacement deviation between adjacent positioning nodes based on the smooth trajectory curve to generate the instantaneous displacement deviation; Calculating the curvature derivatives of the smooth trajectory curve at adjacent interpolation points, determining the direction of curvature change according to the change in the sign of the curvature derivative, and extracting the curvature change of the continuous trajectory; The instantaneous displacement deviation and the continuous trajectory curvature change are combined into a dynamic positioning feature set.

[0005] In a possible implementation of the first aspect, performing environmental interference feature analysis on the target area based on the dynamic positioning feature set to generate the environmental interference feature set includes: Acquire a meteorological monitoring data set of a target area, wherein the meteorological monitoring data set includes wind speed gradient distribution data and atmospheric refractive index distribution data; The wind speed gradient distribution data is vector-decomposed to obtain transverse and longitudinal wind speed components. After the instantaneous displacement deviation is standardized, a nonlinear regression model of the instantaneous displacement deviation and the wind speed component is established. A feature interaction matrix is ​​constructed by introducing quadratic terms and cross terms of the wind speed components. The regression model is trained using a random forest algorithm and feature importance scores are calculated. The dynamic contribution weight of each wind speed component is determined in combination with standardized residual analysis. The nonlinear residual part is compensated and corrected using Gaussian process regression to generate a wind speed interference feature. Performing correlation analysis on the atmospheric refractive index distribution data and the continuous trajectory curvature variation in the dynamic positioning feature set, calculating the covariance matrix of the atmospheric refractive index spatial gradient and the trajectory curvature variation, and extracting the principal component feature as the refractive index interference factor to generate a refractive index interference feature; Constructing an environmental interference feature set according to the wind speed interference feature and the refractive index interference feature; Among them, the environmental interference feature set includes a wind speed contribution weight set and a refractive index interference factor set, and each positioning node corresponds to a set of wind speed contribution weights and refractive index interference factors; the wind speed contribution weight set is used to quantify the degree of influence on the positioning node displacement deviation under different wind speed gradients; the refractive index interference factor set is used to characterize the intensity of the nonlinear effect of the spatial distribution of atmospheric refractive index on the change of trajectory curvature.

[0006] In a possible implementation of the first aspect, determining the real-time azimuth compensation parameter set according to the coupling association relationship between the dynamic positioning feature set and the environmental interference feature set includes: Performing weight distribution processing on the instantaneous displacement deviation in the dynamic positioning feature set and the wind speed contribution weight in the environmental interference feature set to generate a first compensation weight coefficient; Performing weight distribution processing on the continuous trajectory curvature change in the dynamic positioning feature set and the refractive index interference factor in the environmental interference feature set to generate a second compensation weight coefficient; Constructing a dynamic compensation weight set according to the first compensation weight coefficient and the second compensation weight coefficient; the dynamic compensation weight set includes a displacement compensation weight and a curvature compensation weight corresponding to each positioning node; After normalizing the instantaneous displacement deviation and the continuous trajectory curvature change respectively, the fusion weight coefficient of each normalized result is calculated by the entropy weight method based on the displacement compensation weight and the curvature compensation weight in the dynamic compensation weight set, and the weighted fusion processing of each normalized result is performed using the weighted geometric mean algorithm to generate a real-time azimuth compensation parameter set; the real-time azimuth compensation parameter set includes a displacement compensation parameter subset and a curvature compensation parameter subset, each subset is arranged in time series and corresponds one-to-one to a positioning node.

[0007] In a possible implementation of the first aspect, performing azimuth iterative solution processing on the Beidou satellite positioning data set based on the real-time azimuth compensation parameter set to generate real-time azimuth parameters of human shadow operations in the target area includes: Calling the three-dimensional coordinate sequence in the Beidou satellite positioning data set to construct an initial azimuth angle solution model; the initial azimuth angle solution model uses the least squares method to calculate the azimuth angle reference value; Inputting the real-time azimuth compensation parameter set into the initial azimuth solution model to generate a compensated azimuth solution model; the compensated azimuth solution model includes a displacement compensation module and a curvature correction module, which are respectively used to process the displacement compensation parameter subset and the curvature compensation parameter subset; The compensated azimuth angle solution model is iteratively optimized using a gradient descent algorithm until the compensated azimuth angle solution model converges, and then azimuth angle parameters are extracted to generate real-time azimuth angle parameters for human shadow operation.

[0008] In a possible implementation of the first aspect, iteratively optimizing the compensated azimuth angle calculation model using a gradient descent algorithm until the compensated azimuth angle calculation model converges, extracting azimuth angle parameters, and generating real-time azimuth angle parameters for shadow operation include: Get the azimuth angle solution error and error gradient direction of the current iteration cycle; Calculating the product of the gradient direction and the learning rate according to the error gradient direction, and incrementally updating the compensation weight coefficient to adjust the compensation weight coefficient in the real-time azimuth compensation parameter set; Recalculating the outputs of the displacement compensation module and the curvature correction module based on the adjusted compensation weight coefficient, and synchronously updating the internal state variables of the model to update the compensated azimuth angle solution model; Repeat the above operations until the compensated azimuth angle solution model meets the convergence conditions, extract the azimuth angle parameters in the converged azimuth angle solution model, and generate the real-time azimuth angle parameters of the human shadow operation. The convergence conditions include that the error change rate of the azimuth angle solution error is less than the set threshold or reaches the maximum number of iterations.

[0009] In a possible implementation of the first aspect, obtaining the azimuth angle solution error and the error gradient direction of the current iteration cycle includes: Extracting the measured azimuth angle data of the human shadow operation in the current iteration cycle, wherein the measured azimuth angle data of the human shadow operation is collected by the angle measurement equipment of the ground reference station and is synchronized with the Beidou positioning data; Calling the compensated azimuth angle solution model to generate azimuth angle prediction data for a current iteration cycle, the azimuth angle prediction data including an azimuth angle sequence output by the compensated azimuth angle solution model and its confidence interval; Performing time alignment processing on the azimuth measured data and the azimuth predicted data, and calculating the absolute error and relative error point by point to obtain the azimuth solution error; Based on the azimuth angle solution error, the compensation weight coefficients in the real-time azimuth angle compensation parameter set are back-propagated to generate the error gradient direction; wherein, the back-propagation processing includes constructing a calculation graph and deriving the topology structure of the compensated azimuth angle solution model layer by layer to obtain the contribution gradient of each compensation weight coefficient to the total error.

[0010] In a possible implementation of the first aspect, calculating the product of the gradient direction and the learning rate according to the error gradient direction, and incrementally updating the compensation weight coefficient to adjust the compensation weight coefficient in the real-time azimuth compensation parameter set includes: Normalizing the error gradient direction to generate a gradient adjustment step size, wherein the normalization includes calculating the L2 norm of the gradient vector and dividing each gradient component by the norm value; updating the first compensation weight coefficient and the second compensation weight coefficient according to the gradient adjustment step size; The updated first compensation weight coefficient and the second compensation weight coefficient are rewritten into the dynamic compensation weight set. A spatial regularization constraint is introduced during the update so that the difference in the compensation weight coefficients of adjacent positioning nodes does not exceed a set threshold, and the weight records of the most recent N iterations are retained; wherein the value of N is inversely proportional to the spatial distribution density of the positioning nodes.

[0011] In a possible implementation of the first aspect, extracting the azimuth parameters from the converged azimuth solution model to generate real-time azimuth parameters for the human shadow operation includes: Obtaining a parameter space distribution of the converged azimuth solution model, wherein the parameter space distribution includes a weight matrix of a displacement compensation module and a bias vector of a curvature correction module; Calculating the eigenvectors of the parameter covariance matrix of the parameter space distribution, and selecting the first K principal components in order of the eigenvector values ​​to extract the azimuth principal component features; A kernel density estimation method is used to perform non-parametric probability modeling on the principal component characteristics to construct a probability density function of the azimuth parameter; The final real-time azimuth angle parameter of the human shadow operation is determined according to the extreme value point of the probability density function, and the extreme value point is obtained by traversing the definition domain of the probability density function and finding the local maximum point.

[0012] In a possible implementation of the first aspect, performing non-parametric probability modeling on the principal component features using a kernel density estimation method to construct a probability density function of the azimuth parameter includes: Performing kernel density estimation on the azimuth principal component feature to generate an initial probability density distribution; the kernel density estimation process uses a Gaussian kernel function, and the bandwidth parameter is adaptively adjusted according to the variance of the principal component feature; Calibrate the initial probability density distribution according to the three-dimensional coordinate sequence in the Beidou satellite positioning data set to generate a calibrated probability density function; The peak area of ​​the calibrated probability density function is extracted, and the median point of the peak area is used as the real-time azimuth parameter of the human shadow operation; wherein, the peak area is obtained by binarizing the probability density function by setting a density threshold, and the median point is calculated using the regional centroid algorithm.

[0013] On the other hand, an embodiment of the present invention also provides a real-time calculation system for the azimuth angle of human shadow operations based on Beidou satellites, including a processor and a machine-readable storage medium, the machine-readable storage medium being connected to the processor, the machine-readable storage medium being used to store programs, instructions or codes, and the processor being used to execute the programs, instructions or codes in the machine-readable storage medium to implement the above method.

[0014] Based on the above aspects, the embodiment of the present invention obtains Beidou satellite positioning data of the target area and performs dynamic trajectory fitting processing on it, so as to accurately extract dynamic positioning features including instantaneous displacement deviation and continuous trajectory curvature change. Furthermore, by performing environmental interference feature analysis on these dynamic positioning features, it is possible to accurately identify environmental interference factors that affect azimuth calculation. Based on the coupling correlation between dynamic positioning features and environmental interference features, real-time azimuth compensation parameters are determined, which effectively eliminates the influence of environmental interference on azimuth calculation. Finally, through azimuth iterative solution processing, real-time azimuth parameters of human shadow operations in the target area are generated, which significantly improves the real-time and accuracy of azimuth calculation, thereby significantly improving the implementation effect and safety of human shadow operations, and providing a data basis for artificial weather modification operations in the field of meteorology. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 The present invention provides a schematic diagram of the execution flow of a method for real-time calculation of the azimuth angle of a shadow operation based on Beidou satellites.

[0016] Figure 2 It is a schematic diagram of exemplary hardware and software components of a Beidou satellite-based real-time calculation system for human shadow operation azimuth angles provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0017] The present invention will be described in detail below with reference to the accompanying drawings. Figure 1 The present invention provides a flowchart of a method for real-time calculation of the azimuth angle of a human shadow operation based on a Beidou satellite according to an embodiment of the present invention. The method for real-time calculation of the azimuth angle of a human shadow operation based on a Beidou satellite is introduced in detail below.

[0018] Step S110: Acquire a Beidou satellite positioning data set of a target area, wherein the Beidou satellite positioning data set includes a three-dimensional coordinate sequence of multiple positioning nodes within a continuous time period.

[0019] In the scenario of weather modification operations in the field of meteorology, in order to accurately calculate the azimuth angle of weather modification operations, it is first necessary to obtain the Beidou satellite positioning data set of the target area. The target area can be a specific airspace where operations such as artificial rainfall enhancement and hail prevention are planned. Within this target area, multiple positioning nodes can be pre-arranged. These positioning nodes can be Beidou positioning modules installed on meteorological monitoring equipment, unmanned aerial vehicles, or other tools used for operations.

[0020] These positioning nodes continuously record their respective position information within a continuous time period. Assume the time period is represented by T, from time t1 to time t2 (t1 < t2). Each positioning node generates a three-dimensional coordinate (x, y, z) at each moment, where x, y, and z respectively represent the position information on different coordinate axes. The three-dimensional coordinate data recorded by multiple positioning nodes within the time period T constitutes a three-dimensional coordinate sequence. For example, assume there are m positioning nodes, and i represents the number of the positioning node (i = 1, 2,..., m). For each positioning node i, there will be a series of three-dimensional coordinates (xi1, yi1, zi1), (xi2, yi2, zi2),..., (xin, yin, zin) within the time period T, where n represents the number of coordinates recorded by this positioning node within the time period T. The combination of the three-dimensional coordinate sequences of all these positioning nodes forms the Beidou satellite positioning data set.

[0021] Step S120: Perform dynamic trajectory fitting processing on the Beidou satellite positioning data set to obtain a dynamic positioning feature set including instantaneous displacement deviation and continuous trajectory curvature change.

[0022] After obtaining the Beidou satellite positioning data set, due to the fact that the positioning nodes are affected by various factors in the actual environment and their positions change, in order to accurately analyze their movement trajectories and characteristics, it is necessary to perform dynamic trajectory fitting processing on this data set.

[0023] Step S121: Extract the spatio-temporal distribution information of the positioning nodes in the three-dimensional coordinate sequence to generate spatio-temporal correlation features.

[0024] For the three-dimensional coordinate sequence, each coordinate point contains time and space information. The time information is reflected in which moment this coordinate point is recorded, and the space information is the three-dimensional coordinate (x, y, z). First, arrange the three-dimensional coordinate sequences of each positioning node in chronological order. Assume the three-dimensional coordinate sequence of positioning node i is (xi1, yi1, zi1, t1), (xi2, yi2, zi2, t2),..., (xin, yin, zin, tn), where ti represents the moment when this coordinate point is recorded.

[0025] Next, analyze the spatial positional relationships of the positioning nodes at different times. This relationship can be demonstrated by calculating the distance and time interval between coordinate points at adjacent times. For example, for positioning node i, calculate the spatial distance dik = √[(xik+1-xik)² + (yik+1-yik)² + (zik+1-zik)²] between the coordinate points (xik, yik, zik) corresponding to adjacent times tk and tk+1 (the spatial distance dik is the square root of the sum of the squared differences in the x, y, and z directions of the corresponding coordinate points at adjacent times), as well as the time interval Δtk = tk+1-tk.

[0026] By combining this distance and time interval information for all positioning nodes, we can obtain the spatiotemporal distribution information of the positioning nodes. Furthermore, to generate spatiotemporal correlation features, this spatiotemporal distribution information can be encoded. For example, a vector can be formed by sequentially arranging the series of distance and time interval values ​​for each positioning node. Assuming that positioning node i has n coordinate points, this will form a vector of length 2(n-1), where the first n-1 elements are the spatial distances between adjacent coordinate points, and the last n-1 elements are the time intervals between adjacent coordinate points. Combining these vectors for all m positioning nodes constitutes the spatiotemporal correlation feature.

[0027] Step S122: Adaptively perform sliding window segmentation processing on the spatiotemporal correlation features to obtain multiple windowed trajectory segments; wherein, the time span of the sliding window is dynamically adjusted according to the spatial distribution density of the positioning nodes. The lower the spatial density, the smaller the sliding window time span. Kalman filtering is used to smooth the noise of discrete coordinate points in the sliding window.

[0028] After obtaining the spatiotemporal correlation features, in order to further analyze the trajectory characteristics of the positioning nodes, it is necessary to perform adaptive sliding window segmentation processing on them. The sliding window is a fixed-size time period that contains the coordinate information of multiple positioning nodes within this time period.

[0029] First, the time span of the sliding window must be dynamically adjusted based on the spatial distribution density of the positioning nodes. The spatial distribution density of positioning nodes can be measured by counting the number of positioning nodes per unit volume. Assuming there are m positioning nodes in a specific spatial region V, the spatial distribution density ρ = m / V.

[0030] If the spatial distribution density is low, it means that the distance between the positioning nodes is large and their movement changes may be more drastic. In this case, a smaller sliding window time span is required to more accurately capture their trajectory changes. For example, when the spatial distribution density is lower than a certain threshold ρ0, the sliding window time span is set to Δt1; when the spatial distribution density is higher than ρ0, the sliding window time span is set to Δt2, and Δt1<Δt2.

[0031] Then, a sliding window is used to slide across the spatiotemporal correlation features. Assuming the sliding window starts at tstart and ends at tend = tstart + Δt (Δt is the time span of the sliding window), the coordinate information of multiple positioning nodes within this time period is included. This coordinate information constitutes a windowed trajectory segment. As the sliding window continues to slide, multiple windowed trajectory segments are obtained.

[0032] The discrete coordinate points within each windowed trajectory segment may contain noise due to factors such as external environmental interference. Kalman filtering is used to eliminate this noise. Kalman filtering is a recursive optimal estimation algorithm that uses the system's state equation and observation equation, combining the estimated value at the previous moment with the observed value at the current moment, to estimate the current system state. In this scenario, the system state can be considered the true position of the positioning node, and the observed value is the discrete coordinate point within the windowed trajectory segment. Kalman filtering can smooth these discrete coordinate points to obtain more accurate trajectory information.

[0033] Step S123: Adopting an adaptive piecewise interpolation method to perform nonlinear reconstruction processing on the discrete coordinate points in the windowed trajectory segment to generate a smooth trajectory curve.

[0034] After obtaining the windowed trajectory segments and performing noise smoothing, since there may be discontinuities between discrete coordinate points, in order to obtain a more continuous and smooth trajectory, an adaptive piecewise interpolation method is used to perform nonlinear reconstruction on these discrete coordinate points.

[0035] The adaptive piecewise interpolation method divides the windowed trajectory into multiple segments based on the distribution of discrete coordinate points. For each discrete coordinate point within a segment, an appropriate interpolation function is used. For example, a cubic spline interpolation function can be used, which ensures that the function value, first-order derivative, and second-order derivative at the interpolation point are continuous, resulting in a smoother curve.

[0036] When segmenting, adaptive adjustments should be made based on the curvature changes of the discrete coordinate points. If the curvature of a segment of discrete coordinate points changes significantly, it means that the trajectory changes more drastically. In this case, the segment should be divided into smaller segments to improve the accuracy of interpolation. If the curvature change is small, the segment can be divided into larger segments.

[0037] For each discrete coordinate point within each segment (x1, y1, z1), (x2, y2, z2), …, (xk, yk, zk), a cubic spline interpolation function is used to interpolate the x, y, and z coordinates. For the x coordinate, for example, a series of linear equations are solved to determine the coefficients of the cubic spline interpolation function, thus generating the interpolation curve for the x coordinate within that segment. Similarly, the y and z coordinates are processed similarly. By combining the interpolation curves for all segments, a smooth trajectory curve is generated.

[0038] Step S124: Calculating the displacement deviation between adjacent positioning nodes based on the smooth trajectory curve to generate the instantaneous displacement deviation.

[0039] After obtaining the smooth trajectory curve, in order to analyze the relative motion between positioning nodes, it is necessary to calculate the displacement deviation between adjacent positioning nodes to generate the instantaneous displacement deviation.

[0040] For each time instant t on the smooth trajectory curve, each positioning node has a corresponding position. Assume that at time t, the position of positioning node i is (xi(t), yi(t), zi(t)), and the position of positioning node i+1 is (xi+1(t), yi+1(t), zi+1(t)). Then, the displacement deviation di,i+1(t) between adjacent positioning nodes i and i+1 at time t = √[(xi+1(t)-xi(t))²+(yi+1(t)-yi(t))²+(zi+1(t)-zi(t))²] (the displacement deviation is the square root of the sum of the squares of the position differences of adjacent positioning nodes in the x, y, and z directions).

[0041] Combining the displacement deviations at all moments yields the instantaneous displacement deviations between adjacent positioning nodes. This calculation is repeated for all positioning node pairs (i, i+1) (i=1, 2, …, m-1), ultimately yielding the instantaneous displacement deviation for the entire set of positioning nodes.

[0042] Step S125: Calculate the curvature derivatives of the smooth trajectory curve at adjacent interpolation points, determine the direction of curvature change according to the change in the sign of the curvature derivative, and extract the curvature change of the continuous trajectory.

[0043] In order to further analyze the curvature change of the positioning node trajectory, it is necessary to calculate the curvature derivative of the smooth trajectory curve at adjacent interpolation points, determine the direction of curvature change according to its sign change, and then extract the curvature change of the continuous trajectory.

[0044] First, the curvature of the smooth trajectory curve at each interpolation point is calculated. Curvature describes the curvature of a curve and can be calculated using the first and second derivatives of the curve. Calculating the curvature κ(t) of the smooth trajectory curve at the interpolation point (x(t), y(t), z(t)) involves taking the first and second derivatives of x(t), y(t), and z(t), respectively, and then calculating it using conventional technical formulas (the specific formula can be derived from the definition of spatial curve curvature).

[0045] Next, we calculate the curvature derivatives at adjacent interpolation points. Assuming that the curvatures corresponding to adjacent interpolation points t1 and t2 are κ(t1) and κ(t2), respectively, the curvature derivative Δκ = (κ(t2) - κ(t1)) / (t2 - t1).

[0046] The direction of curvature change is determined by the change in the sign of the curvature derivative. If the curvature derivative is greater than 0, it means that the curvature is increasing, that is, the curvature of the curve is increasing; if the curvature derivative is less than 0, it means that the curvature is decreasing, that is, the curvature of the curve is decreasing.

[0047] By combining the curvature derivatives at all adjacent interpolation points and the corresponding curvature change direction information, the curvature change of the continuous trajectory is extracted.

[0048] Step S126: combining the instantaneous displacement deviation and the continuous trajectory curvature change into a dynamic positioning feature set.

[0049] After obtaining the instantaneous displacement deviation and the continuous trajectory curvature change, they are combined into a dynamic positioning feature set. The instantaneous displacement deviation and the continuous trajectory curvature change can be regarded as two vectors, which are then concatenated. For example, assuming the instantaneous displacement deviation vector is D = [d1, d2, …, dn] and the continuous trajectory curvature change vector is C = [c1, c2, …, cn] (where n represents the length of the vector, corresponding to the number of positioning nodes or the number of time points), then the dynamic positioning feature set F can be expressed as F = [D, C]. That is, the two vectors are concatenated in sequence to form a new vector. This new vector contains information about the instantaneous displacement deviation and continuous trajectory curvature change of the positioning node, providing a basis for subsequent environmental interference feature analysis and processing.

[0050] Step S130: performing environmental interference feature analysis processing on the target area based on the dynamic positioning feature set to generate an environmental interference feature set.

[0051] After obtaining the dynamic positioning feature set, since the environmental factors in the target area will interfere with the movement of the positioning node, thereby affecting the calculation of the azimuth angle of the human shadow operation, it is necessary to perform environmental interference feature analysis on the target area to generate an environmental interference feature set.

[0052] Step S131: Acquire a meteorological monitoring data set of a target area, wherein the meteorological monitoring data set includes wind speed gradient distribution data and atmospheric refractive index distribution data.

[0053] To analyze environmental interference characteristics, we first need to obtain a set of meteorological monitoring data for the target area. Multiple meteorological monitoring devices are deployed within the target area, providing real-time meteorological data. This data primarily includes wind speed gradient distribution data and atmospheric refractive index distribution data.

[0054] Wind speed gradient distribution data describes the variation in wind speed at different locations within a target area. Meteorological monitoring equipment measures wind speed at different altitudes and horizontal positions, generating a series of wind speed data. By organizing this wind speed data according to location information, wind speed gradient distribution data can be generated. For example, wind speed values ​​vij (where i represents the altitude and j represents the horizontal position) are measured at different altitudes h1, h2, …, hk, and at different horizontal positions (x1, y1), (x2, y2), …, (xl, yl)). These wind speed values ​​are combined with the corresponding location information to form the wind speed gradient distribution data.

[0055] Atmospheric refractive index distribution data describes the spatial distribution of atmospheric refractive index within a target area. Atmospheric refractive index affects the propagation of electromagnetic waves, thereby interfering with Beidou satellite signal transmission. Meteorological monitoring equipment measures atmospheric parameters such as temperature, humidity, and pressure, and calculates the atmospheric refractive index using conventional formulas. Similarly, atmospheric refractive index data at different locations is organized according to location information to generate atmospheric refractive index distribution data.

[0056] Step S132: Perform vector decomposition on the wind speed gradient distribution data to obtain transverse wind speed components and longitudinal wind speed components, and after standardizing the instantaneous displacement deviation, establish a nonlinear regression model of the instantaneous displacement deviation and the wind speed component, construct a feature interaction matrix by introducing the quadratic terms and cross terms of the wind speed components, use the random forest algorithm to train the regression model and calculate the feature importance score, combine the standardized residual analysis to determine the dynamic contribution weight of each wind speed component, use Gaussian process regression to compensate and correct the nonlinear residual part, and generate the wind speed interference feature.

[0057] After obtaining the wind speed gradient distribution data, it is first vectorized. Wind speed is a vector with magnitude and direction. Assuming the wind speed vector is V, its projection on the horizontal plane can be decomposed into a transverse wind speed component Vx and a longitudinal wind speed component Vy. Based on the wind speed direction angle θ, the transverse wind speed component Vx = |V| * cos(θ) and the longitudinal wind speed component Vy = |V| * sin(θ) can be calculated using trigonometric functions, where |V| represents the magnitude of the wind speed.

[0058] To eliminate the influence of varying magnitudes, the instantaneous displacement deviations in a dynamic positioning feature set require normalization. Common methods, such as Z-score normalization, can be used for this purpose. Assuming the instantaneous displacement deviation is D, its mean is μD, and its standard deviation is σD, then the normalized instantaneous displacement deviation D' = (D - μD) / σD.

[0059] Next, a nonlinear regression model is established between the instantaneous displacement deviation and wind speed components. To capture the complex relationship between them, quadratic and cross terms of the wind speed components are introduced to construct a feature interaction matrix. For example, the transverse wind speed component Vx, the longitudinal wind speed component Vy, Vx², Vy², and Vx*Vy are used as features to establish a regression model with the standardized instantaneous displacement deviation D'.

[0060] The regression model was trained using the random forest algorithm. Random forest is an ensemble learning algorithm composed of multiple decision trees. During training, it randomly samples features and samples, constructs multiple decision trees, and then integrates the results of these decision trees. By training the random forest regression model, we can calculate the importance score of each feature, that is, the degree of influence of each wind speed component and its quadratic and cross terms on the instantaneous displacement deviation.

[0061] The dynamic contribution weights of each wind speed component are determined by combining standardized residual analysis. Standardized residuals refer to the standardized difference between the regression model's predicted value and the actual value. By analyzing standardized residuals, nonlinearities in the regression model can be identified. Gaussian process regression is used to compensate for these nonlinear residuals. Gaussian process regression is a nonparametric regression method that can model arbitrary functions. Modeling and predicting nonlinear residuals using Gaussian process regression yields more accurate predictions.

[0062] The wind speed interference characteristics can be generated by combining the dynamic contribution weights of each wind speed component and the compensated and corrected prediction results.

[0063] Step S133: performing correlation analysis on the atmospheric refractive index distribution data and the continuous trajectory curvature variation in the dynamic positioning feature set, calculating the covariance matrix of the atmospheric refractive index spatial gradient and the trajectory curvature variation, and extracting the principal component feature as the refractive index interference factor to generate the refractive index interference feature.

[0064] In order to analyze the influence of atmospheric refractive index on the trajectory of positioning nodes, it is necessary to perform correlation analysis on the atmospheric refractive index distribution data and the continuous trajectory curvature change in the dynamic positioning feature set.

[0065] First, calculate the spatial gradient of atmospheric refractive index. The spatial gradient of atmospheric refractive index describes the change of atmospheric refractive index in space. The spatial gradient can be calculated by taking partial derivatives of atmospheric refractive index distribution data in different directions. Assuming that the atmospheric refractive index is n(x, y, z), then its gradient components in the x, y, and z directions are respectively 、 、 .

[0066] Then, the covariance matrix of the atmospheric refractive index spatial gradient and the trajectory curvature variation is calculated. The covariance matrix can measure the correlation between the two variables. Assume that the atmospheric refractive index spatial gradient vector is , the trajectory curvature variation vector is C, then the element Cov(Gi, Cj) of the covariance matrix Cov(G, C) represents the covariance between the i-th component of the atmospheric refractive index spatial gradient and the j-th component of the trajectory curvature variation.

[0067] Next, principal component analysis (PCA) is performed on the covariance matrix. PCA is a data dimensionality reduction method that converts high-dimensional data into low-dimensional data while retaining the data's key information. By calculating the eigenvalues ​​and eigenvectors of the covariance matrix, the first k eigenvectors with the largest eigenvalues ​​are selected, and the spatial gradient of the atmospheric refractivity and the change in trajectory curvature are projected onto these eigenvectors to obtain the principal component features.

[0068] These principal component features are used as refractive index interference factors and combined together to generate refractive index interference features.

[0069] Step S134: constructing an environmental interference feature set based on the wind speed interference feature and the refractive index interference feature; wherein the environmental interference feature set includes a wind speed contribution weight set and a refractive index interference factor set, and each positioning node corresponds to a set of wind speed contribution weights and refractive index interference factors; the wind speed contribution weight set is used to quantify the degree of influence on the displacement deviation of the positioning node under different wind speed gradients; the refractive index interference factor set is used to characterize the intensity of the nonlinear effect of the spatial distribution of atmospheric refractive index on the change of trajectory curvature.

[0070] After obtaining the wind speed interference signature and the refractive index interference signature, they are combined to form an environmental interference signature set. The wind speed interference signature includes the wind speed contribution weight, while the refractive index interference signature includes the refractive index interference factor. Combining the wind speed contribution weight and refractive index interference factor corresponding to each positioning node forms the environmental interference signature set.

[0071] Specifically, for each positioning node, a set of wind speed contribution weights is determined based on the wind speed interference signature. These weights are calculated based on a previous analysis of wind speed gradient distribution data and instantaneous displacement deviation. They quantify the impact of wind speed on the displacement deviation of the positioning node under different wind speed gradients. For example, wind speeds of different directions and magnitudes may have varying degrees of impact on the displacement of the positioning node, and the wind speed contribution weights reflect the specific quantitative details of this impact.

[0072] At the same time, for each positioning node, the corresponding refractive index interference factor is extracted from the refractive index interference feature. This factor is obtained by analyzing the correlation between atmospheric refractive index distribution data and the change in continuous trajectory curvature, through steps such as covariance matrix calculation and principal component analysis. It is used to characterize the strength of the nonlinear effect of the spatial distribution of atmospheric refractive index on the change in trajectory curvature of the positioning node. For example, the change in atmospheric refractive index at different spatial locations may cause the propagation path of Beidou satellite signals to bend, thereby affecting the trajectory curvature of the positioning node. The refractive index interference factor reflects the strength of this influence.

[0073] The wind speed contribution weights and refractive index interference factors of all positioning nodes are arranged and combined in order of the positioning nodes to construct an environmental interference feature set. This environmental interference feature set is crucial for subsequent analysis of the impact of environmental factors on the azimuth angle of human shadow operations, as it integrates information on the impact of two important environmental factors, wind speed and atmospheric refractive index, on the motion characteristics of positioning nodes.

[0074] Step S140: determining a real-time azimuth compensation parameter set according to a coupling association relationship between the dynamic positioning feature set and the environmental interference feature set.

[0075] After obtaining the dynamic positioning feature set and the environmental interference feature set, since there is a coupling relationship between the two, that is, environmental factors will affect the dynamic characteristics of the positioning node, and then affect the calculation of the azimuth angle of the human shadow operation, it is necessary to determine the real-time azimuth compensation parameter set based on this coupling relationship to correct the subsequent azimuth angle calculation.

[0076] Step S141: performing weight distribution processing on the instantaneous displacement deviation in the dynamic positioning feature set and the wind speed contribution weight in the environmental interference feature set to generate a first compensation weight coefficient.

[0077] When determining the real-time azimuth compensation parameter set, the instantaneous displacement deviation in the dynamic positioning feature set and the wind speed contribution weight in the environmental interference feature set are first weighted. The instantaneous displacement deviation reflects the position change of the positioning node over a short period of time, while the wind speed contribution weight reflects the degree of wind speed's influence on the positioning node's displacement deviation.

[0078] In order to combine the two, weight distribution is required. A method based on the correlation between the two can be used to distribute the weights. For example, first calculate the correlation coefficient between the instantaneous displacement deviation and the wind speed contribution weight. This correlation coefficient can measure the closeness of the relationship between them. Assuming the correlation coefficient is r, the weight distribution ratio is determined by a function f(r) related to r. Let the weight of the instantaneous displacement deviation be w1, the weight of the wind speed contribution weight be w2, and w1+w2=1. The size of w1 and w2 can be adjusted according to the value of f(r). When r is large, it means that the two are closely related, and it may be necessary to give them relatively balanced weights; when r is small, it may be necessary to adjust the weight distribution according to the specific situation.

[0079] After such weight distribution processing, the instantaneous displacement deviation and the wind speed contribution weight are combined according to their respective weights to generate a first compensation weight coefficient. The first compensation weight coefficient comprehensively considers the impact of the instantaneous displacement deviation and wind speed on the position change of the positioning node.

[0080] Step S142: performing weight distribution processing on the continuous trajectory curvature variation in the dynamic positioning feature set and the refractive index interference factor in the environmental interference feature set to generate a second compensation weight coefficient.

[0081] Similarly, a weighted distribution is performed on the continuous trajectory curvature variation in the dynamic positioning feature set and the refractive index interference factor in the environmental interference feature set. The continuous trajectory curvature variation describes the curvature variation of the positioning node trajectory, while the refractive index interference factor reflects the intensity of the impact of the atmospheric refractive index on the trajectory curvature variation.

[0082] Similarly, first calculate the correlation coefficient between the continuous trajectory curvature change and the refractive index interference factor, denoted as r'. Use a function g(r') related to r' to determine the weight distribution ratio. Let the weight of the continuous trajectory curvature change be w3, and the weight of the refractive index interference factor be w4, with w3 + w4 = 1. Adjust the values ​​of w3 and w4 based on the value of g(r').

[0083] The continuous trajectory curvature change and the refractive index interference factor are combined according to their respective weights to generate a second compensation weight coefficient. This second compensation weight coefficient comprehensively considers the impact of trajectory curvature change and atmospheric refractive index on the positioning node trajectory and is another important parameter for real-time azimuth compensation.

[0084] Step S143: constructing a dynamic compensation weight set according to the first compensation weight coefficient and the second compensation weight coefficient; the dynamic compensation weight set includes a displacement compensation weight and a curvature compensation weight corresponding to each positioning node.

[0085] After obtaining the first and second compensation weight coefficients, they are combined to construct a dynamic compensation weight set. For each positioning node, the first compensation weight coefficient corresponds to the node's displacement compensation weight. This is primarily based on the instantaneous displacement deviation and wind speed contribution weight, reflecting compensation information for the positioning node's displacement. The second compensation weight coefficient corresponds to the node's curvature compensation weight. This is based on the continuous trajectory curvature change and the refractive index interference factor, reflecting compensation information for the positioning node's trajectory curvature.

[0086] The displacement compensation weights and curvature compensation weights of all positioning nodes are arranged and combined in the order of the positioning nodes to form a dynamic compensation weight set. This dynamic compensation weight set provides the basis for accurately compensating the displacement and trajectory curvature of the positioning nodes in the subsequent azimuth calculation.

[0087] Step S144: After normalizing the instantaneous displacement deviation and the continuous trajectory curvature change, the fusion weight coefficients of the normalized results are calculated using the entropy weight method based on the displacement compensation weight and the curvature compensation weight in the dynamic compensation weight set. The normalized results are weightedly fused using the weighted geometric mean algorithm to generate a real-time azimuth compensation parameter set; the real-time azimuth compensation parameter set includes a displacement compensation parameter subset and a curvature compensation parameter subset, each subset is arranged in time series and corresponds one-to-one to a positioning node.

[0088] After constructing the dynamic compensation weight set, to generate the real-time azimuth compensation parameter set, the instantaneous displacement deviation and continuous trajectory curvature change in the dynamic positioning feature set are first standardized. Standardization eliminates the effects of different magnitudes and makes them comparable. For the instantaneous displacement deviation, assuming its mean is μ1 and its standard deviation is σ1, the standardized instantaneous displacement deviation is D'' = (D - μ1) / σ1. For the continuous trajectory curvature change, assuming its mean is μ2 and its standard deviation is σ2, the standardized continuous trajectory curvature change is C'' = (C - μ2) / σ2.

[0089] Then, based on the displacement and curvature compensation weights in the dynamic compensation weight set, the entropy weight method is used to calculate the fusion weight coefficients of each standardized result. The entropy weight method is an objective weighting method that determines weights based on the degree of data dispersion. For the standardized instantaneous displacement deviation D'' and the continuous trajectory curvature change C'', their information entropy is calculated. Information entropy reflects the degree of data disorder. Based on the information entropy, their entropy weights are calculated and denoted as wD and wC, where wD + wC = 1.

[0090] Next, a weighted geometric mean algorithm is used to perform a weighted fusion of the standardized instantaneous displacement deviation D'' and the continuous trajectory curvature change C''. The weighted geometric mean algorithm calculates the geometric mean of each data point according to its corresponding weight. Assuming the fused result is R, then R = (D''^wD)*(C''^wC).

[0091] The fusion results are divided according to positioning nodes and time series, forming displacement compensation parameter subsets and curvature compensation parameter subsets. The displacement compensation parameter subset primarily contains compensation parameters related to the positioning node displacement, while the curvature compensation parameter subset primarily contains compensation parameters related to the curvature of the positioning node trajectory. Each subset is arranged in time series and corresponds to a positioning node one-to-one. Combining these two subsets generates a real-time azimuth compensation parameter set.

[0092] Step S150: performing azimuth iterative solution processing on the Beidou satellite positioning data set based on the real-time azimuth compensation parameter set to generate real-time azimuth parameters of human shadow operations in the target area.

[0093] After obtaining the real-time azimuth compensation parameter set, in order to obtain the real-time azimuth parameters of the human shadow operation in the target area, it is necessary to perform azimuth iterative solution processing on the Beidou satellite positioning data set based on the compensation parameter set.

[0094] Step S151: calling the three-dimensional coordinate sequence in the Beidou satellite positioning data set to construct an initial azimuth angle solution model; the initial azimuth angle solution model uses the least squares method to calculate the azimuth angle reference value.

[0095] First, a three-dimensional coordinate sequence is retrieved from the Beidou satellite positioning data set to construct an initial azimuth angle solution model. The least squares method is a commonly used parameter estimation method that determines the model parameters by minimizing the sum of squared errors. When constructing the initial azimuth angle solution model, the coordinate points in the three-dimensional coordinate sequence are used as observation data, and a functional relationship is assumed between the azimuth angle and these coordinate points. Using the least squares method, the parameters of this functional relationship are estimated to minimize the sum of squared errors between the observation data and the model predictions.

[0096] Specifically, let the coordinate points in the three-dimensional coordinate sequence be (x1, y1, z1), (x2, y2, z2), …, (xn, yn, zn), and assume that the relationship between the azimuth angle α and these coordinate points can be expressed as a function f(x, y, z, α). Using the least squares method, find the value of α that minimizes Σ[f(xi, yi, zi, α) - yi_obs]², where yi_obs is the observed value related to the azimuth angle. The resulting α value is the azimuth reference value. This completes the construction of the initial azimuth angle solution model.

[0097] Step S152: inputting the real-time azimuth compensation parameter set into the initial azimuth solution model to generate a compensated azimuth solution model; the compensated azimuth solution model includes a displacement compensation module and a curvature correction module, which are used to process the displacement compensation parameter subset and the curvature compensation parameter subset respectively.

[0098] After constructing the initial azimuth solution model, the real-time azimuth compensation parameter set is input into the model to generate the compensated azimuth solution model. The real-time azimuth compensation parameter set includes a displacement compensation parameter subset and a curvature compensation parameter subset. To process the information of these two subsets separately, the compensated azimuth solution model includes a displacement compensation module and a curvature correction module.

[0099] The displacement compensation module primarily processes a subset of displacement compensation parameters. In the initial azimuth solution model, the displacement information of the positioning node may not fully account for the impact of environmental factors. The displacement compensation parameter subset contains compensation information that accounts for the effects of environmental factors such as wind speed on the positioning node's displacement. The displacement compensation module combines this subset of displacement compensation parameters with the displacement-related information in the initial azimuth solution model to correct the positioning node's displacement, thereby more accurately calculating the azimuth.

[0100] The curvature correction module primarily processes a subset of curvature compensation parameters. The initial azimuth solution model may not fully account for the impact of environmental factors such as atmospheric refractive index on the trajectory curvature of the positioning node. The curvature compensation parameter subset contains compensation information to account for the impact of these environmental factors on trajectory curvature. The curvature correction module combines this subset of curvature compensation parameters with the trajectory curvature information from the initial azimuth solution model to correct the trajectory curvature of the positioning node, thereby improving the accuracy of the azimuth calculation.

[0101] Through the function of these two modules, the compensated azimuth angle solution model can more comprehensively consider the impact of environmental factors on the movement of positioning nodes, thereby more accurately calculating the azimuth angle of human shadow operation.

[0102] Step S153: Iteratively optimize the compensated azimuth angle calculation model using a gradient descent algorithm until the compensated azimuth angle calculation model converges, extracting azimuth angle parameters to generate real-time azimuth angle parameters for shadow operation.

[0103] After generating the compensated azimuth angle solution model, in order to obtain more accurate real-time azimuth angle parameters for human shadow operation, it is necessary to iteratively optimize the model using the gradient descent algorithm.

[0104] Step S1531: Obtain the azimuth angle solution error and error gradient direction of the current iteration cycle.

[0105] In each iteration, the azimuth solution error and error gradient direction must first be obtained. To obtain the azimuth solution error, it is necessary to extract the measured azimuth data of the human shadow operation during the current iteration. This measured data is collected by the angle measurement equipment of the ground reference station and must be synchronized with the Beidou positioning data to accurately compare the model calculation results with the actual measurement results.

[0106] At the same time, the compensated azimuth solution model is called to generate the azimuth prediction data for the current iteration cycle. The azimuth prediction data includes the azimuth sequence output by the compensated azimuth solution model and its confidence interval. The measured azimuth data and the azimuth prediction data are time-aligned, and the absolute and relative errors are calculated point by point. The azimuth solution error is obtained by comprehensively calculating these errors.

[0107] Based on the azimuth solution error, backpropagation is performed on the compensation weight coefficients in the real-time azimuth compensation parameter set. Backpropagation involves constructing a computational graph and performing layer-by-layer differentiation along the topology of the compensated azimuth solution model. This process determines the gradient of each compensation weight coefficient's contribution to the total error. This gradient direction is the error gradient direction. The error gradient direction indicates in which direction in parameter space, adjusting the compensation weight coefficients will reduce the azimuth solution error.

[0108] Step S1532: Calculate the product of the gradient direction and the learning rate according to the error gradient direction, and incrementally update the compensation weight coefficient to adjust the compensation weight coefficient in the real-time azimuth compensation parameter set.

[0109] After obtaining the error gradient direction, in order to adjust the compensation weight coefficient in the real-time azimuth compensation parameter set, it is necessary to calculate the product of the gradient direction and the learning rate based on the error gradient direction. First, the error gradient direction is normalized. Normalization involves calculating the L2 norm of the gradient vector and dividing each gradient component by the norm value. This ensures that the length of the gradient direction is 1, which facilitates subsequent calculations.

[0110] Assuming the normalized error gradient direction vector is G' and the learning rate is η, then the gradient adjustment step size is ΔG = η * G'. The first compensation weight coefficient and the second compensation weight coefficient are updated based on the gradient adjustment step size. Assuming the first compensation weight coefficient is w1 and the second compensation weight coefficient is w2, the updated first compensation weight coefficient is w1' = w1 + ΔG1, and the updated second compensation weight coefficient is w2' = w2 + ΔG2, where ΔG1 and ΔG2 are the corresponding components of the gradient adjustment step size vector ΔG.

[0111] The updated first and second compensation weight coefficients are rewritten into the dynamic compensation weight set. A spatial regularization constraint is introduced during the update to ensure that the difference in compensation weight coefficients between adjacent positioning nodes does not exceed a set threshold. This ensures the spatial smoothness of the compensation weight coefficients. At the same time, the weight records for the most recent N iterations are retained, with the value of N being inversely proportional to the spatial distribution density of the positioning nodes. That is, the lower the spatial distribution density of the positioning nodes, the smaller the value of N. Through this update operation, the compensation weight coefficients in the real-time azimuth compensation parameter set are adjusted, preparing for the next iterative optimization.

[0112] Step S1533: recalculating the outputs of the displacement compensation module and the curvature correction module based on the adjusted compensation weight coefficient, and synchronously updating the internal state variables of the model to update the compensated azimuth angle solution model.

[0113] After adjusting the compensation weight coefficients in the real-time azimuth compensation parameter set, the outputs of the displacement compensation module and the curvature correction module are recalculated based on these adjusted compensation weight coefficients. The displacement compensation module re-corrects the displacement information of the positioning node based on the adjusted displacement compensation weight coefficients and displacement compensation parameter subset; the curvature correction module re-corrects the trajectory curvature information of the positioning node based on the adjusted curvature compensation weight coefficients and curvature compensation parameter subset.

[0114] At the same time, the model's internal state variables are updated synchronously. These state variables record the intermediate results and parameter information of the model during the calculation process. By updating these state variables, the compensated azimuth solution model can reflect the changes after the compensation weight coefficient adjustment, thereby updating the entire model.

[0115] Step S1534: Repeat the above operations until the compensated azimuth angle solution model meets the convergence conditions, extract the azimuth angle parameters in the converged azimuth angle solution model, and generate real-time azimuth angle parameters for human shadow operation. The convergence conditions include that the error change rate of the azimuth angle solution error is less than the set threshold or the maximum number of iterations is reached.

[0116] Repeat steps S1531-S1533 continuously, adjusting the compensation weight coefficients with each iteration and updating the compensated azimuth solution model to gradually reduce the azimuth solution error. During the iteration process, determine whether the compensated azimuth solution model meets the convergence conditions. Convergence conditions include the error rate of the azimuth solution error being less than a set threshold, meaning the error decreases very slowly, indicating that the model is close to the optimal solution; or reaching the maximum number of iterations. To avoid infinite iterations, a maximum number of iterations is set, and iterations are terminated when this maximum number of iterations is reached.

[0117] When the compensated azimuth solution model meets the convergence criteria, the azimuth parameters in the converged azimuth solution model are extracted. These azimuth parameters are obtained through multiple iterations of optimization and can more accurately reflect the real-time azimuth of human shadow operations in the target area. By combining these azimuth parameters, the real-time azimuth parameters of human shadow operations can be generated.

[0118] Step S1534-11: Obtain the parameter space distribution of the converged azimuth solution model, where the parameter space distribution includes the weight matrix of the displacement compensation module and the bias vector of the curvature correction module.

[0119] After obtaining the converged azimuth solution model, its parameter space distribution needs to be obtained. The parameter space distribution reflects the values ​​of each parameter in the model. For the compensated azimuth solution model, this mainly includes the weight matrix of the displacement compensation module and the bias vector of the curvature correction module.

[0120] The displacement compensation module's weight matrix records the weight distribution of each parameter when processing a subset of displacement compensation parameters. The size of the elements in this weight matrix reflects the degree to which different displacement compensation parameters affect the final azimuth calculation. For example, a larger element in the weight matrix indicates that the corresponding displacement compensation parameter plays a more important role in the azimuth calculation.

[0121] The curvature correction module's offset vector records the bias information used to adjust the model output when processing a subset of curvature compensation parameters. Each component of the offset vector corresponds to a dimension of the curvature correction module's output. By adjusting the offset vector, the curvature correction module's output can be shifted and adjusted to better fit the actual azimuth angle.

[0122] Step S1534-12: Calculate the eigenvectors of the parameter covariance matrix of the parameter space distribution, select the first K principal components by sorting the eigenvector values, and extract the azimuth principal component features.

[0123] After obtaining the parameter space distribution, the parameter covariance matrix is ​​calculated. The parameter covariance matrix reflects the correlation between the parameters in the parameter space. By calculating the parameter covariance matrix, we can understand the mutual influence relationship between different parameters.

[0124] Next, we calculate the eigenvectors of the parameter covariance matrix. An eigenvector is a special type of matrix vector that only scales, not rotates, under matrix transformations. Each eigenvector corresponds to an eigenvalue, the magnitude of which reflects the degree of data dispersion in the direction represented by the eigenvector.

[0125] The eigenvectors are sorted by their values ​​and the principal components corresponding to the first K eigenvectors are selected. Principal component analysis is a data dimensionality reduction method that can reduce the high-dimensional parameter space to K dimensions by selecting the first K principal components while retaining the key information of the data.

[0126] The parameters of the converged azimuth solution model are projected onto these first K principal components. The resulting projections constitute the azimuth principal component features. These principal component features represent a dimensionality-reduced representation of the original parameter space. They retain the most important information in the parameter space and can more concisely describe the azimuth characteristics while reducing the dimensionality of the data and the complexity of subsequent processing.

[0127] Step S1534-13: Use the kernel density estimation method to perform non-parametric probability modeling on the principal component features and construct a probability density function of the azimuth parameters.

[0128] After obtaining the azimuth principal component characteristics, in order to further analyze the distribution of azimuth parameters, the kernel density estimation method is used to perform non-parametric probability modeling on the principal component characteristics and construct the probability density function of the azimuth parameters.

[0129] Kernel density estimation is a non-parametric estimation method that does not make assumptions about the data distribution. Instead, it estimates the probability density of the data using a kernel function. First, kernel density estimation is performed on the principal component features, using a Gaussian kernel function. The Gaussian kernel function is smooth and symmetric, and can fit the data distribution well. When using the Gaussian kernel function for kernel density estimation, it is necessary to determine the bandwidth parameter. The bandwidth parameter controls the width of the kernel function, and its size affects the smoothness of the estimated results. The bandwidth parameter is adaptively adjusted based on the variance of the principal component features. When the variance is large, the bandwidth parameter is appropriately increased to better cover the data; when the variance is small, the bandwidth parameter is appropriately reduced to more precisely characterize the local characteristics of the data.

[0130] After kernel density estimation, an initial probability density distribution is generated. However, this initial probability density distribution may deviate from the actual azimuth angle because it is estimated only based on principal component features. To improve the accuracy of the probability density function, the initial probability density distribution is calibrated based on the three-dimensional coordinate sequence in the Beidou satellite positioning data set. The three-dimensional coordinate sequence contains the actual position information of the positioning node. By associating and adjusting the initial probability density distribution with the three-dimensional coordinate sequence, the probability density function can better reflect the actual azimuth angle distribution, thereby generating a calibrated probability density function.

[0131] Step S1534-14: Determine the final real-time azimuth angle parameter of the human shadow operation according to the extreme point of the probability density function, wherein the extreme point is obtained by traversing the domain of the probability density function and finding the local maximum point.

[0132] After constructing the probability density function for the azimuth parameters, the final real-time azimuth parameters for the human shadow operation are determined based on the function's extreme points. Extreme points are the points within the domain where the function reaches its maximum or minimum value. For a probability density function, the extreme points typically represent locations with a high probability of azimuth occurring.

[0133] To find the extreme values ​​of a probability density function, you need to traverse the domain of the probability density function. The domain is the range of values ​​that the probability density function can take. Within this range, the value of the probability density function is calculated point by point, and the function values ​​of adjacent points are compared. When the function value of a point is greater than that of its adjacent points, that point is a local maximum. By continuously traversing and comparing, all local maximum points can be found.

[0134] Among these local maximum points, the point that best represents the actual azimuth angle is selected as the final real-time azimuth angle parameter for the human shadow operation. This selection process can be based on specific application requirements and actual conditions. For example, the local maximum point with the highest probability density can be selected because it indicates that the azimuth angle has the highest probability of occurring and is most likely to be the actual human shadow operation azimuth angle.

[0135] Step S1534-141: perform kernel density estimation processing on the azimuth principal component feature to generate an initial probability density distribution; the kernel density estimation processing adopts a Gaussian kernel function, and the bandwidth parameter is adaptively adjusted according to the variance of the principal component feature.

[0136] When performing kernel density estimation, a Gaussian kernel function is used to estimate the probability density of the azimuth principal component features. The Gaussian kernel function is based on the normal distribution and can fit the data in a smooth manner. For each principal component feature data point, the Gaussian kernel function is used to calculate its contribution to the probability density of the surrounding area with that point as the center.

[0137] When using the Gaussian kernel function, the choice of the bandwidth parameter is crucial. The bandwidth parameter determines the width of the Gaussian kernel function, that is, the extent to which a data point influences the surrounding area. Here, the bandwidth parameter is adaptively adjusted based on the variance of the principal component features. The variance reflects the degree of dispersion of the principal component feature data. A large variance indicates that the data is dispersed, and a larger bandwidth parameter is required. This allows the Gaussian kernel function to cover a wider area, thereby more comprehensively estimating the probability density. A small variance indicates that the data is concentrated, and a smaller bandwidth parameter is required to more accurately characterize the local characteristics of the data.

[0138] By applying the Gaussian kernel function to each principal component feature data point and summing up the contributions of all data points, we obtain an initial probability density distribution, which describes the probability of the azimuth principal component feature at different values.

[0139] Step S1534-142: Calibrate the initial probability density distribution according to the three-dimensional coordinate sequence in the Beidou satellite positioning data set to generate a calibrated probability density function.

[0140] The initial probability density distribution is based on the principal component characteristics of the azimuth angle, which may deviate from the actual azimuth angle. To make the probability density function more accurately reflect the actual situation, it is necessary to calibrate the initial probability density distribution based on the 3D coordinate sequence in the Beidou satellite positioning data set.

[0141] The 3D coordinate sequence contains the actual location information of the positioning node, which is closely related to the azimuth angle. First, the relationship between each coordinate point in the 3D coordinate sequence and the azimuth angle is analyzed to determine the azimuth angle range corresponding to each coordinate point. Then, the initial probability density distribution is compared with these azimuth angle ranges and adjusted.

[0142] For example, if the frequency of positioning nodes in a certain area is high in the 3D coordinate sequence, it means that the actual probability of the azimuth angle corresponding to that area is high. In this case, the probability value within the corresponding azimuth angle range in the initial probability density distribution needs to be appropriately increased. Conversely, if the frequency of positioning nodes in a certain area is low, the probability value within the corresponding azimuth angle range needs to be appropriately reduced.

[0143] By adjusting the initial probability density distribution according to the three-dimensional coordinate sequence, the probability density function can better reflect the actual azimuth distribution, thereby generating a calibrated probability density function.

[0144] Step S1534-143: Extract the peak area of ​​the probability density function after calibration, and use the median point of the peak area as the real-time azimuth parameter of the human shadow operation; wherein, the peak area is obtained by binarizing the probability density function by setting a density threshold, and the median point is calculated using the regional centroid algorithm.

[0145] After obtaining the calibrated probability density function, in order to determine the final real-time azimuth parameters of the human shadow operation, it is necessary to extract the peak area of ​​the calibrated probability density function. The peak area represents the area with a high probability of azimuth occurrence, which usually contains the most likely azimuth of the human shadow operation.

[0146] The peak region is obtained by binarizing the probability density function by setting a density threshold. The density threshold is a pre-set probability density value. When the value of the probability density function is greater than the threshold, the point is marked as belonging to the peak region; when the value of the probability density function is less than the threshold, the point is marked as not belonging to the peak region. Through this binary segmentation method, the calibrated probability density function is divided into peak and non-peak regions.

[0147] For the extracted peak area, calculate the median point of the peak area, and use the median point as the real-time azimuth parameter of the human shadow operation. The median point is calculated using the regional centroid algorithm. The regional centroid algorithm is a method for calculating the center position of the region, which takes into account the position and weight of each point in the region. In the peak area, the weight of each point can be determined according to its corresponding probability density value. The larger the probability density value, the greater the weight of the point. By performing a weighted average calculation on all points in the peak area, the result is the center of mass of the peak area, that is, the median point. The median point represents the center position of the peak area. Using it as the real-time azimuth parameter of the human shadow operation can more accurately reflect the actual operation azimuth situation.

[0148] Throughout the entire process, data collection, including Beidou satellite positioning data and meteorological monitoring data, strictly complies with relevant laws and regulations to ensure data legitimacy and compliance. Data encryption technology is employed to protect potentially privacy-sensitive data and prevent leakage. During data transmission, secure encryption protocols are used to prevent data theft or tampering during transmission. A secure storage system is used to encrypt data and implement strict access control, ensuring that only authorized personnel can access and process data.

[0149] Figure 2A schematic diagram illustrates exemplary hardware and software components of a BeiDou satellite-based real-time azimuth angle calculation system 100 for shadow work, which can implement the concepts of the present application, according to some embodiments of the present application. For example, processor 120 can be used in BeiDou satellite-based real-time azimuth angle calculation system 100 for shadow work and perform the functions described in the present application.

[0150] The BeiDou satellite-based system 100 for calculating the azimuth angle of a shadow operation can be a general-purpose server or a special-purpose server, both of which can be used to implement the BeiDou satellite-based method for calculating the azimuth angle of a shadow operation. Although only one server is shown in this application, for convenience, the functions described in this application can be implemented in a distributed manner on multiple similar platforms to balance the processing load.

[0151] For example, the Beidou satellite-based real-time azimuth angle calculation system 100 for human shadow operations can include a network port 110 connected to a network, one or more processors 120 for executing program instructions, a communication bus 130, and various forms of storage media 140, such as a disk, ROM, or RAM, or any combination thereof. For example, the Beidou satellite-based real-time azimuth angle calculation system 100 for human shadow operations can also include program instructions stored in ROM, RAM, or other types of non-transitory storage media, or any combination thereof. The method of the present application can be implemented based on these program instructions. The Beidou satellite-based real-time azimuth angle calculation system 100 also includes an I / O interface 150 between the computer and other input and output devices.

[0152] For ease of explanation, only one processor is described in the Beidou satellite-based real-time calculation system 100 for the azimuth angle of human shadow operations. However, it should be noted that the Beidou satellite-based real-time calculation system 100 for the azimuth angle of human shadow operations in this application can also include multiple processors, so the steps performed by one processor described in this application can also be performed jointly or individually by multiple processors. For example, if the processor of the Beidou satellite-based real-time calculation system 100 for the azimuth angle of human shadow operations executes steps A and B, it should be understood that steps A and B can also be performed jointly by two different processors or individually in one processor. For example, the first processor executes step A and the second processor executes step B, or the first processor and the second processor execute steps A and B together.

[0153] In addition, an embodiment of the present invention also provides a readable storage medium, which has computer-executable instructions preset in the readable storage medium. When the processor executes the computer-executable instructions, the above-mentioned method for real-time calculation of the azimuth angle of human shadow operation based on Beidou satellite is implemented.

[0154] It should be noted that in order to simplify the description of the present invention and thus help understand one or more embodiments of the invention, in the foregoing description of the embodiments of the present invention, multiple features are sometimes combined into one embodiment, figure or description thereof.

Claims

1. A real-time calculation method for the azimuth angle of human shadow operation based on Beidou satellite, characterized in that: The method comprises: Acquire a Beidou satellite positioning data set for a target area, the Beidou satellite positioning data set comprising a three-dimensional coordinate sequence of multiple positioning nodes within a continuous time period; Performing dynamic trajectory fitting processing on the Beidou satellite positioning data set to obtain a dynamic positioning feature set including an instantaneous displacement deviation and a continuous trajectory curvature change; Performing environmental interference feature analysis processing on the target area based on the dynamic positioning feature set to generate an environmental interference feature set; Determining a real-time azimuth compensation parameter set according to a coupling association relationship between the dynamic positioning feature set and the environmental interference feature set; An azimuth iterative solution process is performed on the Beidou satellite positioning data set based on the real-time azimuth compensation parameter set to generate a real-time azimuth parameter of the human shadow operation in the target area.

2. The method for real-time calculation of azimuth angle of shadow operation based on Beidou satellite according to claim 1 is characterized in that: The dynamic trajectory fitting process is performed on the Beidou satellite positioning data set to obtain a dynamic positioning feature set including an instantaneous displacement deviation and a continuous trajectory curvature change, including: Extracting spatiotemporal distribution information of positioning nodes in the three-dimensional coordinate sequence to generate spatiotemporal correlation features; Adaptively segmenting the spatiotemporal correlation features using a sliding window to obtain multiple windowed trajectory segments. The time span of the sliding window is dynamically adjusted according to the spatial distribution density of the positioning nodes. The lower the spatial density, the smaller the sliding window time span. Kalman filtering is used to smooth the noise of discrete coordinate points in the sliding window. Adopting an adaptive piecewise interpolation method to perform nonlinear reconstruction processing on the discrete coordinate points in the windowed trajectory segments to generate a smooth trajectory curve; Calculating the displacement deviation between adjacent positioning nodes based on the smooth trajectory curve to generate the instantaneous displacement deviation; Calculating the curvature derivatives of the smooth trajectory curve at adjacent interpolation points, determining the direction of curvature change according to the change in the sign of the curvature derivative, and extracting the curvature change of the continuous trajectory; The instantaneous displacement deviation and the continuous trajectory curvature change are combined into a dynamic positioning feature set.

3. The method for real-time calculation of azimuth angle of shadow operation based on Beidou satellite according to claim 1 is characterized in that: The performing environmental interference feature analysis processing on the target area based on the dynamic positioning feature set to generate the environmental interference feature set includes: Acquire a meteorological monitoring data set of a target area, wherein the meteorological monitoring data set includes wind speed gradient distribution data and atmospheric refractive index distribution data; The wind speed gradient distribution data is vector-decomposed to obtain transverse and longitudinal wind speed components. After the instantaneous displacement deviation is standardized, a nonlinear regression model of the instantaneous displacement deviation and the wind speed component is established. A feature interaction matrix is ​​constructed by introducing quadratic terms and cross terms of the wind speed components. The regression model is trained using a random forest algorithm and feature importance scores are calculated. The dynamic contribution weight of each wind speed component is determined in combination with standardized residual analysis. The nonlinear residual part is compensated and corrected using Gaussian process regression to generate a wind speed interference feature. Performing correlation analysis on the atmospheric refractive index distribution data and the continuous trajectory curvature variation in the dynamic positioning feature set, calculating the covariance matrix of the atmospheric refractive index spatial gradient and the trajectory curvature variation, and extracting the principal component feature as the refractive index interference factor to generate a refractive index interference feature; Constructing an environmental interference feature set according to the wind speed interference feature and the refractive index interference feature; Among them, the environmental interference feature set includes a wind speed contribution weight set and a refractive index interference factor set, and each positioning node corresponds to a set of wind speed contribution weights and refractive index interference factors; the wind speed contribution weight set is used to quantify the degree of influence on the positioning node displacement deviation under different wind speed gradients; the refractive index interference factor set is used to characterize the intensity of the nonlinear effect of the spatial distribution of atmospheric refractive index on the change of trajectory curvature.

4. The method for real-time calculation of azimuth angle of shadow operation based on Beidou satellite according to claim 1 is characterized in that: The determining of a real-time azimuth compensation parameter set according to a coupling association relationship between the dynamic positioning feature set and the environmental interference feature set includes: Performing weight distribution processing on the instantaneous displacement deviation in the dynamic positioning feature set and the wind speed contribution weight in the environmental interference feature set to generate a first compensation weight coefficient; Performing weight distribution processing on the continuous trajectory curvature change in the dynamic positioning feature set and the refractive index interference factor in the environmental interference feature set to generate a second compensation weight coefficient; Constructing a dynamic compensation weight set according to the first compensation weight coefficient and the second compensation weight coefficient; the dynamic compensation weight set includes a displacement compensation weight and a curvature compensation weight corresponding to each positioning node; After normalizing the instantaneous displacement deviation and the continuous trajectory curvature change respectively, the fusion weight coefficient of each normalized result is calculated by the entropy weight method based on the displacement compensation weight and the curvature compensation weight in the dynamic compensation weight set, and the weighted fusion processing of each normalized result is performed using the weighted geometric mean algorithm to generate a real-time azimuth compensation parameter set; the real-time azimuth compensation parameter set includes a displacement compensation parameter subset and a curvature compensation parameter subset, each subset is arranged in time series and corresponds one-to-one to a positioning node.

5. The method for real-time calculation of azimuth angle of shadow operation based on Beidou satellite according to claim 4 is characterized in that: The performing azimuth iterative solution processing on the Beidou satellite positioning data set based on the real-time azimuth compensation parameter set to generate real-time azimuth parameters of human shadow operations in the target area includes: Calling the three-dimensional coordinate sequence in the Beidou satellite positioning data set to construct an initial azimuth angle solution model; the initial azimuth angle solution model uses the least squares method to calculate the azimuth angle reference value; Inputting the real-time azimuth compensation parameter set into the initial azimuth solution model to generate a compensated azimuth solution model; the compensated azimuth solution model includes a displacement compensation module and a curvature correction module, which are respectively used to process the displacement compensation parameter subset and the curvature compensation parameter subset; The compensated azimuth angle solution model is iteratively optimized using a gradient descent algorithm until the compensated azimuth angle solution model converges, and then azimuth angle parameters are extracted to generate real-time azimuth angle parameters for human shadow operation.

6. The method for real-time calculation of azimuth angle of shadow operation based on Beidou satellite according to claim 5 is characterized in that: The method of iteratively optimizing the compensated azimuth angle calculation model using a gradient descent algorithm until the compensated azimuth angle calculation model converges and extracting azimuth angle parameters to generate real-time azimuth angle parameters for shadow operation includes: Get the azimuth angle solution error and error gradient direction of the current iteration cycle; Calculating the product of the gradient direction and the learning rate according to the error gradient direction, and incrementally updating the compensation weight coefficient to adjust the compensation weight coefficient in the real-time azimuth compensation parameter set; Recalculating the outputs of the displacement compensation module and the curvature correction module based on the adjusted compensation weight coefficient, and synchronously updating the internal state variables of the model to update the compensated azimuth angle solution model; Repeat the above operations until the compensated azimuth angle solution model meets the convergence conditions, extract the azimuth angle parameters in the converged azimuth angle solution model, and generate the real-time azimuth angle parameters of the human shadow operation. The convergence conditions include that the error change rate of the azimuth angle solution error is less than the set threshold or reaches the maximum number of iterations.

7. The method for real-time calculation of azimuth angle of shadow operation based on Beidou satellite according to claim 6 is characterized in that: The obtaining of the azimuth angle solution error and the error gradient direction of the current iteration cycle includes: Extracting the measured azimuth angle data of the human shadow operation in the current iteration cycle, wherein the measured azimuth angle data of the human shadow operation is collected by the angle measurement equipment of the ground reference station and is synchronized with the Beidou positioning data; Calling the compensated azimuth angle solution model to generate azimuth angle prediction data for a current iteration cycle, the azimuth angle prediction data including an azimuth angle sequence output by the compensated azimuth angle solution model and its confidence interval; Performing time alignment processing on the azimuth measured data and the azimuth predicted data, and calculating the absolute error and relative error point by point to obtain the azimuth solution error; Based on the azimuth angle solution error, the compensation weight coefficients in the real-time azimuth angle compensation parameter set are back-propagated to generate the error gradient direction; wherein, the back-propagation processing includes constructing a calculation graph and deriving the topology structure of the compensated azimuth angle solution model layer by layer to obtain the contribution gradient of each compensation weight coefficient to the total error.

8. The method for real-time calculation of azimuth angle of shadow operation based on Beidou satellite according to claim 6 is characterized in that: The calculating the product of the gradient direction and the learning rate according to the error gradient direction, and incrementally updating the compensation weight coefficient to adjust the compensation weight coefficient in the real-time azimuth compensation parameter set, includes: Normalizing the error gradient direction to generate a gradient adjustment step size, wherein the normalization includes calculating the L2 norm of the gradient vector and dividing each gradient component by the norm value; updating the first compensation weight coefficient and the second compensation weight coefficient according to the gradient adjustment step size; The updated first compensation weight coefficient and the second compensation weight coefficient are rewritten into the dynamic compensation weight set. A spatial regularization constraint is introduced during the update so that the difference in the compensation weight coefficients of adjacent positioning nodes does not exceed a set threshold, and the weight records of the most recent N iterations are retained; wherein the value of N is inversely proportional to the spatial distribution density of the positioning nodes.

9. The method for real-time calculation of azimuth angle of shadow operation based on Beidou satellite according to claim 6, characterized in that: The step of extracting the azimuth parameters from the converged azimuth solution model to generate real-time azimuth parameters for the shadow operation includes: Obtaining a parameter space distribution of the converged azimuth solution model, wherein the parameter space distribution includes a weight matrix of a displacement compensation module and a bias vector of a curvature correction module; Calculating the eigenvectors of the parameter covariance matrix of the parameter space distribution, and selecting the first K principal components in order of the eigenvector values ​​to extract the azimuth principal component features; A kernel density estimation method is used to perform non-parametric probability modeling on the principal component characteristics to construct a probability density function of the azimuth parameter; The final real-time azimuth angle parameter of the human shadow operation is determined according to the extreme value point of the probability density function, and the extreme value point is obtained by traversing the definition domain of the probability density function and finding the local maximum point.

10. The method for real-time calculation of azimuth angle of shadow operation based on Beidou satellite according to claim 9, characterized in that: The method of using a kernel density estimation method to perform non-parametric probability modeling on the principal component features and constructing a probability density function of the azimuth parameter includes: Performing kernel density estimation on the azimuth principal component feature to generate an initial probability density distribution; the kernel density estimation process uses a Gaussian kernel function, and the bandwidth parameter is adaptively adjusted according to the variance of the principal component feature; Calibrate the initial probability density distribution according to the three-dimensional coordinate sequence in the Beidou satellite positioning data set to generate a calibrated probability density function; The peak area of ​​the calibrated probability density function is extracted, and the median point of the peak area is used as the real-time azimuth parameter of the human shadow operation; wherein, the peak area is obtained by binarizing the probability density function by setting a density threshold, and the median point is calculated using the regional centroid algorithm.

Images