Real-time calculation method of azimuth angle for weather modification operations based on BeiDou satellite

Abstract

The present invention provides a real-time calculation method for the azimuth of a human shadow operation based on a Beidou satellite. The method comprises the following steps: firstly, obtaining a Beidou satellite positioning data set of a target area, comprising a three-dimensional coordinate sequence of multiple positioning nodes within a continuous time period; then, performing dynamic trajectory fitting processing on the Beidou satellite positioning data set to obtain a dynamic positioning feature set; performing environmental interference feature analysis based on the dynamic positioning feature set to generate an environmental interference feature set; determining a real-time azimuth compensation parameter set based on a coupling correlation relationship between the dynamic positioning feature and the environmental interference feature; finally, performing azimuth iterative solution processing on the positioning data set based on the real-time azimuth compensation parameter set to generate real-time azimuth parameters for the human shadow operation in the target area, thereby achieving real-time and accurate calculation of the azimuth of the human shadow operation and improving the operation effect and safety.

Description

Technical Field

[0001] This invention relates to the fields of meteorological technology and satellite navigation and positioning technology, and more specifically, to a method for real-time calculation of the azimuth angle of shadowing operations based on the BeiDou satellite. Background Technology

[0002] In the field of meteorology, weather modification operations (artificial weather modification), such as artificial rain enhancement and hail suppression, are of great significance for agricultural production, water resource management, and disaster prevention. The precise implementation of artificial weather modification operations relies on accurate calculation of the operational azimuth angle to ensure maximum effectiveness. Traditional methods for calculating the azimuth angle of artificial weather modification operations mainly depend on ground observation equipment and empirical formulas, which can provide relatively accurate results in static or low-dynamic environments. However, in complex and variable natural environments, such as mountainous areas, oceans, or regions with drastic weather changes, ground observation equipment may be affected by various factors such as terrain obstruction and meteorological interference, leading to significant errors in azimuth angle calculations. Furthermore, with technological advancements, the requirements for the real-time performance and accuracy of operational azimuth angles are becoming increasingly stringent, and traditional methods are no longer sufficient to meet these demands. Therefore, developing a method capable of calculating the azimuth angle of artificial weather modification operations in real time and accurately has become an urgent problem to be solved. Summary of the Invention

[0003] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a method for real-time calculation of the azimuth angle of cloud seeding operations based on the BeiDou satellite system, the method comprising:

[0004] Acquire a set of BeiDou satellite positioning data for the target area, wherein the set of BeiDou satellite positioning data contains a sequence of three-dimensional coordinates of multiple positioning nodes within a continuous time period;

[0005] The BeiDou satellite positioning data set is subjected to dynamic trajectory fitting processing to obtain a dynamic positioning feature set that includes instantaneous displacement deviation and continuous trajectory curvature change.

[0006] Based on the dynamic positioning feature set, environmental interference feature analysis processing is performed on the target area to generate an environmental interference feature set.

[0007] Based on the coupling relationship between the dynamic positioning feature set and the environmental interference feature set, the real-time azimuth compensation parameter set is determined;

[0008] Based on the real-time azimuth compensation parameter set, the azimuth iterative calculation of the BeiDou satellite positioning data set is performed to generate real-time azimuth parameters for cloud seeding operations in the target area.

[0009] In one possible implementation of the first aspect, the dynamic trajectory fitting processing of the BeiDou satellite positioning data set to obtain a dynamic positioning feature set including instantaneous displacement deviation and continuous trajectory curvature change includes:

[0010] Extract the spatiotemporal distribution information of the positioning nodes in the three-dimensional coordinate sequence to generate spatiotemporal correlation features;

[0011] The spatiotemporal correlation features are subjected to adaptive sliding window segmentation 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 time span of the sliding window. Kalman filtering is used to smooth the noise of discrete coordinate points within the sliding window.

[0012] An adaptive piecewise interpolation method is used to perform nonlinear reconstruction processing on the discrete coordinate points in the windowed trajectory segment to generate a smooth trajectory curve.

[0013] The displacement deviation between adjacent positioning nodes is calculated based on the smooth trajectory curve, and the instantaneous displacement deviation is generated.

[0014] Calculate the curvature derivative of the smooth trajectory curve at adjacent interpolation points, determine the direction of curvature change based on the change in the sign of the curvature derivative, and extract the curvature change of the continuous trajectory.

[0015] The instantaneous displacement deviation and the continuous trajectory curvature change are combined into a dynamic positioning feature set.

[0016] In one possible implementation of the first aspect, the step of performing environmental interference feature analysis on the target area based on the dynamic positioning feature set to generate an environmental interference feature set includes:

[0017] Acquire a set of meteorological monitoring data for the target area, wherein the set of meteorological monitoring data includes wind speed gradient distribution data and atmospheric refractive index distribution data;

[0018] The wind speed gradient distribution data is vector decomposed to obtain the lateral wind speed component and the longitudinal wind speed component. After standardizing the instantaneous displacement deviation, a nonlinear regression model of instantaneous displacement deviation and wind speed component is established. By introducing the quadratic term and cross term of the wind speed component, a feature interaction matrix is ​​constructed. The regression model is trained using the random forest algorithm and the feature importance score is calculated. Combined with standardized residual analysis, the dynamic contribution weight of each wind speed component is determined. The nonlinear residual part is compensated and corrected using Gaussian process regression to generate wind speed interference features.

[0019] The atmospheric refractive index distribution data and the continuous trajectory curvature change in the dynamic positioning feature set are correlated and analyzed. The covariance matrix of the atmospheric refractive index spatial gradient and the trajectory curvature change is calculated, and the principal component features are extracted as refractive index interference factors to generate refractive index interference features.

[0020] Construct an environmental interference feature set based on the wind speed interference features and the refractive index interference features;

[0021] The environmental interference feature set includes a wind speed contribution weight set and a refractive index interference factor set. 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 of different wind speed gradients on the displacement deviation of the positioning node. The refractive index interference factor set is used to characterize the nonlinear effect of the spatial distribution of atmospheric refractive index on the change of trajectory curvature.

[0022] In one possible implementation of the first aspect, determining the real-time azimuth compensation parameter set based on the coupling relationship between the dynamic positioning feature set and the environmental interference feature set includes:

[0023] The instantaneous displacement deviation in the dynamic positioning feature set and the wind speed contribution weight in the environmental disturbance feature set are weighted and processed to generate a first compensation weight coefficient.

[0024] The continuous trajectory curvature change in the dynamic positioning feature set and the refractive index interference factor in the environmental interference feature set are weighted and processed to generate a second compensation weight coefficient.

[0025] A dynamic compensation weight set is constructed based on the first compensation weight coefficient and the second compensation weight coefficient; the dynamic compensation weight set includes the displacement compensation weight and curvature compensation weight corresponding to each positioning node;

[0026] After standardizing the instantaneous displacement deviation and the continuous trajectory curvature change, the fusion weight coefficient of each standardized result is calculated using the entropy weight method based on the displacement compensation weight and curvature compensation weight in the dynamic compensation weight set. The weighted geometric average algorithm is then used to perform weighted fusion processing on each standardized result to generate a real-time azimuth compensation parameter set. The real-time azimuth compensation parameter set includes a subset of displacement compensation parameters and a subset of curvature compensation parameters. Each subset is arranged in a time sequence and corresponds one-to-one with the positioning node.

[0027] In one possible implementation of the first aspect, the step of performing azimuth iterative calculation processing on the BeiDou satellite positioning data set based on the real-time azimuth compensation parameter set to generate real-time azimuth parameters for cloud seeding operations in the target area includes:

[0028] An initial azimuth angle calculation model is constructed by calling the three-dimensional coordinate sequence in the BeiDou satellite positioning data set; the initial azimuth angle calculation model uses the least squares method to calculate the azimuth angle reference value.

[0029] The real-time azimuth compensation parameter set is input into the initial azimuth calculation model to generate the compensated azimuth calculation model; the compensated azimuth calculation 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.

[0030] The compensated azimuth angle calculation model is iteratively optimized using the gradient descent algorithm until the compensated azimuth angle calculation model converges. Then, the azimuth angle parameters are extracted to generate real-time azimuth angle parameters for cloud seeding operations.

[0031] In one possible implementation of the first aspect, the step of iteratively optimizing the compensated azimuth angle calculation model using a gradient descent algorithm until the compensated azimuth angle calculation model converges, and then extracting the azimuth angle parameters to generate real-time azimuth angle parameters for cloud seeding operations, includes:

[0032] Obtain the azimuth calculation error and error gradient direction for the current iteration period;

[0033] The product of the gradient direction and the learning rate is calculated based on the error gradient direction, and the compensation weight coefficients are incrementally updated to adjust the compensation weight coefficients in the real-time azimuth compensation parameter set.

[0034] The outputs of the displacement compensation module and curvature correction module are recalculated based on the adjusted compensation weight coefficients, and the internal state variables of the model are updated synchronously to update the compensated azimuth angle solution model.

[0035] Repeat the above operations until the compensated azimuth angle calculation model meets the convergence condition. Extract the azimuth angle parameters from the converged azimuth angle calculation model and generate the real-time azimuth angle parameters for the cloud seeding operation. The convergence condition includes the error change rate of the azimuth angle calculation error being less than a set threshold or reaching the maximum number of iterations.

[0036] In one possible implementation of the first aspect, obtaining the azimuth calculation error and error gradient direction of the current iteration period includes:

[0037] Extract the measured azimuth angle data of the weather modification operation within the current iteration cycle. The measured azimuth angle data of the weather modification operation is collected by the angle measuring equipment of the ground reference station and synchronized with the Beidou positioning data in time.

[0038] The compensated azimuth angle calculation model is invoked to generate azimuth angle prediction data for the current iteration period. The azimuth angle prediction data includes the azimuth angle sequence and its confidence interval output by the compensated azimuth angle calculation model.

[0039] The measured azimuth data and the predicted azimuth data are time-aligned, and the absolute error and relative error are calculated point by point to obtain the azimuth calculation error.

[0040] Based on the azimuth angle calculation error, the compensation weight coefficients in the real-time azimuth angle compensation parameter set are backpropagated to generate the error gradient direction; wherein, the backpropagation process includes constructing a computation graph and taking derivatives layer by layer along the topology of the compensated azimuth angle calculation model to obtain the contribution gradient of each compensation weight coefficient to the total error.

[0041] In one possible implementation of the first aspect, the step of calculating the product of the gradient direction and the learning rate based on the error gradient direction, and incrementally updating the compensation weight coefficients to adjust the compensation weight coefficients in the real-time azimuth compensation parameter set, includes:

[0042] The error gradient direction is normalized to generate a gradient adjustment step size. The normalization process includes calculating the L2 norm of the gradient vector and dividing each gradient component by the norm value.

[0043] The first compensation weight coefficient and the second compensation weight coefficient are updated according to the gradient adjustment step size.

[0044] The updated first and second compensation weight coefficients are rewritten into the dynamic compensation weight set. During the update, spatial regularization constraints are introduced to ensure that the difference in compensation weight coefficients between adjacent positioning nodes does not exceed a set threshold, and the weight records of the most recent N iterations are retained. The value of N is inversely proportional to the spatial distribution density of the positioning nodes.

[0045] In one possible implementation of the first aspect, the step of extracting the azimuth parameters from the converged azimuth calculation model and generating real-time azimuth parameters for cloud seeding operations includes:

[0046] Obtain the parameter space distribution of the converged azimuth angle solution model, wherein the parameter space distribution includes the weight matrix of the displacement compensation module and the bias vector of the curvature correction module;

[0047] Calculate the eigenvectors of the parameter covariance matrix of the parameter spatial distribution, and select the top K principal components according to the size of the eigenvector values ​​to extract the azimuth principal component features;

[0048] The kernel density estimation method is used to perform non-parametric probability modeling on the principal component features to construct the probability density function of the azimuth parameter;

[0049] The final real-time azimuth angle parameters for the shadowing operation are determined based on the extreme points of the probability density function. These extreme points are obtained by traversing the domain of the probability density function and finding local maximum points.

[0050] In one possible implementation of the first aspect, the step of using a kernel density estimation method to perform non-parametric probabilistic modeling of the principal component features and constructing a probability density function for the azimuth parameter includes:

[0051] The azimuth principal component features are subjected to kernel density estimation to generate an initial probability density distribution; the kernel density estimation uses a Gaussian kernel function, and the bandwidth parameter is adaptively adjusted according to the variance of the principal component features.

[0052] The initial probability density distribution is calibrated based on the three-dimensional coordinate sequence in the BeiDou satellite positioning data set to generate a calibrated probability density function.

[0053] The peak region of the calibrated probability density function is extracted, and the median point of the peak region is used as the real-time azimuth angle parameter for the cloud seeding operation. The peak region is obtained by binarizing the probability density function by setting a density threshold, and the median point is calculated using the region centroid algorithm.

[0054] In another aspect, embodiments of the present invention also provide a real-time calculation system for the azimuth angle of human shadowing operations based on BeiDou satellites, including a processor and a machine-readable storage medium. The machine-readable storage medium is connected to the processor. The machine-readable storage medium is used to store programs, instructions, or code. The processor is used to execute the programs, instructions, or code in the machine-readable storage medium to implement the above-mentioned method.

[0055] Based on the above, this embodiment of the invention acquires BeiDou satellite positioning data of the target area and performs dynamic trajectory fitting processing on it. This allows for the accurate extraction of dynamic positioning features, including instantaneous displacement deviation and continuous trajectory curvature changes. Furthermore, by analyzing the environmental interference features of these dynamic positioning features, environmental interference factors affecting azimuth calculation can be accurately identified. Based on the coupling relationship between dynamic positioning features and environmental interference features, real-time azimuth compensation parameters are determined, effectively eliminating the impact of environmental interference on azimuth calculation. Finally, through iterative azimuth calculation, real-time azimuth parameters for weather modification operations in the target area are generated, significantly improving the real-time performance and accuracy of azimuth calculation. This significantly enhances the implementation effect and safety of weather modification operations, providing data support for weather modification operations in the meteorological field. Attached Figure Description

[0056] Figure 1It is a schematic execution flow diagram of the real-time calculation method for the azimuth angle of cloud seeding operations based on Beidou satellites provided by an embodiment of the present invention.

[0057] Figure 2 It is a schematic diagram of exemplary hardware and software components of a real-time calculation system for the azimuth angle of cloud seeding operations based on Beidou satellites provided by an embodiment of the present invention. Detailed implementation manners

[0058] The present invention will be specifically described below with reference to the accompanying drawings of the specification. Figure 1 It is a schematic flow diagram of a real-time calculation method for the azimuth angle of cloud seeding operations based on Beidou satellites provided by an embodiment of the present invention. The real-time calculation method for the azimuth angle of cloud seeding operations based on Beidou satellites will be introduced in detail below.

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

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

[0061] These positioning nodes will continuously record their position information within a continuous time period. Assume that the time period is represented by T, from time t1 to time t2 (t1 < t2). Each positioning node will generate 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 three-dimensional coordinate sequences of all these positioning nodes are combined together to form a set of Beidou satellite positioning data.

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

[0063] After obtaining the BeiDou satellite positioning data set, the positioning nodes will be affected by various factors in the actual environment, resulting in position changes. In order to accurately analyze their motion trajectory and characteristics, it is necessary to perform dynamic trajectory fitting processing on the data set.

[0064] Step S121: Extract the spatiotemporal distribution information of the positioning nodes in the three-dimensional coordinate sequence and generate spatiotemporal correlation features.

[0065] For a three-dimensional coordinate sequence, each coordinate point contains both temporal and spatial information. The temporal information indicates the exact time the coordinate point was recorded, while the spatial information is the three-dimensional coordinates (x, y, z). First, the three-dimensional coordinate sequence of each positioning node is arranged 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 time the coordinate point was recorded.

[0066] Next, analyze the spatial relationship between the positioning nodes at different times. This relationship can be represented by calculating the distance and time interval between adjacent coordinate points. 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) and (xik+1, yik+1, zik+1) at adjacent times tk and tk+1 (the spatial distance dik is the square root of the sum of the squares of the differences between the corresponding coordinate points at adjacent times in the x, y, and z directions), and the time interval Δtk = tk+1-tk.

[0067] By combining the distance and time interval information of all the positioning nodes, the spatiotemporal distribution information of the positioning nodes can be obtained. Furthermore, to generate spatiotemporal correlation features, this spatiotemporal distribution information can be encoded. For example, it can be represented as a vector, where a series of distance and time interval values ​​for each positioning node are arranged sequentially. Assuming positioning node i has n coordinate points, this would 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 from all m positioning nodes constitutes the spatiotemporal correlation features.

[0068] Step S122: Perform adaptive sliding window segmentation 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 time span of the sliding window. Kalman filtering is used to smooth the noise of the discrete coordinate points within the sliding window.

[0069] After obtaining the spatiotemporal correlation features, adaptive sliding window segmentation is required to further analyze the trajectory features of the positioning nodes. The sliding window is a fixed-size time period that contains the coordinate information of multiple positioning nodes.

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

[0071] If the spatial distribution density is low, it indicates that the distance between the positioning nodes is relatively large, and their movement changes may be more drastic. In this case, a smaller sliding window time span is needed to capture their trajectory changes more precisely. For example, when the spatial distribution density is below a certain threshold ρ0, the sliding window time span is set to Δt1; when the spatial distribution density is above ρ0, the sliding window time span is set to Δt2, and Δt1 < Δt2.

[0072] Then, a sliding window is used to slide across the spatiotemporal correlation features. Assuming the starting time of the sliding window is tstart and the ending time is tend = tstart + Δt (where Δt is the time span of the sliding window), this time period contains the coordinate information of multiple positioning nodes, which constitute a windowed trajectory segment. As the sliding window continues to slide, multiple windowed trajectory segments can be obtained.

[0073] For each discrete coordinate point within a windowed trajectory segment, noise may exist due to factors such as external environmental interference. To eliminate this noise, Kalman filtering is employed. Kalman filtering is a recursive optimal estimation algorithm that estimates the system state at the current moment based on the system's state equation and observation equation, combined with the estimated value from the previous moment and the observed value from the current moment. In this scenario, the system state can be considered as the actual position of the positioning node, and the observed values ​​are the discrete coordinate points within the windowed trajectory segment. Through Kalman filtering, these discrete coordinate points can be smoothed to obtain more accurate trajectory information.

[0074] Step S123: Use 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.

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

[0076] The adaptive piecewise interpolation method divides the windowed trajectory segment into multiple smaller segments based on the distribution of discrete coordinate points. For each segment's discrete coordinate points, an appropriate interpolation function is used for interpolation. For example, a cubic spline interpolation function can be used, which ensures that the function value, first derivative, and second derivative are continuous at the interpolation points, thus obtaining a smoother curve.

[0077] When segmenting, adaptive adjustments should be made based on the curvature changes of the discrete coordinate points. If the curvature change of a certain segment is large, it indicates that the trajectory changes 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.

[0078] For each discrete coordinate point (x1, y1, z1), (x2, y2, z2), ..., (xk, yk, zk) within a small segment, a cubic spline interpolation function is used to interpolate the x, y, and z coordinates respectively. Taking the x-coordinate as an example, the coefficients of the cubic spline interpolation function are determined by solving a series of linear equations, thus obtaining the interpolation curve for the x-coordinate within that small segment. Similarly, the same process is applied to the y and z coordinates. By combining the interpolation curves of all small segments, a smooth trajectory curve can be generated.

[0079] Step S124: Calculate the displacement deviation between adjacent positioning nodes based on the smooth trajectory curve, and generate the instantaneous displacement deviation.

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

[0081] For each time t on the smooth trajectory curve, there is a corresponding position for each positioning node. 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 between adjacent positioning nodes i and i+1 at time t is, i+1(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 between adjacent positioning nodes in the x, y, and z directions).

[0082] By combining the displacement deviations at all times, we obtain the instantaneous displacement deviation between adjacent positioning nodes. This calculation is performed for all positioning node pairs (i, i+1) (i=1, 2, ..., m-1) to finally obtain the instantaneous displacement deviation of the entire set of positioning nodes.

[0083] Step S125: Calculate the curvature derivative of the smooth trajectory curve at adjacent interpolation points, determine the direction of curvature change based on the change in the sign of the curvature derivative, and extract the curvature change amount of the continuous trajectory.

[0084] To further analyze the curvature changes 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 based on its sign change, and then extract the curvature change of the continuous trajectory.

[0085] First, the curvature of the smooth trajectory curve at each interpolation point is calculated. Curvature is a quantity describing the degree of bending of a curve, and it can be calculated using the first and second derivatives of the curve. For the smooth trajectory curve at the interpolation point (x(t), y(t), z(t)), the calculation of its curvature κ(t) involves taking the first and second derivatives of x(t), y(t), and z(t) respectively, and then calculating it using conventional relevant technical formulas (the specific formulas can be obtained from the definition of curvature of space curves).

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

[0087] The direction of curvature change is determined by the sign of the derivative of curvature. If the derivative of curvature is greater than 0, it indicates that the curvature is increasing, meaning the curve is becoming more curved; if the derivative of curvature is less than 0, it indicates that the curvature is decreasing, meaning the curve is becoming less curved.

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

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

[0090] After obtaining the instantaneous displacement deviation and the continuous trajectory curvature change, they are merged into a dynamic positioning feature set. The instantaneous displacement deviation and the continuous trajectory curvature change can be considered 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 time points), then the dynamic positioning feature set F can be represented as F=[D, C]. This means that the two vectors are connected sequentially to form a new vector, which contains information about the instantaneous displacement deviation and continuous trajectory curvature change of the positioning nodes, providing a foundation for subsequent environmental interference feature analysis.

[0091] Step S130: Based on the dynamic positioning feature set, perform environmental interference feature analysis on the target area to generate an environmental interference feature set.

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

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

[0094] To analyze the characteristics of environmental disturbances, the first step is to obtain a set of meteorological monitoring data for the target area. Multiple meteorological monitoring devices will be deployed within the target area to monitor meteorological data in real time. The meteorological monitoring data set mainly includes wind speed gradient distribution data and atmospheric refractive index distribution data.

[0095] Wind speed gradient distribution data describes the variation of wind speed at different locations within a target area. Meteorological monitoring equipment measures wind speed at different altitudes and horizontal positions, obtaining a series of wind speed data. Organizing this wind speed data according to location information yields wind speed gradient distribution data. For example, measuring wind speed at different altitudes h1, h2, ..., hk, and at different horizontal positions (x1, y1), (x2, y2), ..., (xl, yl), yields wind speed values ​​vij (where i represents the altitude number and j represents the horizontal position number). Combining these wind speed values ​​with the corresponding location information constitutes the wind speed gradient distribution data.

[0096] 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, thus interfering with the transmission of BeiDou satellite signals. Meteorological monitoring equipment calculates atmospheric refractive index by measuring parameters such as atmospheric temperature, humidity, and pressure, using conventional formulas. Similarly, by organizing atmospheric refractive index data from different locations according to location information, atmospheric refractive index distribution data is obtained.

[0097] Step S132: Perform vector decomposition on the wind speed gradient distribution data to obtain the lateral wind speed component and the longitudinal wind speed component. After standardizing the instantaneous displacement deviation, establish a nonlinear regression model between the instantaneous displacement deviation and the wind speed component. Construct a feature interaction matrix by introducing the quadratic term and cross term of the wind speed component. Train the regression model using the random forest algorithm and calculate the feature importance score. Combine standardized residual analysis to determine the dynamic contribution weight of each wind speed component. Compensate and correct the nonlinear residual part using Gaussian process regression to generate wind speed interference features.

[0098] After obtaining the wind speed gradient distribution data, it is first decomposed into vectors. Wind speed is a vector with magnitude and direction. Assuming the wind speed vector is V, its projection onto the horizontal plane can be decomposed into a lateral wind speed component Vx and a longitudinal wind speed component Vy. Based on the wind speed direction angle θ, the lateral 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.

[0099] To eliminate the influence of different orders of magnitude, the instantaneous displacement deviation in the dynamic positioning feature set needs to be standardized. Common methods for standardization include Z-score standardization. Assuming the instantaneous displacement deviation is D, with a mean of μD and a standard deviation of σD, then the standardized instantaneous displacement deviation D' = (D - μD) / σD.

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

[0101] The random forest algorithm was used to train the regression model. Random forest is an ensemble learning algorithm consisting of multiple decision trees. During training, random forest randomly samples features and samples data to construct multiple decision trees, and then combines the results of these decision trees. By training the random forest regression model, the importance score of each feature can be calculated, that is, the degree of influence of each wind speed component and its quadratic and cross terms on the instantaneous displacement deviation.

[0102] The dynamic contribution weights of each wind speed component are determined by combining standardized residual analysis. Standardized residuals are the standardized values ​​obtained by standardizing the difference between the predicted and actual values ​​from the regression model. Analyzing the standardized residuals reveals nonlinear components in the regression model. Gaussian process regression is used to compensate for and correct these nonlinear residuals. Gaussian process regression is a nonparametric regression method that can model any function. Modeling and predicting nonlinear residuals using Gaussian process regression yields more accurate prediction results.

[0103] By combining the dynamic contribution weights of each wind speed component with the compensated and corrected prediction results, wind speed interference characteristics can be generated.

[0104] Step S133: Perform correlation analysis on the atmospheric refractive index distribution data and the continuous trajectory curvature change in the dynamic positioning feature set, calculate the covariance matrix of the atmospheric refractive index spatial gradient and the trajectory curvature change, and extract the principal component features as refractive index interference factors to generate refractive index interference features.

[0105] To analyze the impact 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 curvature change of continuous trajectories in the dynamic positioning feature set.

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

[0107] Then, the covariance matrix of the spatial gradient of atmospheric refractive index and the change in trajectory curvature is calculated. The covariance matrix measures the correlation between two variables. Assume the spatial gradient vector of atmospheric refractive index is... If the vector of trajectory curvature change 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 change.

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

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

[0110] Step S134: Construct an environmental interference feature set based on the wind speed interference features and the refractive index interference features; 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 of different wind speed gradients on the displacement deviation of the positioning node; the refractive index interference factor set is used to characterize the nonlinear effect of the spatial distribution of atmospheric refractive index on the change of trajectory curvature.

[0111] After obtaining the wind speed interference features and refractive index interference features, they are combined to construct an environmental interference feature set. The wind speed interference features include the wind speed contribution weight, while the refractive index interference features include 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 feature set.

[0112] Specifically, for each positioning node, a set of wind speed contribution weights is determined from the wind speed interference characteristics. These weights are calculated through previous analysis of wind speed gradient distribution data and instantaneous displacement deviations. They can quantify the degree of influence 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 different degrees of pushing or hindering effects on the displacement of the positioning node, and the wind speed contribution weights reflect the specific quantification of this influence.

[0113] Simultaneously, for each positioning node, its corresponding refractive index interference factor is extracted from the refractive index interference characteristics. This factor is obtained through correlation analysis of atmospheric refractive index distribution data and continuous trajectory curvature changes, followed by steps such as covariance matrix calculation and principal component analysis. It is used to characterize the nonlinear effect of the spatial distribution of atmospheric refractive index on the trajectory curvature change of the positioning node. For example, changes in atmospheric refractive index at different spatial locations may cause bending of the BeiDou satellite signal propagation path, thereby affecting the trajectory curvature of the positioning node; the refractive index interference factor reflects the strength of this effect.

[0114] By arranging and combining the wind speed contribution weights and refractive index interference factors of all positioning nodes in the order of the positioning nodes, an environmental interference feature set is constructed. This environmental interference feature set is crucial for subsequent analysis of the impact of environmental factors on the azimuth angle of weather modification operations, because it integrates the influence information of two important environmental factors, wind speed and atmospheric refractive index, on the motion characteristics of positioning nodes.

[0115] Step S140: Determine the real-time azimuth compensation parameter set based on the coupling relationship between the dynamic positioning feature set and the environmental interference feature set.

[0116] 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 thus affect the calculation of the azimuth angle of the cloud seeding operation, it is necessary to determine the real-time azimuth angle compensation parameter set based on this coupling relationship in order to correct the subsequent azimuth angle calculation.

[0117] Step S141: Perform weight allocation 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 the first compensation weight coefficient.

[0118] When determining the set of real-time azimuth compensation parameters, 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 positional change of the positioning node in a short period of time, while the wind speed contribution weight reflects the degree of influence of wind speed on the displacement deviation of the positioning node.

[0119] To combine these two factors, weight allocation is necessary. A method based on their correlation can be used. For example, first calculate the correlation coefficient between the instantaneous displacement deviation and the wind speed contribution weight; this coefficient measures the degree of their association. Assuming the correlation coefficient is r, a function f(r) related to r is used to determine the weight allocation ratio. Let w1 be the weight of the instantaneous displacement deviation and w2 be the weight of the wind speed contribution weight, with w1 + w2 = 1. The values ​​of w1 and w2 can be adjusted based on the value of f(r). When r is large, it indicates a strong correlation, and a relatively balanced weight may be needed; when r is small, the weight allocation may need to be adjusted according to specific circumstances.

[0120] After this weighting process, the instantaneous displacement deviation and wind speed contribution weights are combined according to their respective weights to generate the first compensation weight coefficient. This first compensation weight coefficient comprehensively considers the impact of instantaneous displacement deviation and wind speed on the position change of the positioning node.

[0121] Step S142: Perform weighted allocation 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.

[0122] Similarly, weighted processing is applied to the continuous trajectory curvature change in the dynamic positioning feature set and the refractive index interference factor in the environmental interference feature set. The continuous trajectory curvature change describes the curvature variation of the positioning node trajectory, while the refractive index interference factor reflects the intensity of the influence of atmospheric refractive index on the trajectory curvature change.

[0123] First, calculate the correlation coefficient between the change in curvature of the continuous trajectory and the refractive index interference factor, denoted as r'. Then, determine the weight allocation ratio using a function g(r') related to r'. Let w3 be the weight of the change in curvature of the continuous trajectory, and w4 be the weight of the refractive index interference factor, with w3 + w4 = 1. Adjust the values ​​of w3 and w4 according to the value of g(r').

[0124] The continuous trajectory curvature change and 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 influence of trajectory curvature change and atmospheric refractive index on the positioning node trajectory, and is another important parameter for real-time azimuth compensation.

[0125] Step S143: Construct a dynamic compensation weight set based on the first compensation weight coefficient and the second compensation weight coefficient; the dynamic compensation weight set includes the displacement compensation weight and curvature compensation weight corresponding to each positioning node.

[0126] 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 displacement compensation weight of that node, as it is mainly based on the instantaneous displacement deviation and wind speed contribution weight, reflecting the compensation information for the displacement of the positioning node; the second compensation weight coefficient corresponds to the curvature compensation weight of that node, as it is based on the continuous trajectory curvature change and refractive index interference factor, reflecting the compensation information for the trajectory curvature of the positioning node.

[0127] By arranging and combining the displacement compensation weights and curvature compensation weights of all positioning nodes in the order of the positioning nodes, a dynamic compensation weight set is formed. This dynamic compensation weight set provides the foundation for accurate compensation of the displacement and trajectory curvature of the positioning nodes in subsequent azimuth calculations.

[0128] Step S144: After standardizing the instantaneous displacement deviation and the continuous trajectory curvature change, based on the displacement compensation weight and curvature compensation weight in the dynamic compensation weight set, the fusion weight coefficient of each standardized result is calculated by the entropy weight method. The weighted geometric average algorithm is used to perform weighted fusion processing on each standardized result to generate a real-time azimuth compensation parameter set. The real-time azimuth compensation parameter set includes a subset of displacement compensation parameters and a subset of curvature compensation parameters. Each subset is arranged in time sequence and corresponds one-to-one with the positioning node.

[0129] After constructing the dynamic compensation weight set, in order to generate the real-time azimuth compensation parameter set, the instantaneous displacement deviation and the continuous trajectory curvature change in the dynamic positioning feature set are first standardized. Standardization can eliminate the influence of different orders of magnitude, making 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.

[0130] Then, based on the displacement compensation weights and curvature compensation weights in the dynamic compensation weight set, the fusion weight coefficients of each standardized result are calculated using the entropy weight method. The entropy weight method is an objective weighting method that determines weights based on the dispersion of the data. 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 disorder in the data. Their entropy weights are calculated based on the information entropy, denoted as wD and wC, where wD + wC = 1.

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

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

[0133] Step S150: Based on the real-time azimuth compensation parameter set, perform azimuth iterative calculation on the BeiDou satellite positioning data set to generate real-time azimuth parameters for cloud seeding operations in the target area.

[0134] After obtaining the set of real-time azimuth compensation parameters, in order to obtain the real-time azimuth parameters of the shadowing operation in the target area, it is necessary to perform azimuth iterative calculation processing on the BeiDou satellite positioning data set based on the set of compensation parameters.

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

[0136] First, a three-dimensional coordinate sequence is retrieved from the BeiDou satellite positioning dataset to construct an initial azimuth angle calculation model. Least squares 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 calculation model, the coordinate points in the three-dimensional coordinate sequence are used as observation data, assuming a functional relationship between the azimuth angle and these coordinate points. The parameters in this functional relationship are estimated using the least squares method to minimize the sum of squared errors between the observed data and the model's predicted values.

[0137] Specifically, let the coordinates in the three-dimensional coordinate sequence be (x1, y1, z1), (x2, y2, z2), ..., (xn, yn, zn). Assume the relationship between the azimuth angle α and these coordinates can be represented by a function f(x, y, z, α). Using the least squares method, find the α value that minimizes Σ[f(xi, yi, zi, α) - yi_obs]², where yi_obs is the actually observed azimuth-related value. The obtained α value is the azimuth reference value. This completes the construction of the initial azimuth angle calculation model.

[0138] Step S152: Input the real-time azimuth compensation parameter set into the initial azimuth calculation model to generate the compensated azimuth calculation model; the compensated azimuth calculation 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.

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

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

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

[0142] Through the functions of these two modules, the compensated azimuth calculation model can more comprehensively consider the impact of environmental factors on the movement of positioning nodes, thereby calculating the azimuth of the shadowing operation more accurately.

[0143] Step S153: The compensated azimuth angle calculation model is iteratively optimized using the gradient descent algorithm until the compensated azimuth angle calculation model converges. Then, the azimuth angle parameters are extracted to generate real-time azimuth angle parameters for cloud seeding operations.

[0144] After generating the compensated azimuth calculation model, in order to obtain more accurate real-time azimuth parameters for cloud seeding operations, the model needs to be iteratively optimized using the gradient descent algorithm.

[0145] Step S1531: Obtain the azimuth calculation error and error gradient direction for the current iteration cycle.

[0146] In each iteration cycle, the azimuth angle calculation error and error gradient direction must first be obtained. To obtain the azimuth angle calculation error, the measured azimuth angle data of the cloud seeding operation within the current iteration cycle needs to be extracted. This measured data is collected by the angle measuring equipment of the ground reference station, and it is necessary to ensure that this measured data is synchronized with the BeiDou positioning data in time. Only in this way can the model calculation results and the actual measurement results be accurately compared.

[0147] Simultaneously, the compensated azimuth angle calculation model is invoked to generate azimuth angle prediction data for the current iteration cycle. This prediction data includes the azimuth angle sequence output by the compensated azimuth angle calculation model and its confidence interval. The measured azimuth angle data and the predicted azimuth angle data are time-aligned, and then the absolute and relative errors are calculated point-by-point. The azimuth angle calculation error is obtained by comprehensively calculating these errors.

[0148] Based on the azimuth calculation error, backpropagation processing is performed on the compensation weight coefficients in the real-time azimuth compensation parameter set. Backpropagation involves constructing a computational graph and calculating derivatives layer by layer along the topology of the compensated azimuth calculation model. This yields the gradient of each compensation weight coefficient's contribution to the total error; the direction of this gradient is the error gradient direction. The error gradient direction indicates in the parameter space which direction adjusting the compensation weight coefficients will reduce the azimuth calculation error.

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

[0150] After obtaining the error gradient direction, in order to adjust the compensation weight coefficients 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. The normalization process includes 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.

[0151] 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 and second compensation weight coefficients are updated based on this 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.

[0152] The updated first and second compensation weight coefficients are rewritten into the dynamic compensation weight set. Spatial regularization constraints are introduced during the update to ensure that the difference in compensation weight coefficients between adjacent positioning nodes does not exceed a set threshold, thus guaranteeing the spatial smoothness of the compensation weight coefficients. Simultaneously, the weight records from the most recent N iterations are retained; the value of N is 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 iteration of optimization.

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

[0154] 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 according to the adjusted displacement compensation weight coefficients and the subset of displacement compensation parameters; the curvature correction module re-corrects the trajectory curvature information of the positioning node according to the adjusted curvature compensation weight coefficients and the subset of curvature compensation parameters.

[0155] Simultaneously, the model's internal state variables are updated. These state variables record intermediate results and parameter information during the calculation process. By updating these state variables, the compensated azimuth angle calculation model can reflect the changes after the adjustment of the compensation weight coefficients, thereby updating the entire model.

[0156] Step S1534: Repeat the above operation until the compensated azimuth angle solution model meets the convergence condition, extract the azimuth angle parameters in the converged azimuth angle solution model, and generate the real-time azimuth angle parameters for the cloud seeding operation. The convergence condition includes the error change rate of the azimuth angle solution error being less than a set threshold or reaching the maximum number of iterations.

[0157] The steps S1531-S1533 are repeated continuously. In each iteration, the compensation weight coefficients are adjusted, and the compensated azimuth calculation model is updated, gradually reducing the azimuth calculation error. During the iteration process, it is determined whether the compensated azimuth calculation model meets the convergence conditions. Convergence conditions include: the rate of change of the azimuth calculation error being less than a set threshold (i.e., the error reduction rate becomes very slow, indicating the model is close to the optimal solution); or reaching the maximum number of iterations. To avoid infinite iteration, a maximum number of iterations is set, and iteration stops when this maximum number of iterations is reached.

[0158] When the compensated azimuth calculation model meets the convergence condition, the azimuth parameters in the converged azimuth calculation model are extracted. These azimuth parameters are obtained after multiple iterations of optimization and can more accurately reflect the real-time azimuth of the weather modification operation in the target area. Combining these azimuth parameters together generates the real-time azimuth parameters for the weather modification operation.

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

[0160] After obtaining the converged azimuth calculation model, it is necessary to acquire its parameter space distribution. The parameter space distribution reflects the values ​​of each parameter in the model. For the compensated azimuth calculation model, it mainly includes the weight matrix of the displacement compensation module and the bias vector of the curvature correction module.

[0161] The weight matrix of the displacement compensation module records the weight distribution of each parameter when processing a subset of displacement compensation parameters. The magnitude of the element values ​​in this weight matrix reflects the degree of influence of different displacement compensation parameters on 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.

[0162] The bias vector of the curvature correction module records the bias information used to adjust the model output when processing a subset of curvature compensation parameters. Each component of the bias vector corresponds to an output dimension of the curvature correction module. By adjusting the value of the bias vector, the output of the curvature correction module can be shifted and adjusted as a whole, thereby better fitting the actual azimuth angle.

[0163] Step S1534-12: Calculate the eigenvectors of the parameter covariance matrix of the parameter spatial distribution, and select the top K principal components according to the size of the eigenvector values ​​to extract the azimuth principal component features.

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

[0165] Next, the eigenvectors of the parametric covariance matrix are calculated. An eigenvector is a special type of vector within a matrix that undergoes only scaling, not rotation, during matrix transformations. Each eigenvector corresponds to an eigenvalue, the magnitude of which reflects the degree of data dispersion in the direction represented by that eigenvector.

[0166] The eigenvectors are sorted according to their values, and the principal components corresponding to the top K eigenvectors are selected. Principal component analysis (PCA) is a data dimensionality reduction method. By selecting the top K principal components, the high-dimensional parameter space can be reduced to K dimensions while retaining the main information of the data.

[0167] Projecting the parameters of the converged azimuth calculation model onto the first K principal components yields projection values ​​that constitute the azimuth principal component features. These principal component features are a dimensionality-reduced representation of the original parameter space. They retain the most important information in the parameter space, can more concisely describe the characteristics of the azimuth, and at the same time reduce the dimensionality of the data, thus reducing the complexity of subsequent processing.

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

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

[0170] Kernel density estimation is a non-parametric estimation method that does not require assumptions about the distribution of the data; instead, it estimates the probability density of the data through a kernel function. First, kernel density estimation is performed on the principal component features, using a Gaussian kernel function. The Gaussian kernel function possesses smoothness and symmetry, enabling it to fit the data distribution well. When using the Gaussian kernel function for kernel density estimation, the bandwidth parameter needs to be determined. The bandwidth parameter controls the width of the kernel function, and its magnitude affects the smoothness of the estimation result. Here, the bandwidth parameter is adaptively adjusted according to the variance of the principal component features. When the variance is large, the bandwidth parameter is appropriately increased to cover the data more broadly; when the variance is small, the bandwidth parameter is appropriately decreased to more finely characterize the local features of the data.

[0171] 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 using a three-dimensional coordinate sequence from the BeiDou satellite positioning dataset. The three-dimensional coordinate sequence contains the actual position information of the positioning nodes. 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, thus generating a calibrated probability density function.

[0172] Step S1534-14: Determine the final real-time azimuth angle parameters of the shadowing operation based on the extreme points of the probability density function. The extreme points are obtained by traversing the domain of the probability density function and finding the local maximum points.

[0173] After constructing the probability density function of the azimuth parameters, the final real-time azimuth parameters for cloud seeding operations are determined based on the extreme points of this function. Extreme points are the points where the function reaches its maximum or minimum value within its domain. For a probability density function, the extreme points typically represent locations where the azimuth angle has a high probability of occurring.

[0174] To find the extreme points of a probability density function, it is necessary to traverse its domain. The domain is the range of values ​​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 ​​at adjacent points are compared. When the function value at a point is found to be greater than the function values ​​at its neighboring points, that point is a local maximum. By continuously traversing and comparing, all local maximum points are found.

[0175] 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 weather modification operation. This selection process can be carried out according to 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 occurrence and is most likely to be the actual azimuth angle for the weather modification operation.

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

[0177] In kernel density estimation, a Gaussian kernel function is used to estimate the probability density of the principal component features of the azimuth angle. The Gaussian kernel function is based on a normal distribution, which allows it to fit the data smoothly. For each principal component feature data point, the Gaussian kernel function is used, centered on that point, to calculate its contribution to the probability density of the surrounding area.

[0178] When using the Gaussian kernel function, the choice of bandwidth parameter is crucial. The bandwidth parameter determines the width of the Gaussian kernel function, that is, the range of influence of the data points on the surrounding area. Here, the bandwidth parameter is adaptively adjusted based on the variance of the principal component features. Variance reflects the dispersion of the principal component feature data. If the variance is large, it indicates that the data is relatively dispersed, requiring a larger bandwidth parameter so that the Gaussian kernel function can cover a wider area, thus providing a more comprehensive estimate of the probability density. If the variance is small, it indicates that the data is relatively concentrated, requiring a smaller bandwidth parameter to more accurately characterize the local features of the data.

[0179] By applying a Gaussian kernel function to each principal component feature data point and summing the contributions of all data points, an initial probability density distribution is obtained. This initial probability density distribution describes the probability of the azimuth principal component feature under different values.

[0180] Step S1534-142: The initial probability density distribution is calibrated based on the three-dimensional coordinate sequence in the BeiDou satellite positioning data set to generate a calibrated probability density function.

[0181] The initial probability density distribution is obtained based on the principal component features of the azimuth angle, which may deviate from the actual azimuth angle. To make the probability density function more accurately reflect the actual situation, the initial probability density distribution needs to be calibrated based on the three-dimensional coordinate sequence in the BeiDou satellite positioning dataset.

[0182] The three-dimensional coordinate sequence contains the actual position information of the positioning nodes, which is closely related to the azimuth angle. First, the relationship between each coordinate point in the three-dimensional 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 and adjusted with these azimuth angle ranges.

[0183] For example, if a location node appears frequently in a certain area in a 3D coordinate sequence, it indicates that the actual probability of the corresponding azimuth angle appearing in that area is relatively 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 a location node appears frequently in a certain area, the probability value within the corresponding azimuth angle range needs to be appropriately decreased.

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

[0185] Steps S1534-143: Extract the peak region of the calibrated probability density function, and use the median point of the peak region as the real-time azimuth angle parameter for the cloud seeding operation; wherein, the peak region is obtained by binarizing the probability density function by setting a density threshold, and the median point is calculated using the region centroid algorithm.

[0186] After obtaining the calibrated probability density function, in order to determine the final real-time azimuth parameters for the weather modification operation, it is necessary to extract the peak region of the calibrated probability density function. The peak region represents the area where the azimuth angle has a high probability of occurrence, and usually contains the most likely azimuth angle for the weather modification operation.

[0187] Peak regions are obtained by binarizing the probability density function by setting a density threshold. The density threshold is a pre-defined 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 binarization segmentation, the calibrated probability density function is divided into peak and non-peak regions.

[0188] For the extracted peak regions, the median point of the peak region is calculated and used as the real-time azimuth parameter for weather modification operations. The median point is calculated using the regional centroid algorithm. The regional centroid algorithm is a method for calculating the center position of a region, considering the position and weight of each point within the region. In the peak region, the weight of each point can be determined based on its corresponding probability density value; the higher the probability density value, the greater the weight. The result of calculating a weighted average of all points within the peak region is the centroid of the peak region, i.e., the median point. This median point represents the center position of the peak region, and using it as the real-time azimuth parameter for weather modification operations can more accurately reflect the actual azimuth situation of the operation.

[0189] Throughout the process, data collection, including BeiDou satellite positioning data and meteorological monitoring data, strictly adhered to relevant laws and regulations to ensure the legality and compliance of the data. For potentially privacy-sensitive data, data encryption technology was employed for privacy protection and leak prevention. During data transmission, secure encryption protocols were used to encrypt the data, preventing it from being stolen or tampered with during transmission. For data storage, a secure storage system was used to encrypt and store the data, with strict access control measures in place, ensuring that only authorized personnel can access and process the data.

[0190] Figure 2 The illustration shows exemplary hardware and software components of a BeiDou satellite-based real-time azimuth angle calculation system 100 for weather modification operations, which can implement the ideas of this application, according to some embodiments of this application. For example, a processor 120 can be used in the BeiDou satellite-based real-time azimuth angle calculation system 100 for weather modification operations and to perform the functions in this application.

[0191] The BeiDou satellite-based real-time azimuth angle calculation system 100 for weather modification operations can be a general-purpose server or a special-purpose server; both can be used to implement the BeiDou satellite-based real-time azimuth angle calculation method for weather modification operations described in this application. 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 load.

[0192] For example, the BeiDou satellite-based real-time calculation system azimuth angle for weather modification operations 100 may 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. Exemplarily, the BeiDou satellite-based real-time calculation system azimuth angle for weather modification operations 100 may also include program instructions stored in ROM, RAM, or other types of non-transitory storage media, or any combination thereof. The methods of this application can be implemented according to these program instructions. The BeiDou satellite-based real-time calculation system azimuth angle for weather modification operations 100 also includes an I / O interface 150 between the computer and other input / output devices.

[0193] For ease of explanation, only one processor is described in the BeiDou satellite-based real-time calculation system for azimuth angle of weather modification operations. However, it should be noted that the BeiDou satellite-based real-time calculation system for azimuth angle of weather modification operations 100 may also include multiple processors. Therefore, the steps executed by one processor as described in this application may also be executed jointly by multiple processors or individually. For example, if the processor of the BeiDou satellite-based real-time calculation system for azimuth angle of weather modification operations 100 executes steps A and B, it should be understood that steps A and B may also be executed jointly by two different processors or individually by one processor. For example, the first processor executes step A, the second processor executes step B, or the first processor and the second processor jointly execute steps A and B.

[0194] Furthermore, this embodiment of the invention also provides a readable storage medium, which contains pre-set computer-executable instructions. When the processor executes the computer-executable instructions, the above-mentioned method for real-time calculation of azimuth angle for artificial weather modification operations based on BeiDou satellites is implemented.

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

Claims

1. A method for real-time calculation of azimuth angle for weather modification operations based on BeiDou satellites, characterized in that, The method includes: Acquire a set of BeiDou satellite positioning data for the target area, wherein the set of BeiDou satellite positioning data contains a sequence of three-dimensional coordinates of multiple positioning nodes within a continuous time period; The BeiDou satellite positioning data set is subjected to dynamic trajectory fitting processing to obtain a dynamic positioning feature set that includes instantaneous displacement deviation and continuous trajectory curvature change. Based on the dynamic positioning feature set, the target area is subjected to environmental interference feature analysis processing to generate an environmental interference feature set, 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. Based on the coupling relationship between the dynamic positioning feature set and the environmental interference feature set, the real-time azimuth compensation parameter set is determined; Based on the real-time azimuth compensation parameter set, the azimuth iterative calculation of the BeiDou satellite positioning data set is performed to generate real-time azimuth parameters for cloud seeding operations in the target area.

2. The method for real-time calculation of azimuth angle for weather modification operations based on BeiDou satellite as described in claim 1, characterized in that, The dynamic trajectory fitting process performed on the BeiDou satellite positioning data set yields a dynamic positioning feature set containing instantaneous displacement deviation and continuous trajectory curvature change, including: Extract the spatiotemporal distribution information of the positioning nodes in the three-dimensional coordinate sequence to generate spatiotemporal correlation features; The spatiotemporal correlation features are subjected to adaptive sliding window segmentation 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 time span of the sliding window. Kalman filtering is used to smooth the noise of discrete coordinate points within the sliding window. An adaptive piecewise interpolation method is used to perform nonlinear reconstruction processing on the discrete coordinate points in the windowed trajectory segment to generate a smooth trajectory curve. The displacement deviation between adjacent positioning nodes is calculated based on the smooth trajectory curve, and the instantaneous displacement deviation is generated. Calculate the curvature derivative of the smooth trajectory curve at adjacent interpolation points, determine the direction of curvature change based on the change in the sign of the curvature derivative, and extract 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 for shadowing operations based on BeiDou satellite as described in claim 1, characterized in that, The step of performing environmental interference feature analysis on the target area based on the dynamic positioning feature set to generate an environmental interference feature set includes: Acquire a set of meteorological monitoring data for the target area, wherein the set of meteorological monitoring data includes wind speed gradient distribution data and atmospheric refractive index distribution data; The wind speed gradient distribution data is vector decomposed to obtain the lateral wind speed component and the longitudinal wind speed component. After standardizing the instantaneous displacement deviation, a nonlinear regression model of instantaneous displacement deviation and wind speed component is established. By introducing the quadratic term and cross term of the wind speed component, a feature interaction matrix is ​​constructed. The regression model is trained using the random forest algorithm and the feature importance score is calculated. Combined with standardized residual analysis, the dynamic contribution weight of each wind speed component is determined. The nonlinear residual part is compensated and corrected using Gaussian process regression to generate wind speed interference features. The atmospheric refractive index distribution data and the continuous trajectory curvature change in the dynamic positioning feature set are correlated and analyzed. The covariance matrix of the atmospheric refractive index spatial gradient and the trajectory curvature change is calculated, and the principal component features are extracted as refractive index interference factors to generate refractive index interference features. Construct an environmental interference feature set based on the wind speed interference features and the refractive index interference features; The wind speed contribution weight set is used to quantify the degree of influence of different wind speed gradients on the displacement deviation of the positioning node; 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 trajectory curvature change.

4. The method for real-time calculation of azimuth angle for weather modification operations based on BeiDou satellite as described in claim 1, characterized in that, The step of determining the real-time azimuth compensation parameter set based on the coupling relationship between the dynamic positioning feature set and the environmental interference feature set includes: The instantaneous displacement deviation in the dynamic positioning feature set and the wind speed contribution weight in the environmental disturbance feature set are weighted and processed to generate a first compensation weight coefficient. The continuous trajectory curvature change in the dynamic positioning feature set and the refractive index interference factor in the environmental interference feature set are weighted and processed to generate a second compensation weight coefficient. A dynamic compensation weight set is constructed based on the first compensation weight coefficient and the second compensation weight coefficient; the dynamic compensation weight set includes the displacement compensation weight and curvature compensation weight corresponding to each positioning node; After standardizing the instantaneous displacement deviation and the continuous trajectory curvature change, the fusion weight coefficient of each standardized result is calculated using the entropy weight method based on the displacement compensation weight and curvature compensation weight in the dynamic compensation weight set. The weighted geometric average algorithm is then used to perform weighted fusion processing on each standardized result to generate a real-time azimuth compensation parameter set. The real-time azimuth compensation parameter set includes a subset of displacement compensation parameters and a subset of curvature compensation parameters. Each subset is arranged in a time sequence and corresponds one-to-one with the positioning node.

5. The method for real-time calculation of azimuth angle for shadowing operations based on BeiDou satellite as described in claim 4, characterized in that, The step of performing azimuth iterative calculation on the BeiDou satellite positioning data set based on the real-time azimuth compensation parameter set to generate real-time azimuth parameters for cloud seeding operations in the target area includes: An initial azimuth angle calculation model is constructed by calling the three-dimensional coordinate sequence in the BeiDou satellite positioning data set; the initial azimuth angle calculation model uses the least squares method to calculate the azimuth angle reference value. The real-time azimuth compensation parameter set is input into the initial azimuth calculation model to generate the compensated azimuth calculation model; the compensated azimuth calculation 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. The compensated azimuth angle calculation model is iteratively optimized using the gradient descent algorithm until the compensated azimuth angle calculation model converges. Then, the azimuth angle parameters are extracted to generate real-time azimuth angle parameters for cloud seeding operations.

6. The method for real-time calculation of azimuth angle for shadowing operations based on BeiDou satellite as described in claim 5, characterized in that, The compensated azimuth angle calculation model is iteratively optimized using a gradient descent algorithm until it converges. Then, azimuth angle parameters are extracted to generate real-time azimuth angle parameters for cloud seeding operations, including: Obtain the azimuth calculation error and error gradient direction for the current iteration period; The product of the gradient direction and the learning rate is calculated based on the error gradient direction, and the compensation weight coefficients are incrementally updated to adjust the compensation weight coefficients in the real-time azimuth compensation parameter set. The outputs of the displacement compensation module and curvature correction module are recalculated based on the adjusted compensation weight coefficients, and the internal state variables of the model are updated synchronously to update the compensated azimuth angle solution model. Repeat the above operations until the compensated azimuth angle calculation model meets the convergence condition. Extract the azimuth angle parameters from the converged azimuth angle calculation model and generate the real-time azimuth angle parameters for the cloud seeding operation. The convergence condition includes the error change rate of the azimuth angle calculation error being less than a set threshold or reaching the maximum number of iterations.

7. The method for real-time calculation of azimuth angle for shadowing operations based on BeiDou satellite as described in claim 6, characterized in that, The step of obtaining the azimuth calculation error and error gradient direction for the current iteration period includes: Extract the measured azimuth angle data of the weather modification operation within the current iteration cycle. The measured azimuth angle data of the weather modification operation is collected by the angle measuring equipment of the ground reference station and synchronized with the Beidou positioning data in time. The compensated azimuth angle calculation model is invoked to generate azimuth angle prediction data for the current iteration period. The azimuth angle prediction data includes the azimuth angle sequence and its confidence interval output by the compensated azimuth angle calculation model. The measured azimuth data and the predicted azimuth data are time-aligned, and the absolute error and relative error are calculated point by point to obtain the azimuth calculation error. Based on the azimuth angle calculation error, the compensation weight coefficients in the real-time azimuth angle compensation parameter set are backpropagated to generate the error gradient direction; wherein, the backpropagation process includes constructing a computation graph and taking derivatives layer by layer along the topology of the compensated azimuth angle calculation model to obtain the contribution gradient of each compensation weight coefficient to the total error.

8. The method for real-time calculation of azimuth angle for shadowing operations based on BeiDou satellite as described in claim 6, characterized in that, The step of calculating the product of the gradient direction and the learning rate based on the error gradient direction, and incrementally updating the compensation weight coefficients to adjust the compensation weight coefficients in the real-time azimuth compensation parameter set, includes: The error gradient direction is normalized to generate a gradient adjustment step size. The normalization process includes calculating the L2 norm of the gradient vector and dividing each gradient component by the norm value. The first compensation weight coefficient and the second compensation weight coefficient are updated according to the gradient adjustment step size. The updated first and second compensation weight coefficients are rewritten into the dynamic compensation weight set. During the update, spatial regularization constraints are introduced to ensure that the difference in compensation weight coefficients between adjacent positioning nodes does not exceed a set threshold, and the weight records of the most recent N iterations are retained. 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 for shadowing operations based on BeiDou satellite as described in claim 6, characterized in that, The step of extracting the azimuth parameters from the converged azimuth calculation model to generate real-time azimuth parameters for weather modification operations includes: Obtain the parameter space distribution of the converged azimuth angle solution model, wherein the parameter space distribution includes the weight matrix of the displacement compensation module and the bias vector of the curvature correction module; Calculate the eigenvectors of the parameter covariance matrix of the parameter spatial distribution, and select the top K principal components according to the size of the eigenvector values ​​to extract the azimuth principal component features; The kernel density estimation method is used to perform non-parametric probability modeling on the principal component features to construct the probability density function of the azimuth parameter; The final real-time azimuth angle parameters for the shadowing operation are determined based on the extreme points of the probability density function. These extreme points are obtained by traversing the domain of the probability density function and finding local maximum points.

10. The method for real-time calculation of azimuth angle for shadowing operations based on BeiDou satellite as described in claim 9, characterized in that, The step of using kernel density estimation to perform non-parametric probability modeling of the principal component features and constructing the probability density function of the azimuth parameter includes: The azimuth principal component features are subjected to kernel density estimation to generate an initial probability density distribution; the kernel density estimation uses a Gaussian kernel function, and the bandwidth parameter is adaptively adjusted according to the variance of the principal component features. The initial probability density distribution is calibrated based on the three-dimensional coordinate sequence in the BeiDou satellite positioning data set to generate a calibrated probability density function. The peak region of the calibrated probability density function is extracted, and the median point of the peak region is used as the real-time azimuth angle parameter for the cloud seeding operation. The peak region is obtained by binarizing the probability density function by setting a density threshold, and the median point is calculated using the region centroid algorithm.

Images