Ground laser positioning mountainous area multi-machine lifting cooperative control method and system

Abstract

The invention provides a ground laser positioning mountainous area multi-machine lifting cooperative control method and system, and relates to the technical field of cooperative control, and the method comprises the steps: 1, calculating a first real-time pose of an aircraft through a laser positioning signal transmitted by a laser positioning base station; according to the image data of the three optical identification points, calculating spatial three-dimensional coordinates of the three optical identification points to define a spatial geometric curved surface patch, and decomposing the curved surface patch into a strong signal manifold and a weak signal homotopy group; 2, taking the optical identification point of the central point of the upper edge of the front end surface of the load structure as a spatial reference origin, and connecting the two optical identification points of the left vertex and the right vertex of the lower edge of the rear end surface to form a first spatial vector and a second spatial vector; through laser positioning signal partitioning, load pose multi-source fusion optimization and dynamic compensation control of the Lyapunov potential function, the method adapts to signal shielding and attenuation scenes of mountainous complex terrains, inhibits load swing, reduces the risk of obstacle scraping and collision, and improves the safety of multi-machine lifting.

Description

Technical Field

[0001] This invention relates to the field of collaborative control technology, and in particular to a ground-based laser positioning method and system for collaborative control of multi-machine hoisting in mountainous areas. Background Technology

[0002] When conducting multi-aircraft collaborative hoisting operations in complex mountainous terrain environments, the reliability of existing ground-based collaborative control methods sometimes faces challenges.

[0003] For example, when hoisting a large box-shaped structure in a canyon with a U-shaped bend, as the load moves along the canyon's direction with the flight formation to the inside of the bend, the direct line of sight between the optical markers mounted on its side and a ground-based laser positioning base station located on the outside of the bend may be gradually obstructed by the protruding rock walls. This localized, dynamic signal obstruction caused by terrain may occur in actual operations. Existing methods typically design controllers based on the assumption that positioning signals are continuously available throughout the entire operating space. When critical markers on the load enter areas where such signals are partially obstructed or attenuated by the terrain, the accuracy of the system's calculation of its three-dimensional coordinates can easily decrease, and may even introduce jump errors. If the control system fails to effectively distinguish this change in measurement quality caused by the environment and treats it as equivalent to the measurement values ​​in areas with high-quality signals, then the real-time pose of the load calculated based on these data may not accurately reflect its true state. This could lead to the cooperative control commands generated based on this not being able to optimally suppress the load's swaying, posing a potential risk of collision with the mountain in the bend area, affecting the smoothness and safety of the flight. Summary of the Invention

[0004] The technical problem to be solved by this invention is to provide a ground laser positioning method and system for collaborative control of multi-machine hoisting in mountainous areas, which solves the problem of positioning signal obstruction / attenuation caused by mountainous terrain and improves the safety and accuracy of multi-machine hoisting.

[0005] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: Firstly, a ground-based laser positioning method for coordinated control of multi-machine hoisting operations in mountainous areas, the method comprising: Step 1: The laser positioning signal emitted by the laser positioning base station is used to calculate the first real-time pose of the aircraft; based on the image data of the three optical markers, the three-dimensional spatial coordinates of the three optical markers are calculated to define a spatial geometric surface patch, and the surface patch is decomposed into a strong signal manifold and a weak signal homotopy group. Step 2: Using the optical marker point at the center of the upper edge of the front face of the load structure as the spatial reference origin, connect the two optical marker points at the left and right vertices of the lower edge of the rear face to form the first spatial vector and the second spatial vector; Step 3: Construct a local coordinate system for the load based on the spatial reference origin and the first and second spatial vectors, and calculate the initial pose of the load; fuse the ranging data with the first real-time pose to optimize and calibrate the initial pose of the load, and obtain the optimized real-time pose of the load. Step 4: Establish a collaborative hoisting kinematics model based on the first real-time pose and the load-optimized real-time pose to calculate the initial collaborative control commands for each aircraft. Step 5: Determine the line-of-sight visibility of the load to each laser positioning base station based on the load's real-time pose optimization; construct the control Lyapunov potential function and the equipotential surface of the control Lyapunov potential function based on the topology of the strong signal manifold and the weak signal homotopy group, and divide the operating space where the load's centroid is located into a central operating region and a boundary warning region; when the load's centroid coordinates are located within the boundary warning region, generate a dynamic attitude compensation command to correct the initial cooperative control command, and output the final control command to each aircraft.

[0006] Secondly, the ground-based laser positioning multi-machine hoisting collaborative control system for mountainous areas includes: The calculation module is used to calculate the first real-time pose of the aircraft based on the laser positioning signal emitted by the laser positioning base station; based on the image data of the three optical markers, it calculates the three-dimensional spatial coordinates of the three optical markers to define the spatial geometric surface patch, and decomposes the surface patch into a strong signal manifold and a weak signal homotopy group; The decomposition module is used to connect the two optical markers at the center point of the upper edge of the front face of the load structure as the spatial reference origin, and form a first spatial vector and a second spatial vector by connecting the two optical markers at the left and right vertices of the lower edge of the rear face. The module is used to construct a local coordinate system for the load based on the spatial reference origin and the first and second spatial vectors, and to calculate the initial pose of the load; the initial pose of the load is optimized and calibrated by fusing ranging data with the first real-time pose to obtain the optimized real-time pose of the load. The calibration module is used to establish a collaborative hoisting kinematics model based on the first real-time pose and the load-optimized real-time pose, so as to calculate the initial collaborative control commands for each aircraft. The control module is used to determine the line-of-sight visibility of the load to each laser positioning base station based on the load's real-time pose optimization. Based on the topology of the strong signal manifold and the weak signal homotopy group, it constructs a control Lyapunov potential function and an equipotential surface of the control Lyapunov potential function to divide the operating space where the load's centroid is located into a central operating region and a boundary warning region. When the load's centroid coordinates are located within the boundary warning region, it generates a dynamic attitude compensation command to correct the initial cooperative control command and outputs the final control command to each aircraft.

[0007] The above-described solution of the present invention has at least the following beneficial effects: By partitioning the spatial geometric surface patch constructed from optical markers to differentiate between strong signal regions and signal obstruction / attenuation regions, the assumption of continuous availability of positioning signals throughout space is overcome. This effectively adapts to dynamically changing signal transmission environments in complex terrains such as mountainous valleys, avoiding the decrease in positioning accuracy due to terrain obstruction or signal attenuation, and ensuring the stability and effectiveness of positioning data. A local coordinate system for the load is constructed based on non-coplanar optical markers, and multi-source information is fused using real-time aircraft pose and ranging data. An iterative optimization and calibration of the load pose is achieved through a recursive state estimation algorithm, reflecting the true motion state of the load. This solves the problems of susceptibility to interference from single data sources and deviations in pose calculation, providing a basis for collaborative control. The system provides a high-precision state feedback foundation; based on the topology structure of signal partitioning, a control Lyapunov potential function is constructed to realize the regional division of the operating space and dynamic early warning. When the load approaches the weak signal area, the initial control command is corrected through dynamic attitude compensation command, which effectively suppresses load swing and reduces the potential risk of the load colliding with obstacles such as mountains due to positioning deviation or environmental interference, ensuring the flight smoothness and operational safety of multi-aircraft hoisting operations; the collaborative hoisting kinematics model fully considers the translation of the aircraft, the rotation of the load and the constraint relationship of the sling, and forms a closed-loop control with the dynamic compensation mechanism, which can adapt to the multi-aircraft formation hoisting needs in complex mountainous terrain and enhance the robustness of the system to complex working conditions such as terrain changes and signal fluctuations. Attached Figure Description

[0008] Figure 1 This is a flowchart illustrating the ground laser positioning multi-machine hoisting collaborative control method for mountainous areas provided by an embodiment of the present invention; Figure 2 This is a schematic diagram of a ground-based laser positioning multi-machine hoisting collaborative control system for mountainous areas, provided by an embodiment of the present invention. Detailed Implementation

[0009] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0010] like Figure 1 As shown, an embodiment of the present invention proposes a ground-based laser positioning method for multi-machine hoisting and collaborative control in mountainous areas. The method includes the following steps: Step 1: The laser positioning signal emitted by the laser positioning base station is used to calculate the first real-time pose of the aircraft; based on the image data of the three optical markers, the three-dimensional spatial coordinates of the three optical markers are calculated to define a spatial geometric surface patch, and the surface patch is decomposed into a strong signal manifold and a weak signal homotopy group. Step 2: Using the optical marker point at the center of the upper edge of the front face of the load structure as the spatial reference origin, connect the two optical marker points at the left and right vertices of the lower edge of the rear face to form the first spatial vector and the second spatial vector; Step 3: Construct a local coordinate system for the load based on the spatial reference origin and the first and second spatial vectors, and calculate the initial pose of the load; fuse the ranging data with the first real-time pose to optimize and calibrate the initial pose of the load, and obtain the optimized real-time pose of the load. Step 4: Establish a collaborative hoisting kinematics model based on the first real-time pose and the load-optimized real-time pose to calculate the initial collaborative control commands for each aircraft. Step 5: Determine the line-of-sight visibility of the load to each laser positioning base station based on the load's real-time pose optimization; construct the control Lyapunov potential function and the equipotential surface of the control Lyapunov potential function based on the topology of the strong signal manifold and the weak signal homotopy group, and divide the operating space where the load's centroid is located into a central operating region and a boundary warning region; when the load's centroid coordinates are located within the boundary warning region, generate a dynamic attitude compensation command to correct the initial cooperative control command, and output the final control command to each aircraft.

[0011] In this embodiment of the invention, by partitioning the spatial geometric surface patch constructed from optical markers to differentiate between strong signal regions and signal obstruction / attenuation regions, the assumption of continuous availability of positioning signals throughout space is broken. This effectively adapts to the dynamically changing signal transmission environment in complex terrains such as mountainous valleys, avoiding the decrease in positioning accuracy due to terrain obstruction or signal attenuation, and ensuring the stability and effectiveness of positioning data. Furthermore, a local coordinate system for the load is constructed based on non-coplanar optical markers. Multi-source information fusion is performed by combining real-time aircraft pose and ranging data. An iterative optimization and calibration of the load pose is achieved through a recursive state estimation algorithm, reflecting the true motion state of the load and solving the problems of susceptibility to interference from single data sources and deviations in pose calculation. This provides a high-precision state feedback basis for collaborative control; the control Lyapunov potential function is constructed based on the topology structure of signal partitioning to realize the regional division of the operating space and dynamic early warning. When the load approaches the weak signal area, the initial control command is corrected through dynamic attitude compensation command, which effectively suppresses load swing and reduces the potential risk of the load colliding with obstacles such as mountains due to positioning deviation or environmental interference, ensuring the flight smoothness and operational safety of multi-aircraft hoisting operations; the collaborative hoisting kinematics model fully considers the translation of the aircraft, the rotation of the load and the constraint relationship of the sling, and forms a closed-loop control with dynamic compensation mechanism, which can adapt to the multi-aircraft formation hoisting needs in complex mountainous terrain and enhance the robustness of the system to complex working conditions such as terrain changes and signal fluctuations.

[0012] In another preferred embodiment of the present invention, before step 1, laser positioning base stations are deployed on the ground in the canyon area, and three non-coplanar optical markers are set on the load-carrying structure: the center point of the upper edge of the front face of the load structure, the left vertex of the lower edge of the rear face, and the right vertex of the lower edge of the rear face. Each aircraft is equipped with an airborne positioning and sensing unit to acquire aircraft position data, image data containing the three optical markers, and aircraft-load ranging data. Specifically, this includes: firstly, conducting a canyon topographic survey to clarify the multi-aircraft hoisting operation route, determining the starting point, ending point, and key turning points along the route, and simultaneously measuring the slope of the mountains on both sides of the canyon, and the location and height of the protruding rock walls. Investigate potential terrain obstacles that could obstruct the signal. Based on the survey results, plan the deployment locations of laser positioning base stations. Site selection should prioritize open, unobstructed areas, and the base station deployment points must cover the entire operational flight path and the area where the load might move. The theoretical coverage radius of a single laser positioning base station is determined based on equipment parameters. The actual effective coverage distance is then calculated by considering the obstruction angle of the canyon terrain. If the obstruction angle is α and the theoretical coverage radius of a single base station is r, then the actual effective coverage distance is the product of r and cosα. Let the actual effective coverage distance be D (D = r × cosα). Based on the total length L of the operational area, calculate the required number of base stations N using the following formula: ,in This indicates rounding up, and the coverage areas of adjacent base stations must retain a certain overlap area, with an overlap rate of no less than 20%, to avoid signal blind spots. After the base stations are deployed, the parameters of each base station are calibrated, and the transmission power and signal transmission angle are adjusted to ensure that the signal can stably cover the target area.

[0013] Next, optical markers were set on the load structure. First, the front and rear faces of the load structure were defined. The front face faces the direction of flight during hoisting, and the rear face faces away from the direction of flight. The length of the upper edge of the front face was measured, and the midpoint between the left and right endpoints of the upper edge was taken as the center point of the upper edge of the front face. The coordinates of this point were calculated by adding the coordinates of the left endpoint (x_left, y_left, z_left) and the right endpoint (x_right, y_right, z_right) to the corresponding coordinate axes, then dividing by 2 to obtain the center point coordinates ((x_left + x_right) / 2, (y_left + y_right) / 2, (z_left + z_right) / 2). The length of the lower edge of the rear face was measured to determine the left and right vertices of the lower edge, and the spatial coordinates of these two points were directly read. To ensure that the three optical markers are not coplanar, spatial geometry verification was performed, and the front... The first vector from the center point of the upper edge of the face to the left vertex of the lower edge of the rear face, and the second vector from the center point of the upper edge of the front face to the right vertex of the lower edge of the rear face are cross-producted to obtain the cross-product vector. If the magnitude of this vector is not zero, it proves that the three points are not coplanar. If the magnitude is zero, the position of one of the marker points is adjusted (usually the left or right vertex of the lower edge of the rear face is selected, and small adjustments are made along the z-axis, with each adjustment controlled within 1 to 5 millimeters. The adjustment can be made upward first, and if it is still coplanar, it can be adjusted downward, or vice versa). The cross-product result is calculated again until the non-coplanar requirement is met. Finally, the three optical marker points are fixed in the determined position. When fixing, a bracket made of high-strength wear-resistant material is used to ensure that the surface of the marker points is flat, the reflective performance is uniform, and the position will not shift due to vibration or collision during the hoisting process.

[0014] Finally, each aircraft is equipped with an airborne positioning and sensing unit, which integrates three core components: a laser signal receiver, a visual sensor, and a ranging module. The laser signal receiver is installed on the top of the aircraft in an unobstructed position to receive laser positioning signals emitted by the ground laser positioning base station, thereby acquiring the aircraft's position data. The visual sensor is installed on the bottom of the aircraft, with the installation angle adjusted to ensure that its field of view can completely cover the three optical markers on the load, for acquiring image data containing the three optical markers. The ranging module is installed near the sling connection position of the aircraft, measuring the straight-line distance between the aircraft and the load using the laser ranging principle to acquire ranging data. After installation, the airborne positioning and sensing unit is calibrated to ensure that the laser signal receiver is synchronized with the ground base station signal, that the intrinsic parameters of the visual sensor (focal length, principal point coordinates, etc.) are accurate, and that the measurement error of the ranging module is controlled within the allowable range of ±2mm + measurement distance × 0.1%. At the same time, a data transmission link is established between the unit and the aircraft's flight control system to ensure that all types of data are transmitted to the flight control system in real time and stably.

[0015] This embodiment, through surveying the canyon terrain and scientifically deploying laser positioning base stations, coupled with the setting of optical markers resistant to ambient light interference, effectively improves the coverage and transmission stability of laser positioning signals in complex mountainous canyon terrain, reducing interference from mountain obstruction and terrain undulations on positioning signals. The symmetrical layout and spatial geometric design of the three non-coplanar optical markers construct a stable and clearly meaningful spatial geometric reference structure, ensuring a reliable geometric benchmark during pose calculation and improving the accuracy of load pose calculation. The airborne positioning and sensing unit integrates three major functions: laser positioning, visual capture, and distance measurement, realizing multi-source synchronous acquisition of aircraft position data, optical marker image data, and aircraft and load ranging data, enriching the dimensions of key data acquisition. The specific installation location design of the optical markers, combined with the precise acquisition capabilities of the airborne positioning and sensing unit, ensures the effective acquisition and transmission of key data required for collaborative control, reducing the risk of hoisting due to insufficient data acquisition and data deviation from the source of operation, and improving the basic reliability of multi-aircraft collaborative hoisting operations in mountainous areas.

[0016] In a preferred embodiment of the present invention, step 1 includes: Step 1.1: Process the laser positioning signal acquired by the airborne positioning and sensing unit. Solve the stochastic differential equation containing measurement noise terms to obtain the first real-time pose of the aircraft in the global coordinate system. Process the image data containing three optical markers acquired by the airborne positioning and sensing unit, extract the sub-pixel level image coordinates of each optical marker, and calculate the spatial three-dimensional coordinates of each optical marker in the camera coordinate system based on the image coordinates of the optical markers and the intrinsic parameters of the visual sensor. Specifically, this includes: first, processing the laser positioning signal received by the airborne positioning and sensing unit (including raw data such as laser propagation time and phase difference); to improve the pose calculation accuracy, establishing a stochastic differential equation containing measurement noise terms, specifically defined as follows; constructing a 9-dimensional vector covering the three-dimensional position (x, y, z) and corresponding directional velocity in the global coordinate system. , , And the three attitude angles: roll, pitch, and yaw; state equations The state transition function f(X, t) is used to quantify the dynamic evolution characteristics of the aircraft and describe the variation of each component in the state vector X with time t. Its output is a 9-dimensional vector (with the same dimension as the state vector X), and each component corresponds to the time change rate of each parameter in the state vector X. W(t) is a zero-mean Gaussian white noise vector. The noise variance is obtained through ground static calibration. The aircraft is fixed on a test platform with known coordinates, 1000 sets of data are collected, the residual between the measured value and the true value is calculated, the arithmetic mean of the residual is taken, and then the variance is solved by dividing the sum of the squares of the difference between the residual and the mean by 999. The observation equation Z(t) = H(X, t) + V(t) is used, where Z(t) is the laser measurement value (including three-dimensional position and three attitude angles), H(X, t) is the unit observation matrix (the main diagonal elements are 1 and the rest are 0, directly observing the state vector), and V(t) is the zero-mean Gaussian observation noise vector. The Euler-Markov numerical method is used to solve this stochastic differential equation. The core of this method is to achieve accurate calculation of the aircraft's attitude through iterative processes of initialization, prediction, covariance estimation, gain calculation, and state update. The specific process is as follows: First, the control cycle is set to 0.01 seconds (the standard cycle for flight control systems, balancing real-time performance and calculation accuracy). Two core parameters are initialized simultaneously: one is the state vector, which uses the initial measurements from the laser positioning signal during the first calculation, including the raw data of three-dimensional position, three-dimensional velocity, and three attitude angles; the other is the covariance matrix, based on the factory measurements from the laser sensor. The accuracy is specifically preset. It is assumed that the factory parameters specify the position measurement error as ±0.01 meters, the velocity measurement error as ±0.05 meters / second, and the attitude angle measurement error as ±0.1°. Since the correlation between the initial state variables is negligible, the covariance matrix is ​​set as a 9×9 diagonal matrix. The diagonal elements are the squares of the measurement errors of each state variable (i.e., position variance, velocity variance, and attitude angle variance, with the attitude angle unit converted to radians for uniformity of dimensions). The size of the matrix elements is positively correlated with the uncertainty of the initial state estimation. The larger the sensor measurement error, the larger the corresponding element value, thereby quantifying the reliability of the initial estimation.

[0017] After initialization, the prediction step is performed using the formula. Deducing the predicted state vector at the current moment, where, It is the current moment (the first time) The predicted state vector (period), and the optimal estimate obtained from the state update of the previous period. , refers to the final state vector output after the state update step is executed in the previous control cycle. There is no historical update result during the first solution, so the initialized state vector is used directly instead. It is the result of the state transition function calculation (substituting the optimal estimate of the previous period into the state transition function to obtain the state change rate of the previous period, such as the rate of change of position = velocity, the rate of change of velocity = acceleration, etc.). The control period is 0.01 seconds. After the prediction step is completed, covariance prediction is performed simultaneously using the formula. Quantifying the uncertainty of the predicted state, among which It is the prediction covariance matrix at the current moment, and the partial derivative Jacobian matrix. The acquisition of the values ​​must be performed according to a fixed procedure. First, the nine output elements of the state transition function f(X, t) must be identified. Then, the partial derivatives of the nine state variables are calculated for each output element. Finally, all the partial derivatives are arranged into a 9×9 matrix according to the rule that rows = output elements of the state transition function and columns = state variables. This is the covariance matrix of the previous cycle, continuing the uncertainty quantification results from the previous round; It is the transpose of the partial derivative Jacobian matrix. The process noise variance matrix is ​​constructed from the noise variance of W(t) in the state equation. The diagonal elements are the process noise variances of each state variable (such as position noise caused by airflow disturbances and velocity noise caused by dynamic system fluctuations), and the remaining elements are 0, quantifying the influence of uncontrollable noise during motion. After the covariance prediction is completed, an update step operation is performed to calculate the Kalman gain. ,in It is Kalman gain. It is the unit observation matrix (consistent with the initial observation matrix, with 1s on the main diagonal and 0s elsewhere; its function is to directly map the parameters in the state vector to the measured values, i.e., state variable = measured value, such as position state = position measured value). It is the transpose of the unit observation matrix. The observation noise variance matrix is ​​constructed from the observation noise V(t) in the observation equation, with diagonal elements representing the observation noise variance of each measurement value and the remaining elements being 0. The values ​​are determined through ground static calibration; for example, the position measurement noise variance of laser positioning is ±0.01 meters, reflecting the inherent uncertainty of the measurement data. Finally, a state update is performed using the formula... To obtain the optimal state estimate at the current moment, It is the optimal state estimate at the current moment (the final output state vector, used to calculate the aircraft's pose). All the above steps are executed in a loop in each control cycle. Through continuous iterative optimization, the prediction results and real-time measurement data are continuously integrated, and finally the three-dimensional position and three attitude angles of the aircraft in the global coordinate system are output, which is the first real-time pose.

[0018] While calculating the first real-time pose of the aircraft, image data containing three optical markers acquired by the airborne vision sensor is processed simultaneously. The specific process is as follows: First, the image is converted to grayscale using a weighted average method. The grayscale value of each pixel is equal to the sum of the red channel value multiplied by 0.299, the green channel value multiplied by 0.587, and the blue channel value multiplied by 0.114. These three coefficients, 0.299, 0.587, and 0.114, are set according to the sensitivity of the human eye to red, green, and blue colors. The human eye is most sensitive to green, so the green channel coefficient is the largest, followed by red, and then blue is the smallest. Next, interference and noise are removed by using a 5×5 convolution kernel and a Gaussian filter with a standard deviation of 1.0. The filtering process involves covering the target pixel area on the image with the convolution kernel, multiplying each element in the convolution kernel by the corresponding image pixel value, and summing all the product results to obtain the filtered grayscale value of the target pixel. The values ​​of the convolution kernel elements are calculated based on the Gaussian function to ensure that the filtering process is smooth and can effectively suppress noise.

[0019] Next, sub-pixel level coordinate extraction is performed. The Otsu method is used to determine an adaptive threshold. This method automatically finds an optimal grayscale threshold (typically 127, corresponding to a grayscale value range of 0 to 255) that maximizes the ratio of inter-class variance to intra-class variance between the foreground region (optical markers) and the background region. Based on this threshold, preliminary pixel coordinates of the optical markers are extracted, and then corrected using bilinear interpolation. The correction process involves taking four adjacent integer pixels around the preliminary coordinates (top left, top right, bottom left, and bottom right, respectively) and obtaining the grayscale values ​​of these four pixels; calculating the decimal differences (i.e., the fractional part of the preliminary coordinates) between the preliminary coordinates and the top left integer pixel; and performing horizontal interpolation, i.e., calculating the horizontal decimal difference between the grayscale values ​​of the top left and top right pixels. Weighted summation yields the upper center grayscale value; the grayscale values ​​of the lower left and lower right pixels are weighted and summed according to their horizontal decimal differences to obtain the lower center grayscale value; vertical interpolation is performed by weighted summation of the upper and lower center grayscale values ​​according to their vertical decimal differences to obtain the grayscale value at the sub-pixel position; the coordinates corresponding to the extreme grayscale values ​​are found, which are the sub-pixel level image coordinates of the optical marker; combined with the intrinsic parameters of the visual sensor (including horizontal / vertical focal length, principal point coordinates, and calibrated distortion coefficients, i.e., radial k1=-0.04, k2=0.002, k3=-0.001, tangential p1=0.001, p2=-0.0005), distortion correction is performed on the sub-pixel coordinates. First, the squared distance from the sub-pixel coordinates (u, v) to the principal point coordinates (cx, cy) is calculated. Substituting into the radial distortion correction formula ur=u×(1+k1) +k2 +k3 vr=v×(1+k1) +k2 +k3 The original coordinates are scaled and corrected using k1, k2, and k3 to eliminate radial distortion deviation, where ur and vr represent the horizontal and vertical image coordinates after radial distortion correction, respectively; then, based on the original coordinates and... Substituting into the tangential distortion compensation formula uc=ur+[2p1uv+p2( +2 )]、vc=vr+[p1( +2 [)+2p2uv], calculate the coordinate compensation amount of installation offset and manufacturing error through p1 and p2, and superimpose it on the radially corrected coordinates to finally obtain the accurate image coordinates (uc, vc) to eliminate distortion.

[0020] The physical diameter of the optical marker is preset to 0.1 meters. Its pixel diameter is obtained as 20 pixels using an image measurement tool. Dividing 0.1 meters by 20 pixels yields the ratio of physical size to pixel size. Combined with the lateral focal length fx of the visual sensor, the actual distance from the marker to the lens is calculated using the formula d = (physical diameter × fx) / (pixel diameter × ratio), meaning the actual distance equals the lateral focal length. Finally, substituting these values ​​into the inverse perspective projection formula, the camera coordinate system's lateral coordinate = (corrected lateral coordinate - principal point lateral coordinate) × actual distance ÷ lateral focal length, vertical coordinate = (corrected vertical coordinate - principal point vertical coordinate) × actual distance ÷ vertical focal length, and depth coordinate = actual distance. This allows for the calculation of the three-dimensional spatial coordinates of each optical marker in the camera coordinate system.

[0021] Step 1.2a: Using the three-dimensional spatial coordinates of the three optical markers as vertices, construct a spatial geometric surface patch. Within the spatial geometric surface patch, establish a two-dimensional parametric mesh with the first optical marker as the starting point, the second optical marker as the first side direction, and the third optical marker as the second side direction. On the two-dimensional parametric mesh, select grid points with integer x and y coordinates at intervals of 0.05 meters as sampling points to obtain twenty-one sampling points. Specifically, this includes: using the three-dimensional camera coordinates of the three optical markers obtained in Step 1.1 as vertices, construct a triangular spatial geometric surface patch (forming a uniquely determined triangular plane because the three points are not coplanar). With the first optical marker (center point of the upper edge of the front face) as the origin, and the vector pointing from the first optical marker to the second optical marker (left vertex of the lower edge of the rear face)... Using the U-axis direction as the reference direction and the vector pointing from the first optical marker point to the third optical marker point (the right vertex of the lower edge of the rear face) as the V-axis direction, a two-dimensional parametric mesh is established. The values ​​of parameters U and V are both in the range of 0 to 1. Based on the actual physical size of the mesh and the parameter range, the parameter step size is calculated. Since the sampling interval is set to 0.05 meters, the physical length in the U-axis direction (the straight-line distance from the first optical marker point to the second optical marker point) and the physical length in the V-axis direction (the straight-line distance from the first optical marker point to the third optical marker point) are calculated first. Then, the parameter step size of the U-axis and V-axis is obtained by dividing the physical length by 0.05 meters respectively. On the two-dimensional parametric mesh, grid points are selected according to the rule that the horizontal and vertical coordinates (U and V) are both integer step sizes, and finally 21 uniformly distributed sampling points are obtained to ensure that the sampling points fully cover the entire spatial geometric surface patch.

[0022] Step 1.2b: For each sampling point, calculate the geometric line vector from the sampling point to each laser positioning base station, as the theoretical line-of-sight vector; based on the first real-time pose of the aircraft and the spatial three-dimensional coordinates of the three optical markers, determine the real-time three-dimensional envelope range of the aircraft in space, and determine the spatial range occupied by the load based on the spatial positional relationship of the three optical markers. Specifically, this includes: first, obtaining the three-dimensional coordinates of the sampling points in the global coordinate system. These coordinates are obtained in the following way: based on the spatial three-dimensional coordinates of the three optical markers in the camera coordinate system obtained in step 1.1, and the first real-time pose of the aircraft (including its position and attitude in the global coordinate system), the coordinates of the three markers are uniformly transformed to the global coordinate system, and a triangular spatial geometric surface patch is constructed accordingly. Sampling points are directly generated on the two-dimensional parametric mesh defined by this surface patch.

[0023] After completing the coordinate transformation of the sampling points, the real-time three-dimensional envelope of the aircraft is constructed. The position coordinates in the first real-time pose of the aircraft are used as the center coordinates of the envelope. Combined with the preset geometric dimensions of the aircraft (2.5 meters long, 1.8 meters wide, and 1.2 meters high; these parameters are set based on the typical dimensions of a medium-sized multi-rotor aircraft, comprehensively considering the balance between payload capacity and flight maneuverability of commonly used aircraft in the industry), the corresponding dimensions are expanded by half in both positive and negative directions along the three coordinate axes of the global coordinate system. Specifically, the expansion is 1.25 meters (half of 2.5 meters) in both positive and negative directions along the X-axis, 0.9 meters (half of 1.8 meters) in both positive and negative directions along the Y-axis, and 0.6 meters (half of 1.2 meters) in both positive and negative directions along the Z-axis. This expansion method is used to construct the initial axis pairs. The rectangular prism envelope is constructed, and the coordinates of its eight vertices can be directly calculated by adding or subtracting the corresponding size offset from the center coordinates (for example, the X coordinate of the top left vertex of the front end = center X coordinate + 1.25 meters, Y coordinate = center Y coordinate + 0.9 meters, Z coordinate = center Z coordinate + 0.6 meters, and the remaining vertices are derived using the same logic). If any of the aircraft's roll angle, pitch angle, or yaw angle is not 0 (i.e., the aircraft has attitude deflection), then the 3x3 rotation matrix corresponding to the attitude angle generated above is applied to the coordinates of the eight vertices of the initial envelope, and the three-dimensional coordinates of each vertex are rotated and transformed to finally obtain a real-time three-dimensional envelope range that is consistent with the actual attitude of the aircraft. This range can accurately represent the actual area occupied by the aircraft in space.

[0024] Finally, based on the global coordinates of the three optical markers and the preset load geometric dimensions (length 1.0 m, width 0.8 m, height 0.6 m, these parameters refer to the actual size range of common mounted equipment to ensure that the envelope can completely cover the load body and connection structure), the spatial occupancy range of the load is determined. The first step is to calculate the spatial geometric center coordinates of the three optical markers. The calculation method is to add the X coordinates of the three markers and divide by 3 to obtain the X coordinate of the geometric center; add the Y coordinates of the three markers and divide by 3 to obtain the Y coordinate of the geometric center; add the Z coordinates of the three markers and divide by 3 to obtain the Z coordinate of the geometric center. These geometric center coordinates serve as the center coordinates of the load's spatial range, accurately reflecting the overall spatial position of the load. The second step is to extend the preset load geometric dimensions by half in both the forward and reverse directions along the three coordinate axes of the global coordinate system, that is, by 0.5 meters (1 / 3 of 1.0 m) in both the forward and reverse directions along the X-axis. 2) Extend the load's initial axis-aligned cuboid by 0.4 meters (half of 0.8 meters) in both directions along the Y-axis and 0.3 meters (half of 0.6 meters) in both directions along the Z-axis. The third step involves constructing a spatial plane using three marker points and calculating the normal vector of this plane (specifically, constructing two vectors using any two marker points and performing a cross product to obtain the normal vector). The direction of this normal vector directly determines the load's pitch and roll directions. The yaw direction of the load is determined by the direction of the first spatial vector (the vector pointing from the center point of the front face's upper edge to the left vertex of the rear face's lower edge), ensuring the load's attitude matches the actual lifting state. Based on this initial attitude, a corresponding 3x3 rotation matrix is ​​generated. This rotation matrix is ​​applied to the coordinates of the eight vertices of the initial axis-aligned cuboid. After completing the rotation transformation, the real-time spatial occupancy range of the load is finally determined. This range characterizes the actual area occupied by the load in space.

[0025] Suppose the three-dimensional coordinates of a vertex in the global coordinate system are (1.0 m, 0.8 m, 0.6 m). The corresponding 3x3 rotation matrix is ​​a standard rotation matrix that rotates the vertex by 30° around the Z-axis of the global coordinate system (the elements in the first row of the matrix are cos30°, -sin30°, 0; the elements in the second row are sin30°, cos30°, 0; and the elements in the third row are 0, 0, 1). During the calculation, the three-dimensional coordinates of the vertex are used as a 1x3 row vector and multiplied by the 3x3 rotation matrix. Specifically, the new X-coordinate = 1.0 m × cos30° + 0.8 m. ×(-sin30°)+0.6 m×0, new Y coordinate = 1.0 m×sin30°+0.8 m×cos30°+0.6 m×0, new Z coordinate = 1.0 m×0+0.8 m×0+0.6 m×1. The new coordinates of the vertex after attitude deflection are obtained through this calculation. The new coordinates of each vertex after attitude deflection are obtained through calculation. After all vertices have completed the rotation transformation, a cuboid structure that perfectly matches the actual attitude of the load is formed. Finally, the real-time spatial occupancy range of the load is determined. This range represents the actual area occupied by the load in space.

[0026] Step 1.2c: Determine whether each theoretical line-of-sight vector intersects with the real-time 3D envelope of the aircraft or the spatial range of the load. If the theoretical line-of-sight vector intersects with the spatial range of the aircraft or the load, it is determined that the corresponding theoretical line-of-sight vector is occluded; otherwise, it is determined that the corresponding theoretical line-of-sight vector is not occluded, and the occlusion judgment result is obtained. Specifically, the ray-cuboid intersection detection algorithm is used to determine whether each theoretical line-of-sight vector is occluded. The specific operation is as follows: the theoretical line-of-sight vector is defined as a ray from the laser positioning base station to the sampling point (the starting point of the ray is the global coordinate of the base station, and the ending point is the global coordinate of the sampling point). At the same time, the real-time 3D envelope of the aircraft and the load are regarded as two independent cuboid obstacles. First, the intersection point of the ray with the six faces of the 3D envelope of the aircraft is calculated. The plane equation of the six faces of the 3D envelope of the aircraft is obtained first. Based on the vertex coordinates of the envelope, each face is determined by four vertices. The plane equation is in the form of ax + by + cz + d = 0, where a, b, and c are the three components of the plane normal vector. By taking any three non-collinear vertices in the face to construct two vectors, the cross product operation of the vectors can be performed to obtain the result. For the plane intercept parameter, substitute the coordinates (x0, y0, z0) of any vertex in the plane into the equation, and then... =-(ax0+by0+cz0) is calculated, and its value is related to the distance from the plane to the origin of the coordinate system. Then, the ray parameter equation (let the ray origin be P0, the direction vector be e, the parameter t≥0, and the ray equation be P(t)=P0+t×e) is substituted into the plane equation of each surface to solve for the value of parameter t, where P(t) represents the three-dimensional coordinate point of the ray in space corresponding to parameter t. After calculating the coordinates of the corresponding intersection point based on the t value, it is determined whether the intersection point falls within the current boundary range (i.e., within the valid interval of the x, y, z coordinates of the surface), and the t value satisfies 0≤t≤1 (ensuring that the intersection point is located within the ray segment from the base station to the sampling point), satisfying the above conditions. The intersection of the ray and the aircraft's envelope is the valid interior point. If such a valid interior point exists, the theoretical line-of-sight vector is directly determined to be occluded by the aircraft, and the detection process for the next vector begins. If the ray does not have a valid intersection with the aircraft's envelope, the load occlusion judgment continues. The intersection of the ray with the six faces of the load's three-dimensional envelope is calculated using the same logic. If a valid intersection exists within the ray segment, the theoretical line-of-sight vector is determined to be occluded by the load. If the ray does not have a valid intersection with either the aircraft's or the load's three-dimensional envelope, the theoretical line-of-sight vector is determined to be unoccluded. Finally, the occlusion status (occluded / unoccluded) of all theoretical line-of-sight vectors is recorded, and a complete set of occlusion judgment results is compiled.

[0027] Step 1.2d: For the theoretical line-of-sight vector that is not obstructed, based on the transmission power of the laser positioning signal, the atmospheric attenuation characteristics of the signal along the transmission path, and the receiver noise figure, evaluate the signal-to-noise ratio (SNR) level of the received signal along the corresponding path to obtain the SNR evaluation result. Specifically, for the theoretical line-of-sight vector determined to be unobstructed in Step 1.2c, evaluate the SNR level according to the following process and collect the key parameters required for evaluation. The rated transmission power of the laser positioning base station is preset to 5 watts (set according to the signal coverage requirements of the complex environment of mountainous canyons, ensuring sufficient power redundancy to ensure signal transmission while also balancing equipment energy consumption). The atmospheric attenuation coefficient of the laser signal in the mountainous canyon environment is taken as 0.5 dB / km. This value is calibrated through experimental measurement, using the typical range of 0.2 to 0.8 dB / km for the aerosol extinction coefficient in mountainous areas as a benchmark. Refer to the atmospheric attenuation coefficient standard to calibrate environmental parameters (25℃, 50% humidity, 0.1 mg / m³ particulate matter concentration), and include the temperature and humidity, dust and vegetation debris in the target area with medium humidity and low particulate matter concentration. The measured data was compared with the actual data; combined with the characteristics of the 1.55μm near-infrared band commonly used in laser positioning (this band is easily affected by water vapor absorption and aerosol scattering in mountainous areas, and the attenuation law conforms to industry calibration specifications), the coefficient was quantitatively adjusted according to a fixed ratio, that is, for every 5% increase in humidity above the standard value, the coefficient was increased by 0.05dB / km; for every 0.02mg / m³ decrease in particulate matter concentration below the standard value, the coefficient was decreased by 0.05dB / km, and the final value was determined; the inherent loss of the system was preset to 3dB (including antenna insertion loss and link matching loss at the transmitter and receiver, determined with reference to the typical loss value of similar laser positioning systems); the receiver noise figure was preset to 2dB (selecting the conventional indicators of industrial-grade laser receivers to ensure low-noise reception performance); the absolute temperature of the receiver's operating environment was set to 300 Kelvin (corresponding to a normal temperature environment of 27 degrees Celsius, which is in line with the temperature range of most outdoor operations); the signal bandwidth was preset to 10 MHz (matching the transmission rate of the laser positioning signal, while improving anti-interference capability); the Boltzmann constant was a fixed constant, with a value of 1.38 × Joules / Kelvin.

[0028] The signal transmission path length is determined by the theoretical line-of-sight vector length. First, the atmospheric attenuation (the product of the atmospheric attenuation coefficient and the transmission path length) is calculated. Then, the received signal power is calculated using the formula: Received signal power = 5 watts - Atmospheric attenuation coefficient × Transmission path length - 3 dB. This yields the final received signal power. For noise power calculation, the thermal noise power is first calculated using the formula (Thermal noise power = Boltzmann constant × Absolute temperature of receiver operating environment × Signal bandwidth). Then, the total noise power is calculated by combining this with the receiver noise figure. The specific formula is: Total receiver noise power = 1.38 × Joules / Kelvin × 300 Kelvin × Hertz × 2, the total noise power of the receiver is obtained through this formula; the obtained received signal power and the total noise power of the receiver are substituted into the signal-to-noise ratio (SNR) calculation formula (SNR = 10 × 2). (Received signal power / Total noise power of receiver)) Calculate the signal-to-noise ratio (SNR) value for the corresponding transmission path, and finally summarize the SNR results of all unobstructed vectors to form a complete SNR evaluation result set.

[0029] Step 1.2e: Summarize the occlusion judgment results and signal-to-noise ratio (SNR) evaluation results of all sampling points for all laser positioning base stations, as the occlusion judgment results and SNR level evaluation results of the theoretical line-of-sight vector. Specifically, this includes: integrating the associated data of each sampling point with all laser positioning base stations, classifying them by sampling point number and base station number, and storing the corresponding occlusion judgment results (occluded / unoccluded) and SNR evaluation results (specific values) in association, forming a complete theoretical line-of-sight vector occlusion and SNR evaluation dataset. This dataset fully records the transmission quality information of each sampling point under the signal coverage of different base stations.

[0030] Step 1.3: Based on the occlusion judgment results of the theoretical line-of-sight vector and the evaluation results of the signal-to-noise ratio (SNR), the continuous point set regions on the spatial geometric surface patch with two or more unobstructed lines of sight and an average SNR higher than a preset threshold are classified as strong signal manifolds; the point set regions with line-of-sight occlusion or an average SNR lower than a preset threshold are classified as weak signal homotopy groups. Specifically, this includes: classifying strong signal manifolds and weak signal homotopy groups based on the evaluation results. The specific operation is as follows: the preset SNR threshold is calibrated to 20 dB, and this value is determined through joint ground static testing and actual flight experiments. The specific process is as follows: In a flat, open area simulating the signal propagation characteristics of mountain valleys, a laser positioning base station and an experimental aircraft platform equipped with a high-precision GPS (positioning accuracy ±0.01 meters) and a laser receiver are deployed. The distance between the base station and the aircraft is fixed (50 to 500 meters), and different signal-to-noise ratio scenarios are simulated by adjusting the signal attenuator. After receiving the laser signal emitted by the base station, the aircraft synchronously collects the time-of-flight (TOA) data of the laser signal. Combined with the known coordinates of the base station, the real-time position of the aircraft is calculated using a triangulation algorithm, and then integrated with the inertial measurement unit on board the aircraft. Attitude angles (roll, pitch, yaw) were calculated using IMU data, and 100 sets of attitude calculation data were collected under different signal-to-noise ratios (SNRs). Test results showed that when the SNR was ≥20dB, the average position calculation error was ≤0.08 meters and the average attitude calculation error was ≤0.3°, both exceeding the core performance requirements for aircraft attitude control (position error ≤0.1 meters, attitude error ≤0.5°). When the SNR was <20dB, the calculation error increased significantly, reaching an average position error of 0.15 meters at 18dB, exceeding the control accuracy threshold. In a target mountain canyon scenario (sea...),... Flight tests were conducted with an elevation drop of 300 meters and a vegetation coverage of 60%. The aircraft flew along a preset route (altitude 50 to 200 meters and speed 3 to 8 m / s). Laser signals from multiple base stations were collected in real time through a laser receiver, and signal-to-noise ratio data and pose calculation process were recorded simultaneously (based on TOA+IMU fusion algorithm). The results verified that the 20dB threshold can effectively resist the influence of background noise such as airflow interference and vegetation scattering in mountainous areas. During the continuous flight of 2 hours, the time period with a signal-to-noise ratio ≥20dB accounted for 85%, and there were no frame drops or jumps in the calculated data, which met the requirements for continuous control.

[0031] Comparing the three candidate thresholds of 18dB, 20dB, and 22dB, 18dB improves signal coverage but suffers from insufficient stability in computational accuracy; 22dB further reduces errors but reduces signal coverage by 15%, resulting in some areas lacking effective signals; 20dB achieves the optimal balance between accuracy and coverage, and was ultimately determined as the preset threshold. For all sampling points on the spatial geometric surface patch, the number of corresponding unobstructed theoretical line-of-sight vectors is counted, and the arithmetic mean of the signal-to-noise ratio of all unobstructed lines of sight is calculated to determine whether the condition is met. For sampling points that meet the conditions of having ≥2 unobstructed lines of sight and an average signal-to-noise ratio >20dB, their spatial continuity (physical distance between adjacent sampling points <0.1 meters) is assessed, and the continuous point set region that meets the conditions is classified as a strong signal manifold. For sampling points that do not meet the conditions of having ≥2 unobstructed lines of sight or an average signal-to-noise ratio ≤20dB, regardless of their spatial continuity, their point set region is classified as a weak signal homotopy group. The spatial coordinate ranges of the strong signal manifold and the weak signal homotopy group are recorded, and a signal quality partition map is generated.

[0032] This embodiment effectively improves the accuracy and stability of pose and coordinate calculation by solving stochastic differential equations that include measurement noise terms, and combines sub-pixel-level image coordinate extraction and visual intrinsic parameter fusion to calculate the three-dimensional coordinates of optical markers, reducing the impact of environmental interference and measurement errors. The construction and uniform sampling of spatial geometric surface patches enable comprehensive coverage detection of signal transmission paths in key areas of the load. The occlusion judgment and signal-to-noise ratio evaluation of theoretical line-of-sight vectors can identify signal attenuation scenarios caused by terrain or self-occlusion, providing a quantitative basis for signal quality zoning. The division of strong signal manifolds and weak signal homotopy groups enables dynamic perception and classification of signal quality in complex mountainous environments, helping to avoid pose calculation deviations caused by low-quality signals, reducing the risk of load swaying and collisions, and improving the safety and smoothness of multi-aircraft hoisting operations.

[0033] In a preferred embodiment of the present invention, step 2 includes: Step 2.1: From the three-dimensional spatial coordinates of the three optical markers, select the coordinates of the optical marker located at the center of the edge of the front face of the load structure, and define it as the coordinates of the spatial reference origin. Specifically, this includes: First, using a laser positioning system, acquiring the spatial coordinates of the three preset optical markers on the load structure to ensure that each marker can obtain complete three-dimensional coordinate data, including specific measurement values ​​in the x, y, and z axes. Simultaneously, verifying the accuracy of the measurement data to achieve millimeter-level measurement standards is required to meet the high-precision spatial positioning requirements of the load structure. Next, according to the preset load structure marking rules, verify the attributes of the optical markers. These preset rules clearly define the installation position parameters of the three optical markers, including the mounting surface (front face / rear face), edge position (upper edge / lower edge), and vertex. The type (center point / left vertex / right vertex) is determined by comparing the installation drawings of the load structure, the physical markings on the surface of the marker points, or the marker information stored by the measuring equipment. The actual physical position of each optical marker point is confirmed one by one. Finally, the three marker points are identified as the center point of the upper edge of the front face of the load structure, the left vertex of the lower edge of the rear face of the load structure, and the right vertex of the lower edge of the rear face of the load structure, respectively. Finally, from the three optical marker points whose attributes have been verified, the marker point whose physical position is the center point of the upper edge of the front face of the load structure is selected. The x-axis coordinate measurement value of the marker point is directly assigned to the x-axis coordinate of the spatial reference origin, the y-axis coordinate measurement value is assigned to the y-axis coordinate of the spatial reference origin, and the z-axis coordinate measurement value is assigned to the z-axis coordinate of the spatial reference origin. The three-dimensional coordinate data of the marker point is completely used to complete the definition of the coordinates of the spatial reference origin.

[0034] Step 2.2: Calculate the difference between the spatial 3D coordinates of the optical marker point located at the left vertex of the lower edge of the rear face and the coordinates of the spatial reference origin to obtain the first spatial vector; calculate the difference between the spatial 3D coordinates of the optical marker point located at the right vertex of the lower edge of the rear face and the coordinates of the spatial reference origin to obtain the second spatial vector. Specifically, this includes: calculating the first spatial vector by extracting the spatial 3D coordinates of the determined optical marker point at the left vertex of the lower edge of the rear face, and recording its x-axis, y-axis, and z-axis coordinate values; extracting the x-axis, y-axis, and z-axis coordinate values ​​of the spatial reference origin defined in Step 2.1; subtracting the x-axis coordinate value of the spatial reference origin from the x-axis coordinate value of the left vertex of the lower edge of the rear face to obtain the x-axis component of the first spatial vector; subtracting the y-axis coordinate value of the spatial reference origin from the y-axis coordinate value of the left vertex of the lower edge of the rear face to obtain the y-axis component of the first spatial vector; and subtracting the z-axis coordinate value of the spatial reference origin from the z-axis coordinate value of the left vertex of the lower edge of the rear face. The z-axis coordinate value is used to obtain the z-axis component of the first spatial vector. The calculated x-axis, y-axis, and z-axis components are combined to form the complete first spatial vector. For the calculation of the second spatial vector, the spatial three-dimensional coordinates of the determined optical marker point at the right vertex of the lower edge of the rear face are extracted, and its x-axis, y-axis, and z-axis coordinate values ​​are recorded respectively. Using the x-axis, y-axis, and z-axis coordinate values ​​of the spatial reference origin defined in step 2.1, the x-axis coordinate value of the right vertex of the lower edge of the rear face is subtracted from the x-axis coordinate value of the spatial reference origin to obtain the x-axis component of the second spatial vector. The y-axis coordinate value of the right vertex of the lower edge of the rear face is subtracted from the y-axis coordinate value of the spatial reference origin to obtain the y-axis component of the second spatial vector. The z-axis coordinate value of the right vertex of the lower edge of the rear face is subtracted from the z-axis coordinate value of the spatial reference origin to obtain the z-axis component of the second spatial vector. The calculated x-axis, y-axis, and z-axis components are combined to form the complete second spatial vector.

[0035] This embodiment avoids calculation errors caused by fuzzy reference by clearly defining the selection criteria for the spatial reference origin; the standardized coordinate component subtraction operation process ensures the calculation accuracy of the first spatial vector and the second spatial vector, reflects the relative spatial relationship between key marker points of the load structure, and helps to improve the spatial positioning accuracy and control stability of the aircraft load system.

[0036] In a preferred embodiment of the present invention, step 3 includes: Step 3.1: Normalize the first spatial vector and use it as the first coordinate axis of the load local coordinate system. Calculate the dot product of the second spatial vector and the first coordinate axis. Subtract the product of the dot product and the first coordinate axis from the second spatial vector to obtain an intermediate vector orthogonal to the first coordinate axis. Normalize this intermediate vector and use it as the second coordinate axis of the load local coordinate system. Specifically, this includes: Normalizing the first spatial vector: First, calculate the magnitude of the first spatial vector. The magnitude is calculated as the square root of the sum of the squares of the x-axis component, the y-axis component, and the z-axis component. Then, divide the x-axis component, y-axis component, and z-axis component of the first spatial vector by the calculated magnitude to obtain three standardized components. These three components together constitute the normalized first spatial vector, which is then used as the first coordinate axis of the load local coordinate system. Calculate the dot product of the second spatial vector and the first coordinate axis (the normalized first spatial vector). The calculation process involves multiplying the x-axis component of the second spatial vector by the first coordinate axis. The x-axis component of the coordinate axis is multiplied by the y-axis component of the second spatial vector multiplied by the y-axis component of the first coordinate axis, and then multiplied by the z-axis component of the second spatial vector multiplied by the z-axis component of the first coordinate axis, resulting in a scalar dot product. Each component of the second spatial vector is then subtracted from the product of the dot product and the corresponding component of the first coordinate axis. Specifically, the x-axis component of the second spatial vector is subtracted from the dot product multiplied by the first x-axis component, the y-axis component of the second spatial vector is subtracted from the dot product multiplied by the first y-axis component, and the z-axis component of the second spatial vector is subtracted from the dot product multiplied by the first z-axis component, resulting in three new components. These are combined to form an intermediate vector orthogonal to the first coordinate axis. Using the same method as for normalizing the first spatial vector, the magnitude of the intermediate vector is first calculated (the square root of the sum of the squares of the x-axis, y-axis, and z-axis components of the intermediate vector). Then, each of the three components of the intermediate vector is divided by this magnitude to obtain the normalized intermediate vector, which is then used as the second coordinate axis of the local coordinate system.

[0037] Step 3.2: Calculate the cross product of the first and second coordinate axes to obtain the third coordinate axis of the load local coordinate system. The load local coordinate system is defined by the spatial reference origin, the first coordinate axis, the second coordinate axis, and the third coordinate axis. Specifically, when calculating the third coordinate axis, the cross product of the first and second coordinate axes is calculated. The specific calculation process for the cross product is as follows: the x-axis component of the third coordinate axis equals the y-axis component of the first coordinate axis multiplied by the z-axis component of the second coordinate axis, minus the z-axis component of the first coordinate axis multiplied by the y-axis component of the second coordinate axis; the y-axis component of the third coordinate axis equals the z-axis component of the first coordinate axis multiplied by the x-axis component of the second coordinate axis, minus the first... The x-axis component of the first coordinate axis is multiplied by the z-axis component of the second coordinate axis; the z-axis component of the third coordinate axis is equal to the x-axis component of the first coordinate axis multiplied by the y-axis component of the second coordinate axis, minus the y-axis component of the first coordinate axis multiplied by the x-axis component of the second coordinate axis. Through the above calculation, the three components of the third coordinate axis are obtained. This coordinate axis is naturally orthogonal to the first two coordinate axes. Taking the spatial reference origin determined in step 2.1 as the origin of the coordinate system, the first coordinate axis obtained in step 3.1 as the x-axis of the local coordinate system, the second coordinate axis obtained in step 3.1 as the y-axis of the local coordinate system, and the third coordinate axis calculated in this step as the z-axis of the local coordinate system, the three together constitute the complete load local coordinate system, which is a right-handed orthogonal coordinate system.

[0038] Step 3.3: Based on the rotation matrix and translation vector between the load's local coordinate system and the global coordinate system, calculate the initial position and initial attitude of the load in the global coordinate system to obtain the initial pose of the load. Specifically, this includes: determining the rotation matrix and translation vector; the rotation matrix between the load's local coordinate system and the global coordinate system is composed of the unit vectors of the three coordinate axes of the load's local coordinate system. The first column of the rotation matrix corresponds to the unit vector component of the first coordinate axis, the second column corresponds to the unit vector component of the second coordinate axis, and the third column corresponds to the unit vector component of the third coordinate axis. The translation vector is the three-dimensional coordinate of the spatial reference origin defined in Step 2.1 in the global coordinate system, that is, the x-axis, y-axis, and z-axis coordinates of the spatial reference origin are directly used as the corresponding components of the translation vector; initial position calculation, that is, the three-dimensional coordinate data contained in the translation vector are directly used as the initial position of the load in the global coordinate system, that is, the x-axis coordinate of the load's initial position is equal to the x-axis component of the translation vector, the y-axis coordinate is equal to the y-axis component of the translation vector, and the z-axis coordinate is equal to the z-axis component of the translation vector; based on the preceding text... The constructed rotation matrix (composed of three orthogonal unit vectors in the local coordinate system of the load) is used to solve the initial attitude of the load using the standard ZYX sequential Euler angle transformation algorithm. First, the pitch angle is solved using the element in the third row of the first column of the rotation matrix (this element is equal to the sine of the pitch angle, which is obtained by arcsine calculation). Then, the roll angle is solved by combining the elements in the first row and the third row of the first column of the rotation matrix (using the cosine of the pitch angle as the denominator, and dividing the difference between the cosine of the element in the first row of the first column and the sine of the element in the third row of the first column by the denominator). Then, the roll angle is obtained through arctangent calculation; finally, the yaw angle is solved by combining the elements of the third row and the third column of the rotation matrix and the third row of the second column (using the cosine of the pitch angle as the denominator, the difference between the cosine of the third row and the sine of the second row is divided by the denominator and then the yaw angle is obtained through arctangent calculation). Through the above steps, the specific values ​​of the three attitude angles are solved in sequence, accurately representing the attitude state of the load in the global coordinate system. Then, it is integrated with the initial position of the load calculated above to finally form complete initial position and pose data of the load.

[0039] Step 3.4a: Receive the initial pose of the load as the initial state vector of the recursive state estimation algorithm; receive the first real-time pose as the known input condition for the recursive state estimation algorithm; receive the ranging data between the aircraft and the load measured by the airborne positioning and sensing unit as the observations for the recursive state estimation algorithm. Specifically, this includes: using the initial pose of the load obtained in step 3.3 as the initial state vector of the recursive state estimation algorithm. This vector contains the x-axis, y-axis, and z-axis coordinate components of the initial position of the load, as well as the roll angle, pitch angle, and yaw angle components of the initial attitude, ensuring that the algorithm starts... It has an initial state reference; it receives the first real-time pose as the known input condition for the recursive state estimation algorithm. The first real-time pose is the real-time pose data of the aircraft, including the three-dimensional coordinates of the real-time position of the aircraft and the three Euler angles of the real-time attitude, providing an external motion reference for the algorithm to predict the load state; it receives the distance measurement data between the aircraft and the load measured by the airborne positioning and sensing unit as the observation of the recursive state estimation algorithm. This distance measurement data is the straight-line distance value between the aircraft and the load collected in real time. The collection frequency is consistent with the algorithm running frequency to ensure the timeliness of the observation data.

[0040] Step 3.4b: In the prediction step of the recursive state estimation algorithm, the motion of the aircraft and the load itself are described by the first real-time pose. In the update step of the recursive state estimation algorithm, the residual between the predicted value and the actual measured value of the relative ranging observation at the current moment is calculated, and the Kalman gain is calculated based on the state covariance matrix and the observation noise covariance matrix. Specifically, to ensure the timing synchronization between the recursive state estimation algorithm and the aircraft flight control system, and to balance real-time performance and solution accuracy, the control period of the algorithm is set to 0.01 seconds. Based on this control period, the load state vector is first defined and the state transition function is constructed, and then the state change rate is calculated. The load state vector is defined as a 6-dimensional vector, denoted by the symbol E, i.e. ,in, This represents the matrix transpose operation. This represents the horizontal position of the load in the global coordinate system. Represents vertical position. φ represents the depth position; φ represents the roll angle of the load; θ represents the pitch angle; and ψ represents the yaw angle. These six components together completely describe the spatial state of the load. The state transition function, denoted as f(X, t), is used to quantify the dynamic characteristics of the load. Its mathematical expression is the set of rates of change of each component of the load state vector, i.e., f(X, t) = In the formula, Corresponding to the rate of change of lateral position, Corresponding to the rate of change of longitudinal position, The corresponding depth position change rate; φ̇ corresponds to the roll angle change rate, θ̇ corresponds to the pitch angle change rate, and ψ̇ corresponds to the yaw angle change rate; specifically, the load position change rate is equal to the load's velocity, which is derived from the tension and sway characteristics of the sling. Specifically, it is obtained by dividing the load position change at the previous moment by the control cycle to get the velocity at the previous moment, and then combining it with the aircraft's motion trend correction to obtain the current position change rate, i.e., lateral position change rate = (lateral position at the previous moment - lateral position at the moment before that) ÷ control cycle + lateral velocity correction, longitudinal position change rate = (longitudinal position at the previous moment - longitudinal position at the moment before that) ÷ control cycle + longitudinal velocity correction, depth position change rate = (depth position at the previous moment - depth position at the moment before that) ÷ control cycle + depth velocity correction (where the lateral position, longitudinal position, and depth position at the previous moment are the values ​​of the previous moment's lateral position, longitudinal position, and depth position). The three-dimensional position of the load at a given moment, the lateral position of the moment before that, the longitudinal position of the moment before that, and the depth position of the moment before that are the three-dimensional positions of the load at the moment before that. The control period is 0.01 seconds. The lateral velocity correction, longitudinal velocity correction, and depth velocity correction are velocity corrections based on the aircraft's motion trend. The rate of change of the load's attitude angle is equal to the angular velocity of the load. This angular velocity is calculated by the change in the spatial position of three optical markers. That is, the rotation angle is constructed by the coordinate difference of the markers at adjacent moments, and then the rotation angle is divided by the control period. That is, the roll angle change rate = roll angle change ÷ control period, the pitch angle change rate = pitch angle change ÷ control period, and the yaw angle change rate = yaw angle change ÷ control period (where the roll angle change, pitch angle change, and yaw angle change are attitude angle changes constructed by the coordinate difference of the optical markers at adjacent moments, and the control period is 0.01 seconds).

[0041] After obtaining the state change rate, the load state vector at the current moment is predicted. Each component of the predicted state vector at the current moment is calculated using the formula: the corrected component from the previous moment + the corresponding component's state change rate × control period. For example, the predicted load lateral position at the current moment = the corrected load lateral position from the previous moment + the load lateral position change rate × 0.01 seconds; the predicted load roll angle at the current moment = the corrected load roll angle from the previous moment + the load roll angle change rate × 0.01 seconds. The remaining components (longitudinal position, depth position, pitch angle, yaw angle) are calculated one by one using the same logic to obtain the complete predicted state vector. After predicting the state vector, the partial derivative Jacobian matrix is ​​calculated. This matrix is ​​a 6x6 matrix (the state vector contains 6 components: 3D position and 3 attitude angles). Each element in the matrix is ​​the partial derivative of the corresponding output component of the state transition function with respect to the corresponding input component of the state vector. For example, the element in the first row and first column is negative. The partial derivative of the load lateral position change rate with respect to the load lateral position; the element in the fourth row and first column = the partial derivative of the load roll angle change rate with respect to the load lateral position; calculate all 36 elements in sequence according to the above rules, and then arrange them in the order of the row corresponding to the state transition function output and the column corresponding to the state vector input to form a 6x6 partial derivative Jacobian matrix; after completing the construction of the partial derivative Jacobian matrix, further predict the state covariance matrix. The calculation process is as follows: prediction covariance matrix = partial derivative Jacobian matrix × previous corrected state covariance matrix × transpose of partial derivative Jacobian matrix + process noise variance matrix; where the process noise variance matrix is ​​a 6x6 diagonal matrix, the diagonal elements are the process noise variance of each state component (obtained through ground static calibration, i.e., 1000 sets of static data are collected under a fixed load, the residual between each set of data and the true value is calculated, and then the variance of the residual is calculated), and the off-diagonal elements are all 0, which is used to quantify the influence of uncontrollable noise in motion.

[0042] After calculating the load state vector and state covariance matrix in the prediction step, the update step is performed to calculate the residual and solve the Kalman gain. First, the ranging prediction value is calculated. Based on the load position coordinates in the predicted state vector and the aircraft position coordinates in the first real-time pose, the calculation is performed step-by-step using the distance formula between two points in space. The first step calculates the difference between the aircraft's lateral position and the load's lateral position, and then squares this difference. The second step calculates the difference between the aircraft's longitudinal position and the load's longitudinal position, and then squares this difference. The third step calculates the difference between the aircraft's depth position and the load's depth position, and then squares this difference. The fourth step adds the three squared results and performs a square root operation on the sum to obtain the final ranging prediction value. Based on the obtained ranging prediction value, the residual is further calculated. The residual = airborne positioning and... The actual ranging value measured by the sensing unit is divided by the predicted ranging value. The sign of the residual reflects the direction of the deviation between the predicted and actual values, and the absolute value reflects the magnitude of the deviation. Then, an observation matrix is ​​constructed. This matrix is ​​a 1x6 matrix, and the matrix elements are the partial derivatives of the predicted ranging value with respect to each component of the state vector. Specifically, the partial derivative of the predicted ranging value with respect to the lateral position of the load is (lateral position of the aircraft - lateral position of the load) ÷ predicted ranging value; the partial derivative of the predicted ranging value with respect to the longitudinal position of the load is (longitudinal position of the aircraft - longitudinal position of the load) ÷ predicted ranging value; and the partial derivative of the predicted ranging value with respect to the depth position of the load is (depth position of the aircraft - depth position of the load) ÷ predicted ranging value. Since the attitude angle does not directly affect the straight-line distance between two points, the partial derivatives of the predicted ranging value with respect to the three attitude angles are all 0. The above results are arranged to form a 1x6 observation matrix.

[0043] Based on the constructed observation matrix, the observation noise covariance matrix is ​​determined. This matrix is ​​a 1x1 matrix, and its elements are the variance of the ranging measurement noise (obtained by acquiring the known distance between the fixed aircraft and the payload through ground static calibration, collecting 1000 sets of ranging data, calculating the residual between each set of data and the true distance, and then calculating the variance of the residual). Finally, the Kalman gain is calculated in three steps. The first step is to calculate intermediate matrix 1 = predicted state covariance matrix × transpose of observation matrix (resulting in a 6x1 matrix). The second step is to calculate intermediate matrix 2 = observation matrix × predicted state covariance matrix × transpose of observation matrix + observation noise covariance matrix (resulting in a 1x1 matrix). The third step is to obtain the Kalman gain = intermediate matrix 1 ÷ intermediate matrix 2 based on the calculation results of the first two steps (resulting in a 6x1 matrix, where each element corresponds to the correction weight of a state component).

[0044] Step 3.4c: Use Kalman gain to correct and update the load state vector and state covariance matrix obtained in the prediction step, obtaining the corrected load state vector and state covariance matrix at the current time. Output the position and attitude information contained in the corrected load state vector as the load-optimized real-time pose. Specifically, after obtaining the Kalman gain, use the Kalman gain to correct and update the load state vector and state covariance matrix obtained in the prediction step. First, correct the load state vector. The correction value of each state component is calculated by the predicted value + the corresponding component of the Kalman gain × the residual. For example, the corrected load lateral position = the predicted value of the load lateral position + the component corresponding to the lateral position in the Kalman gain × the residual; the corrected load roll angle = the predicted value of the load roll angle + the component corresponding to the roll angle in the Kalman gain × the residual; the other four components (longitudinal position, depth position, pitch angle, yaw angle) are calculated one by one according to the same logic. Finally, the corrected load state vector at the current time is obtained. This vector eliminates... This process reduces some prediction errors and more closely approximates the actual load state. After correcting the load state vector, the state covariance matrix is ​​further corrected in two steps. The first step calculates the correction matrix = 6x6 identity matrix - Kalman gain × observation matrix (diagonal elements of the identity matrix are 1, and off-diagonal elements are 0). The second step calculates the corrected state covariance matrix = correction matrix × predicted state covariance matrix. The values ​​of the diagonal elements of this matrix reflect the uncertainty of the corresponding state components after correction; the smaller the value, the higher the estimation accuracy. After correcting the load state vector and state covariance matrix, the load-optimized real-time pose is output. Steps 3.4b and 3.4c are repeated in each control cycle to continuously iterate and optimize the pose estimation accuracy. The three-dimensional position information (lateral coordinates, longitudinal coordinates, and depth coordinates) and three attitude angle information (roll angle, pitch angle, and yaw angle) in the load state vector after each correction are integrated to form the load-optimized real-time pose, which is then output to the subsequent cooperative control stage in real time.

[0045] This embodiment ensures the accuracy and orthogonality of the load's local coordinate system construction by normalizing and orthogonalizing the spatial vectors, providing a reliable coordinate system foundation for load pose estimation. The initial pose estimation method based on rotation matrices and translation vectors achieves accurate determination of the load's initial position and attitude, providing high-quality initial input for the recursive state estimation algorithm. Utilizing real-time aircraft pose and ranging observation data, dynamic estimation of the load pose is achieved through iterative calculations of prediction and update steps. Kalman gain enables optimal fusion of predicted state and observation data, effectively suppressing the influence of measurement noise and model errors, and improving the accuracy and stability of load pose estimation.

[0046] In a preferred embodiment of the present invention, step 4 includes: Step 4.1: Describe the motion state of each aircraft using the first real-time pose, and describe the rigid body motion state of the load using the load-optimized real-time pose. Combined with the geometric relationship between the sling length and connection point, establish a collaborative hoisting kinematics model that includes aircraft translation, load rotation, and sling constraints. Specifically, to accurately characterize the dynamics of the multi-aircraft collaborative hoisting system, the model is established using the first real-time pose of each aircraft (including lateral position, longitudinal position, depth position, roll angle, pitch angle, and yaw angle in the global coordinate system) as the aircraft motion state input, and the load-optimized real-time pose (including lateral position, longitudinal position, depth position, roll angle, pitch angle, and yaw angle in the global coordinate system) as the input. The input is the rigid body motion state of the load. It also incorporates the preset sling length (assumed to be 2.0m, obtained through ground calibration, based on a load weight of 50kg, dimensions of 1m×1m×0.5m, and a 3m aircraft formation spacing; the distance between the three sets of aircraft connection points and the load connection points is measured using a laser rangefinder, and the average value is taken as the fixed distance) and the geometric relationships of the connection points (fixed coordinates of each aircraft connection point in the body coordinate system and fixed coordinates of the load connection point in the load's center of mass coordinate system). The system is logically structured in layers: aircraft translational dynamics, load translational dynamics, load rotational dynamics, and sling constraint relationships. The specific derivation and implementation of each part are as follows: When constructing the translational dynamics relationship of an aircraft, the translational dynamics relationship is used to describe the translational law under the combined action of thrust decomposition, gravity, and air resistance. It is derived through three steps: coordinate system transformation, acceleration calculation, and state update, forming a complete recursive relationship. The first step is thrust coordinate system transformation (dynamics relationship input preprocessing). Based on the roll, pitch, and yaw angles in the first real-time pose, a 3x3 rotation matrix is ​​generated in the order of Z-axis (yaw) - Y-axis (pitch) - X-axis (roll) (achieving attitude mapping). The thrust vector forward along the central axis in the aircraft body coordinate system is multiplied by this rotation matrix to obtain the lateral, longitudinal, and depth components of the thrust in the global coordinate system, completing the thrust coordinate system transformation and providing force input parameters for the translational dynamics relationship. The second step is translational acceleration calculation (core dynamics relationship). Based on the principle of force balance, the core relationship is constructed: Lateral translational acceleration = ( The thrust lateral component in the global coordinate system (the lateral component of the aircraft's gravity - the lateral component of air resistance) ÷ the aircraft's mass; the dynamic relationship between longitudinal translational acceleration and depth translational acceleration is constructed using the same logic, only replacing the force components in the corresponding directions. Among them, the gravity component is determined according to the direction of the global coordinate system (when the depth direction is opposite to the gravity direction, the gravity component is positive), and the air resistance component is proportional to the aircraft's velocity and opposite in direction, ensuring that the relationship conforms to physical laws; the third step is to update the translational state (iterative output of the dynamic relationship), constructing a recursive relationship with a control cycle (0.01 seconds) as the step size. The translational velocity at the current moment = the translational velocity at the previous moment + the translational acceleration × the control cycle; the position at the current moment = the position at the previous moment + the translational velocity at the current moment × the control cycle. The updated position data is checked in a closed loop with the position data of the first real-time pose to ensure the accuracy of the aircraft's translational dynamic relationship.

[0047] When constructing the load translational dynamics relationship, it is used to describe the load translational law under the combined action of cable tension, gravity, and air resistance. Through tension conversion, resultant force calculation, and state association derivation, a force transmission relationship is formed with the aircraft translational dynamics relationship. The first step is cable tension coordinate system transformation (dynamic relationship force input processing). For each cable, a coordinate transformation relationship is constructed. First, the fixed coordinates of the load connection point in the load centroid coordinate system are converted into global coordinate system coordinates (load connection point global coordinates) through the rotation matrix of the load optimization real-time pose. Simultaneously, the fixed coordinates of the aircraft connection point in the body coordinate system are converted into global coordinate system coordinates (aircraft connection point global coordinates) through the rotation matrix of the first real-time pose. The difference vector between the two is calculated and divided by the preset cable length (2.0m). The first step is to obtain the unit direction vector of the tension (from the load to the aircraft along the cable axis). Multiplying the tension magnitude by this vector yields the lateral, longitudinal, and depth tension components of each cable in the global coordinate system, providing the resultant force input for the load translational dynamics. The second step is to calculate the translational acceleration (core dynamics relationship). Based on the resultant force of all cables, the relationship is constructed: lateral translational acceleration of the load = (sum of lateral tension components of all cables - lateral component of load gravity - lateral component of load air resistance) ÷ load mass. The translational acceleration dynamics relationships in the longitudinal and depth directions are constructed using the same logic, only replacing the force components in the corresponding directions. The load gravity direction is opposite to the depth direction in the global coordinate system, and the air resistance component is proportional to the load's velocity and opposite in direction, ensuring that the load translational dynamics relationship conforms to the laws of rigid body motion.

[0048] When constructing the load rotation dynamics relationship, it is used to describe the load attitude change law under the action of torque. It is constructed in three steps: torque calculation, acceleration derivation, and attitude update. Together with the load translational dynamics relationship, it completely represents the rigid body motion of the load. The first step is to calculate the tension torque (torque input of the dynamics relationship). For each cable, a torque calculation relationship is constructed. First, the position vector of the global coordinates of the load connection point and the global coordinates of the load's center of mass (position coordinates in the real-time pose optimization of the load) is calculated (representing the spatial position of the connection point relative to the center of mass). This position vector is cross-multiplied with the cable tension vector to obtain the torque vector of a single cable about the load's center of mass (containing roll, pitch, and yaw components). The torque vectors of all cables are added together according to their components to obtain the total tension torque driving the load rotation, which serves as the core input of the load rotation dynamics relationship. The second step is to calculate the rotational acceleration (core dynamics relationship), which is combined with the balance relationship between the moment of inertia and the torque. The relationship is constructed as follows: Load roll acceleration = (Roll component of total tension torque - Air damping roll torque) ÷ Load roll moment of inertia. The dynamic relationship of rotational acceleration in the pitch and yaw directions is constructed using the same logic. The air damping roll torque is determined through ground dynamic calibration. The load is fixed on a controllable speed test bench and rotates uniformly at different roll angular velocities (5° / s, 10° / s, 15° / s). The corresponding damping torque is measured using a torque sensor, and the ratio of damping torque to angular velocity is calculated and averaged as the air damping roll coefficient. Finally, the air damping roll torque is calculated as: Air damping roll coefficient × Load roll angular velocity, with its direction opposite to the angular velocity. The moment of inertia is obtained through static calibration using a three-wire pendulum method. A three-wire pendulum device is constructed, and the load is placed horizontally at the center of the suspension plate. Parameters such as suspension length, suspension plate radius, and load mass are measured, and the simple harmonic oscillation period is recorded and averaged. This averaged value is then substituted into the moment of inertia formula: (Load mass × Gravitational acceleration × ... × ) ÷ (4 × The rotational inertia in the roll, pitch, and yaw directions is calculated using the length of the suspension line to ensure the accuracy of the relationship parameters. The third step is to update the rotational state (dynamic relationship iterative output). A recursive relationship is constructed with the control cycle (0.01 seconds) as the step size. The rotational angular velocity at the current moment = the rotational angular velocity at the previous moment + the rotational acceleration × the control cycle; the attitude angle at the current moment = the attitude angle at the previous moment + the rotational angular velocity at the current moment × the control cycle. The updated attitude angle data is then checked in a closed loop with the attitude angle data of the load optimization real-time pose to ensure the reliability of the load rotational dynamic relationship.

[0049] When constructing the sling constraint relationship, the sling constraint relationship is used to ensure the physical consistency of the model. Based on the characteristic of the fixed length of the rigid sling, it serves as a constraint condition for the motion of the associated aircraft and the load. First, the actual length calculation relationship of the sling is constructed, that is, the spatial distance between the global coordinates of the aircraft connection point and the global coordinates of the load connection point corresponding to each sling is calculated. The calculation process is the square root of the sum of the square of the difference in the lateral coordinates + the square of the difference in the longitudinal coordinates + the square of the difference in the depth coordinates. Then, the calculation result is set to be equal to the preset sling length (2.0m), forming a constraint equation, which forces the sling length in the model to remain unchanged, accurately restoring the mechanical characteristics of the rigid sling. The above-mentioned translational dynamics relationship of the aircraft, translational dynamics relationship of the load, rotational dynamics relationship of the load, and sling constraint relationship are integrated. Each part forms a closed loop through the relationship between force and coordinate, and finally a complete collaborative hoisting kinematics model is constructed to comprehensively characterize the dynamic behavior of the multi-aircraft collaborative hoisting system.

[0050] Step 4.2: Based on the collaborative crane kinematics model, with the control objectives of suppressing load swaying and tracking the preset formation trajectory, a collaborative crane kinematics model predictive controller is designed. In each control cycle, based on the current first real-time attitude of the aircraft and the optimized real-time attitude of the load, the optimization problem of minimizing the attitude tracking error between the aircraft and the load within the finite time domain is iteratively solved. The thrust and torque commands of each aircraft at the current moment are calculated to obtain the preliminary collaborative control commands. Specifically, based on the collaborative crane kinematics model established in Step 4.1, with the core control objectives of suppressing load swaying and tracking the preset formation trajectory, a model predictive controller is designed. Dynamic control is achieved through real-time optimization within each control cycle. The specific process is as follows: First, clarify... The dual core control objectives are: first, to suppress load sway by minimizing the fluctuation amplitude and rate of change of load roll angle, pitch angle, and yaw angle to maintain load attitude stability; and second, to track the preset formation trajectory by minimizing the deviations between the positions of each aircraft and the preset formation trajectory, and between the load position and the preset hoisting trajectory, to ensure formation flight accuracy and hoisting trajectory accuracy. Based on these objectives, a weighted sum objective function is constructed, specifically including three types of error terms. The aircraft trajectory tracking error term is calculated in each control cycle, taking into account the lateral, longitudinal, and depth positions of each aircraft and the aircraft trajectory tracking error term. The preset trajectory is determined before the multi-aircraft collaborative hoisting mission, taking into account the hoisting scenario requirements (avoiding equipment obstacles and maintaining a formation spacing of 3m). The planned aircraft formation reference flight path is generated by an offline path planning algorithm. It includes the lateral, longitudinal, and depth position reference values ​​of each aircraft in the global coordinate system within each control cycle. Within each control cycle, the squared differences between the lateral position and the corresponding lateral position reference value of each aircraft in the preset trajectory cycle, the squared differences between the longitudinal position and the corresponding longitudinal position reference value of each aircraft in the preset trajectory cycle, and the squared differences between the depth position and the corresponding depth position reference value of each aircraft in the preset trajectory cycle are calculated. These differences are then multiplied by preset weights (lateral 1.0, longitudinal 1.0, depth 1.2, with weight allocation based on higher priority for lifting accuracy in the depth direction) and summed to quantify the tracking deviation of the aircraft formation trajectory. The load trajectory tracking error term is calculated for each control cycle. Within each cycle, the squared differences between the load's lateral, longitudinal, and depth positions and the corresponding positions on the preset hoisting trajectory are calculated. These squared differences are then multiplied by preset weights (lateral 1.5, longitudinal 1.5, depth 2.0, with the load hoisting trajectory accuracy as the core indicator) and summed to quantify the tracking deviation of the load hoisting trajectory. The load swing suppression error term is maintained. Within each control cycle, the squared differences between the load's roll angle, pitch angle, yaw angle, and desired attitude angle (usually 0 degrees, i.e., horizontally stable attitude), as well as the squared changes in each attitude angle, are calculated. These squared differences are then multiplied by preset weights (attitude angle deviation weight 2.5, attitude angle change rate weight 1.8, with the weight allocation based on prioritizing attitude deviation suppression before reducing the swing rate) and summed to quantify the load swing amplitude.

[0051] The objective function is optimized to minimize the weighted sum of the three types of error terms, achieving synergistic optimization of the dual control objectives. To ensure the safe and stable operation of the system, the following constraints are set during the optimization process to form boundary limits on the control variables and system state. Within each control cycle, the actual length of each sling (calculated using the coordinates of the aircraft connection point and the load connection point) must be strictly equal to the preset sling length (2.0m), consistent with the constraint equations in the dynamic model, ensuring the effectiveness of the rigid sling assumption. The thrust of each aircraft is limited to 120N (minimum) to 180N (maximum). The minimum thrust must offset the self-weight + 10% redundancy, and the maximum thrust is 80% of the rated output of the power system (to avoid long-term full-load overload), ensuring both flight can be maintained and structural damage to the aircraft is prevented. The roll, pitch, and yaw moments of each aircraft are limited to ±20N·m. This range is determined by the torque carrying capacity limit of the quadcopter blade layout (exceeding this limit will cause the motor to stall), while also matching the ±30 degree attitude angle adjustment requirement, balancing attitude response speed and power system safety; the aircraft attitude angle is controlled between -0 degrees and 30 degrees, because exceeding this range will cause a sharp drop in the vertical component of thrust (e.g., when rolling 30 degrees, the vertical component is only 86.6% of the thrust), making it unable to stably bear the load; the load attitude angle is controlled between -15 degrees and 15 degrees, which is the safety threshold under rigid sling constraints, exceeding this range will cause a sudden change in the sling tension torque, leading to the risk of swinging and loss of control; the translational speed of each aircraft must be limited to a preset safety range (maximum speed 5m / s, minimum speed 0.5m / s, set according to the safety requirements of the hoisting scenario, avoiding high-speed inertial impact in indoor hoisting environments), and the translational speed of the load must be limited to a preset safety range (maximum speed 3m / s, minimum speed 0.3m / s, set according to the load weight and sling strength) to avoid excessive speed leading to excessive inertia.

[0052] Within each control cycle (0.01 seconds), the optimization problem is iteratively solved according to the following steps to generate preliminary cooperative control commands and realize rolling optimization control. First, the first real-time pose (position, attitude angle) of each aircraft and the optimized real-time pose (position, attitude angle) of the load at the current moment are collected and used as the initial input state of the controller to ensure that the optimization solution is based on the current real operating state of the system. Second, based on the cooperative crane kinematics model, the system state in the future finite time domain (such as the next 10 control cycles, i.e., 0.1 seconds) is predicted, including the position, attitude angle of each aircraft and the position, attitude angle of the load at each prediction moment. The prediction process involves iterative calculation of the model's state equations. These state equations are the core relationships of the collaborative crane kinematics model constructed in step 4.1, including the translational dynamics of the aircraft, the translational dynamics of the load, and the rotational dynamics of the load. The state at each prediction moment is determined by the state at the previous prediction moment, along with the thrust and torque commands. The recursive calculation of the aforementioned state equations allows for the prediction of the system's future dynamic behavior. The third step involves using the thrust and torque commands of each aircraft within a finite future time domain as optimization variables, minimizing the objective function as the optimization objective, and substituting the aforementioned constraints to construct a constrained nonlinear optimization problem, clarifying the optimization direction and boundary limitations. The fourth step uses numerical optimization algorithms (such as the interior point method) to iteratively solve the optimization problem, first setting the thrust and torque commands... The initial feasible values ​​(all within the preset safety constraints) are determined. In each iteration, the predicted future system state (predicted pose of the aircraft and load) corresponding to the current optimization variable within the finite time domain is substituted into the objective function. The objective function value is calculated according to the weighted sum rule of three types of error terms (i.e., the aircraft trajectory tracking error, load trajectory tracking error, and load sway suppression error are multiplied by their respective preset weights and then summed). Then, based on the negative gradient direction of the objective function (which guides the variable to adjust in the direction of decreasing objective function), an adaptive step size strategy is used to update the variable. First, the initial step size is set according to the variable magnitude (e.g., thrust step size 0.5N, torque step size 0.1N·m). The initial update value is calculated by adding the current variable value to the negative gradient and then checking against the constraint conditions. The thrust needs to be between 120 and 180N, and the torque needs to be between 120 and 180N. If the constraint of 20 N·m is met, the updated value is retained; otherwise, the step size is reduced to half of the original, the updated value is recalculated and verified, until the updated thrust and torque commands meet all constraints; the optimal solution is continuously approximated through multiple rounds of iteration until the difference between the objective functions of two consecutive rounds is less than the preset convergence threshold (e.g., ...). Finally, the thrust and torque command sequence for the next 10 control cycles that satisfies all constraints and minimizes the objective function is obtained. In the fifth step, the control command is output, and the thrust and torque command corresponding to the first control cycle in the optimal command sequence is taken as the initial coordinated control command for each aircraft at the current moment. The above steps are repeated for each control cycle to continuously update the control command and ensure that the controller can dynamically adapt to changes in the system state and achieve real-time optimized control.

[0053] In this embodiment, the collaborative hoisting kinematics model comprehensively considers the aircraft's translation, load rotation, and sling constraints, characterizing the dynamic characteristics of the multi-aircraft collaborative hoisting system. This provides reliable model support for control strategy design and avoids control deviations caused by model simplification. The model predictive controller achieves both preset formation trajectory tracking and load sway suppression through rolling optimization within a finite time domain. This ensures the accuracy of multi-aircraft formation flight while effectively reducing load sway amplitude, improving the stability and accuracy of collaborative hoisting. The controller performs dynamic optimization based on the first real-time pose acquired in real time and the optimized real-time pose of the load, enabling it to quickly adapt to state changes and external disturbances during motion, enhancing the robustness of the control strategy. Strict constraint settings ensure that the motion of the aircraft and load remains within a safe range, effectively avoiding dangerous situations such as overload, instability, and excessive sway, thus improving the safety and reliability of multi-aircraft collaborative hoisting. Control commands are generated through real-time iterative solutions, keeping synchronized with the control cycle of the aircraft's flight control system. This balances real-time control with optimization effects, meeting the dynamic control requirements of multi-aircraft collaborative hoisting.

[0054] In a preferred embodiment of the present invention, step 5 includes: Step 5.1: Based on the load position coordinates in the real-time pose of the load optimization, update the spatial distribution of the strong signal manifold and the weak signal homotopy group to obtain the updated topology of the strong signal manifold and the weak signal homotopy group. Specifically, the positioning reliability of multi-machine collaborative hoisting depends on the spatial topology of the laser signal coverage. The real-time movement of the load will inevitably lead to changes in the relative position of the signal coverage. Therefore, it is necessary to dynamically adjust the signal topology based on the real-time pose of the load. First, the topology is defined as follows: a strong signal manifold is a spatial region with dense laser positioning signal coverage and positioning accuracy meeting preset requirements (positioning error ≤ 0.05 meters); a weak signal homotopy group is a transitional region outside the strong signal manifold where laser signal intensity attenuates but is still detectable (positioning error ≤ 0.2 meters). Next, real-time input data is acquired, collecting the global coordinates (X-axis, Y-axis, Z-axis coordinates) of the load centroid in the real-time pose optimization. These coordinates are measured and filtered in real-time by the laser positioning system to ensure data accuracy. Then, the updated boundary is calculated, using the current coordinates of the load centroid as a reference, combined with preset signal coverage parameters (1.5 meters half-width coverage in the X / Y axis direction and 0.8 meters half-width coverage in the Z-axis direction). These parameters are determined through a complete performance calibration-scenario testing process, specifically as follows: an industrial-grade laser positioning module with a rated power of 5 watts is selected, and calibration is first performed in an unobstructed darkroom environment. The basic performance was determined, confirming the theoretical coverage radius with a positioning error ≤ 0.05 meters. Subsequently, in a test site simulating the signal propagation characteristics of mountain valleys, typical terrain obstructions (such as rock walls and vegetation) were used as references to determine the actual coverage half-width through actual measurements, ensuring that the parameters both met the module's performance and avoided the influence of actual terrain obstructions. Based on this, the updated boundary of the strong signal manifold was obtained, i.e., the X-axis direction is ±1.5 meters of the load centroid X coordinate, and the Y-axis and Z-axis boundaries are similarly determined. Then, based on the strong signal manifold boundary, a preset transition width was extended outward (0.5 meters in the X / Y axis direction and 0.3 meters in the Z-axis direction). The determination of this width also relied on actual measurement data. Starting from the strong signal boundary, the receiving module was moved horizontally at 0.1-meter intervals, and the changes in positioning error were recorded. It was found that for every 0.5-meter outward movement, the error steadily increased by 0.1 meters; for every 0.3-meter vertical movement, the error also increased by 0.1 meters. Based on the requirement that the positioning error of the weak signal homotopy group should be ≤0.2 meters, after expanding outward by 0.5 meters (horizontally) and 0.3 meters (vertically) from the strong signal boundary (error 0.05 meters), the error reaches exactly 0.15 meters, which is within the safe threshold. Considering the 1m x 1m bottom surface size of the load, this width avoids the transition area being too narrow, causing the signal to be triggered only when the load edge touches the warning boundary, thus reserving sufficient reaction space. Based on this, the boundary of the weak signal homotopy group is obtained, that is, the X-axis direction is the X-boundary of the strong signal manifold ±0.5 meters, and the Y-axis and Z-axis boundaries follow the same logic. Finally, the updated topology is determined, that is, the spatial range and relative position of the two types of regions are clarified according to the above boundary coordinates, ensuring that the topology dynamically adjusts synchronously with the load's centroid, always conforming to the load's lifting trajectory.

[0055] Step 5.2: Based on the updated topology of the strong signal manifold and the weak signal homotopy group, construct the control Lyapunov potential function. The input data of the control Lyapunov potential function is the coordinates of the load centroid in the global coordinate system. The function value of the control Lyapunov potential function monotonically increases as the shortest spatial distance between the load centroid coordinates and the weak signal homotopy group region decreases. Specifically, this includes: Step 5.1: The updated topology provides the basis for spatial relationship quantification. The control Lyapunov potential function is the core bridge connecting the topology and control decisions. By converting the load spatial location into a computable function value, the quantitative determination of the safety state is achieved. First, the function input is clarified. The core logic involves inputting real-time global coordinates of the load centroid (X, Y, and Z). The core logic quantifies the spatial relationship between the load and the weak signal homotopy group. The function value must monotonically increase as the shortest distance between the load and the weak signal homotopy group decreases; the closer the distance (closer to the signal edge), the larger the function value, thus intuitively reflecting the safety risk level. Next, the shortest spatial distance is calculated. First, the boundary range of the weak signal homotopy group updated in step 5.1 is extracted (X1 to X2, Y1 to Y2, and Z1 to Z2). Then, the distance from the load centroid coordinates (X0, Y0, Z0) to the boundary in each direction is calculated. If the centroid is within the boundary range... If the distance in that direction is 0, then the distance to the nearest boundary is taken; otherwise, the distance to the nearest boundary is taken. Finally, the shortest spatial distance is calculated using the spatial distance formula to characterize the closest positional relationship between the load and the weak signal homotopy group. Finally, a potential function is constructed, using a monotonically increasing function form to ensure the effectiveness of risk quantization. The baseline constant is set to 1.0. The core purpose is to avoid quantization failure due to an excessively small potential function value when the load is in a safe region (the shortest spatial distance to the weak signal homotopy group is relatively far). 1.0 serves as the basic anchor point for the function value, ensuring that the function values ​​corresponding to different positions within the safe region maintain reasonable differentiation, clearly reflecting the positional differences of the load within the safe range. If the value is less than 1.0, the load in... At the far end of the safe zone, the function value may approach zero, making effective quantification within the safe zone impossible. The scaling factor is set to 5.0 to match the transition width of the weak signal homotopy group, ensuring that the potential function value increases smoothly with the risk level when the load moves within the weak signal zone. A value of 5.0 allows for a moderate range of function value changes, preventing both missed risk detection due to slow changes and control oscillations due to sudden increases. If the scaling factor is too large, a short-distance movement will cause a surge in the function value, resulting in excessive sensitivity. If the value is too small, the function value changes slowly, making it impossible to effectively distinguish different risk levels. The formula for calculating the potential function value is: scaling factor ÷ shortest spatial distance + reference constant.

[0056] Step 5.3: Select the equipotential surface where the control Lyapunov potential function value equals a preset threshold as the interface, and divide the set of load centroid coordinates into two regions. The region where the potential function value is lower than the preset threshold is defined as the central operating region, and the region where the potential function value is greater than or equal to the preset threshold is defined as the boundary warning region. Specifically, this includes: conducting the calibration of the preset threshold of the control Lyapunov potential function. The calibration needs to be performed in conjunction with the actual operation scenario of multi-machine hoisting in mountainous canyons and the system control accuracy indicators; in the test site simulating mountainous canyon terrain, deploy laser positioning base stations, multiple aircraft, and standard hoisting loads according to actual operation requirements, and plan a complete hoisting trajectory covering the central operating region and the weak signal homotope group; control the load to move gradually from the central operating region to the weak signal homotope group along the trajectory, and simultaneously collect three types of data: load optimization... The system includes the three-dimensional coordinates of the centroid in the pose, the corresponding control Lyapunov potential function value, and the real-time pose tracking error of the load (including position error and attitude error). The system control accuracy threshold is set to position error ≤ 0.1 meters and attitude error ≤ 0.5 degrees. When the load moves to a certain position, if the pose tracking error reaches the critical value of this accuracy threshold for the first time (i.e., the position error is close to 0.1 meters or the attitude error is close to 0.5 degrees), the control Lyapunov potential function value at this time is recorded as a single set of critical potential function values. This test process is repeated 50 times. After removing abnormal data with deviations exceeding 30% of the average of all data, the remaining 40 sets of valid critical potential function values ​​are calculated using an arithmetic mean. Specifically, the 40 sets of critical potential function values ​​are summed sequentially to obtain a total, and then this sum is divided by 40. The final result is the preset threshold.

[0057] An equipotential surface whose control Lyapunov potential function value equals a preset threshold is selected as the region boundary. This boundary is a spatially closed surface that dynamically matches the topology of the strong signal manifold and the weak signal homotopy group. When the control Lyapunov potential function value corresponding to the load centroid is lower than the preset threshold, this region is defined as the central operating region. Within this region, the laser positioning signal coverage is stable, and the load pose calculation accuracy meets the requirements of cooperative control, requiring no additional attitude compensation. When the control Lyapunov potential function value corresponding to the load centroid is greater than or equal to the preset threshold, this region is defined as the boundary warning region. Within this region, the laser positioning signal is attenuated or blocked, and the load pose calculation accuracy may decrease. Dynamic attitude compensation is required to correct the control commands to maintain hoisting stability.

[0058] Step 5.4: Real-time determination of the region where the load centroid coordinates are located in the real-time pose optimization. When the load centroid coordinates are located within the boundary warning area, dynamic attitude compensation commands are calculated based on the positional relationship between the load centroid coordinates and the interface. The dynamic attitude compensation commands are then vector-superimposed with the preliminary cooperative control commands to generate the final control commands, which are then sent to the corresponding aircraft. Specifically, this includes: First, calculating the shortest spatial distance from the load centroid to the interface. The interface is an equipotential surface where the control Lyapunov potential function value equals a preset threshold. Its spatial equation is derived based on the topological boundary of the weak signal homotopy group and the definition of the potential function. Let the boundary range of the weak signal homotopy group in the global coordinate system be [x1, x2] horizontally, [y1, y2] vertically, and [z1, z2] in depth. The real-time coordinates of the load centroid are (x0, y0, z0), the potential function scaling factor is (value 5.0), the reference constant is (value 1.0), and the preset threshold is λ. Then, the spatial equation of the interface is... = proportionality coefficient / (λ - reference constant), where, Let x0 be the distance from the load centroid to the boundary of the weak signal homotopy group in the transverse direction, when x0∈[x1, x2]. =0, otherwise =min(|x0-x1|,|x0-x2|); Let be the distance from the load centroid in the longitudinal direction to the boundary of the weak signal homotopy group, when y0∈[y1, y2] =0, otherwise =min(|y0-y1|,|y0-y2|); Let z0 be the distance from the centroid of the load along the depth direction to the boundary of the homotopy group of the weak signal, when z0∈[z1, z2]. =0, otherwise =min(|z0-z1|,|z0-z2|); Based on the above interface spatial equation and the real-time coordinates of the load centroid, the method of finding the shortest distance from a spatial point to the surface is adopted. First, the coordinates of the perpendicular foot point on the interface, which is perpendicular to the interface and connected to the load centroid, are determined by the topological equation of the interface. Then, the spatial distance between the coordinates of the load centroid and the coordinates of the perpendicular foot point is calculated. The specific calculation process is as follows: the horizontal coordinate, vertical coordinate, and depth coordinate of the load centroid are subtracted from the corresponding coordinates of the perpendicular foot point to obtain the distance components in three directions. The distance components in each direction are squared, and the three squared results are added together and the square root is taken to finally obtain the shortest spatial distance from the load centroid to the interface.

[0059] The second step is to calculate the difference between the potential function value and the preset threshold. This is done by subtracting the preset threshold from the control Lyapunov potential function value calculated at the current moment. This difference directly reflects the risk level of the load's centroid within the boundary warning area. The larger the difference, the closer the load is to the core region of the weak signal homotopy group, the worse the signal quality, and the greater the required attitude compensation. The third step is to calculate the attitude angle compensation (in radians). The attitude angle compensation includes roll angle compensation, pitch angle compensation, and yaw angle compensation, all calculated uniformly. Logical calculation; taking roll angle compensation as an example, determine the roll angle compensation ratio coefficient (unit: radians / (meter·unit difference)). Through ground experiment calibration, within the boundary warning area, when the difference between the potential function value and the preset threshold is 1.0 and the shortest distance from the load centroid to the interface is 0.1 meters, the roll angle compensation needs to reach 0.3° (approximately 0.00524 radians) to offset the pose error. Based on this, the roll angle compensation ratio coefficient is calculated. First, calculate the distance correction coefficient: distance correction coefficient = preset reference distance ÷ The shortest spatial distance from the load's center of mass to the interface is calculated, with a preset reference distance of 0.1 meters. The roll angle compensation (in radians) is calculated as: roll angle compensation ratio coefficient × difference between potential function difference and preset threshold × distance correction coefficient. Following the same logic, the pitch angle compensation and yaw angle compensation (both in radians) are calculated separately. These three components together constitute the complete attitude angle compensation, ensuring targeted compensation for all three attitude dimensions of the load. The fourth step is to convert the attitude angle compensation into torque compensation (in Newton-meters). The calculation of the quantities is based on the attitude angle compensation, load moment of inertia, and control cycle, following the basic relationships of rotational dynamics. Taking the rolling torque compensation as an example, the load rolling moment of inertia is obtained by static calibration using the three-line pendulum method. Assuming the calibration value is 0.5 kg·m², the rolling angular velocity compensation (radians / second) = rolling angle compensation (radians) / control cycle (0.01 seconds). According to the rotational dynamics formula, torque = moment of inertia × angular acceleration, and the rolling torque compensation = load rolling moment of inertia × rolling angular velocity compensation, which is 0.Multiplying 5 kg·m² by the roll velocity compensation (radians / second) yields the roll torque compensation (N·m). Following the same logic, the pitch torque compensation and yaw torque compensation are calculated separately. These three together constitute the dynamic attitude compensation command, ensuring that the compensation command matches the control dimension of the aircraft's flight control system. In the fifth step, the dynamic attitude compensation command is vector-superimposed with the preliminary coordinated control command. The superposition process strictly follows the principle of corresponding torque component superposition. The final roll torque command = roll torque in the preliminary coordinated control command + roll torque compensation; the final pitch torque command = pitch torque in the preliminary coordinated control command + pitch torque compensation; the final yaw torque command = yaw torque in the preliminary coordinated control command + yaw torque compensation. The thrust command remains unchanged because its core function is to maintain the lift and translational state of the aircraft; attitude compensation does not require adjustment of the thrust parameters. After superposition, the final control command undergoes boundary verification to ensure that the final roll, pitch, and yaw moment commands are all within a preset constraint range of ±20 N·m. If the moment command in any dimension exceeds the constraint range, the critical value (maximum or minimum value) of that constraint range is taken as the final moment command for that dimension to avoid overload of the power system. After the verification passes, the final control command is sent to the corresponding aircraft flight control system in real time through a stable data transmission link to achieve dynamic compensation of the load attitude, ensuring that even within the boundary warning area, the load can still maintain a stable lifting state, avoiding the risk of load swaying or collision due to signal quality degradation.

[0060] This embodiment identifies the spatial position of the load by updating the signal topology and constructing a potential function in real time, preventing the load from exceeding the reliable range of laser positioning and reducing trajectory deviations caused by signal instability. By dividing the central operating area and the boundary warning area, the attitude compensation mechanism is triggered in advance, effectively preventing the load from touching the lifting boundary or obstacles and improving the safety of collaborative lifting. The dynamic attitude compensation command is accurately calculated based on the positional relationship between the load and the interface. After being vector-superimposed with the initial collaborative control command, it can achieve smooth and rapid attitude correction, ensuring the stability of the load attitude. Combining the characteristics of the laser positioning signal with the real-time dynamic adjustment control strategy of the load's posture, the system can adapt to positional changes during lifting, enhancing its resistance to external interference. By quantifying the positional relationship and precise compensation control through the potential function, the system ensures that the load always moves within the preset safety range, maintaining the accuracy of formation flight and the accuracy of the load lifting trajectory.

[0061] like Figure 2 As shown, embodiments of the present invention also provide a ground-based laser positioning multi-machine hoisting collaborative control system for mountainous areas, including: The calculation module is used to calculate the first real-time pose of the aircraft based on the laser positioning signal emitted by the laser positioning base station; based on the image data of the three optical markers, it calculates the three-dimensional spatial coordinates of the three optical markers to define the spatial geometric surface patch, and decomposes the surface patch into a strong signal manifold and a weak signal homotopy group; The decomposition module is used to connect the two optical markers at the center point of the upper edge of the front face of the load structure as the spatial reference origin, and form a first spatial vector and a second spatial vector by connecting the two optical markers at the left and right vertices of the lower edge of the rear face. The module is used to construct a local coordinate system for the load based on the spatial reference origin and the first and second spatial vectors, and to calculate the initial pose of the load; the initial pose of the load is optimized and calibrated by fusing ranging data with the first real-time pose to obtain the optimized real-time pose of the load. The calibration module is used to establish a collaborative hoisting kinematics model based on the first real-time pose and the load-optimized real-time pose, so as to calculate the initial collaborative control commands for each aircraft. The control module is used to determine the line-of-sight visibility of the load to each laser positioning base station based on the load's real-time pose optimization. Based on the topology of the strong-signal manifold and the weak-signal homotopy group, it constructs a control Lyapunov potential function and an equipotential surface of the control Lyapunov potential function, dividing the operating space where the load's centroid is located into a central operating region and a boundary warning region. When the load's centroid coordinates are within the boundary warning region, it generates dynamic attitude compensation commands to correct the initial cooperative control commands and outputs the final control commands to each aircraft. It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0062] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for ground laser positioning mountainous area multi-machine hoisting cooperative control, characterized in that, The method comprises: Step 1, laser positioning base station emits laser positioning signal, and the first real-time pose of the aircraft is solved; according to the image data of the three optical markers, the spatial three-dimensional coordinates of the three optical markers are solved to define the spatial geometric curved surface sheet, and the curved surface sheet is decomposed into strong signal manifold and weak signal homotopy group; Step 2, taking the optical marker of the center point of the upper edge of the front end face of the load structure as the spatial reference origin, connecting the two optical markers of the left top point and the right top point of the lower edge of the rear end face to form the first space vector and the second space vector; Step 3, constructing the local coordinate system of the load according to the spatial reference origin, the first and second space vectors, and solving the initial pose of the load; the initial pose of the load is optimized and calibrated by fusing the ranging data and the first real-time pose to obtain the optimized real-time pose of the load; Step 4, a cooperative hoisting motion mechanics model is established according to the first real-time pose and the optimized real-time pose of the load to calculate the preliminary cooperative control command of each aircraft; Step 5, the line-of-sight visibility of the load to each laser positioning base station is determined according to the optimized real-time pose of the load; based on the topological structure of the strong signal manifold and the weak signal homotopy group, a control Lyapunov potential function is constructed, and the equipotential surface of the control Lyapunov potential function divides the operation space where the center of mass of the load is located into a central operation area and a boundary warning area; when the center of mass of the load is located in the boundary warning area, a dynamic attitude compensation command is generated to correct the preliminary cooperative control command, and a final control command is output to each aircraft.

2. The method according to claim 1, wherein, Before step 1, the laser positioning base station is arranged at a ground position in the canyon area, and three non-coplanar optical markers are arranged on the load structure, which are located at the center point of the upper edge of the front end face, the left top point of the lower edge of the rear end face, and the right top point of the lower edge of the rear end face; an airborne positioning and perception unit is configured for each aircraft to obtain aircraft position data, image data containing three optical markers, and aircraft and load ranging data.

3. The method according to claim 2, wherein, The step 1 comprises: The laser positioning signal obtained by the airborne positioning and perception unit is processed, and the first real-time pose of the aircraft in the global coordinate system is solved by solving a stochastic differential equation containing a measurement noise term; the image data containing three optical markers obtained by the airborne positioning and perception unit is processed, the sub-pixel level image coordinates of each optical marker are extracted, and the spatial three-dimensional coordinates of each optical marker in the camera coordinate system are solved according to the image coordinates of the optical marker and the intrinsic parameters of the vision sensor; The spatial three-dimensional coordinates of the three optical markers are taken as the vertices to construct the spatial geometric curved surface sheet, for each point on the spatial geometric curved surface sheet, the theoretical line-of-sight vector between each point and each laser positioning base station is calculated, and according to the first real-time pose, it is judged whether the theoretical line-of-sight vector is blocked by the terrain or the load itself, and the signal-to-noise ratio level of the received signal is evaluated to obtain the blocking judgment result of the theoretical line-of-sight vector and the evaluation result of the signal-to-noise ratio level; According to the occlusion judgment result of the theoretical line-of-sight vector and the evaluation result of the signal-to-noise ratio level, the continuous point set region with more than two unoccluded lines of sight and the average signal-to-noise ratio higher than the preset threshold on the spatial geometric curved surface sheet is divided into a strong signal manifold; and the point set region with occluded line of sight or the average signal-to-noise ratio lower than the preset threshold is divided into a weak signal homotopy group.

4. The method according to claim 3, wherein, The spatial geometric curved surface sheet is constructed with the spatial three-dimensional coordinates of the three optical marker points as vertices, for each point on the spatial geometric curved surface sheet, the theoretical line-of-sight vector between each point and each laser positioning base station is calculated, and whether the theoretical line-of-sight vector is occluded by the terrain or the load itself is judged according to the first real-time pose, and the signal-to-noise ratio level of the received signal is evaluated, to obtain the occlusion judgment result of the theoretical line-of-sight vector and the evaluation result of the signal-to-noise ratio level, including: The spatial geometric curved surface sheet is constructed with the spatial three-dimensional coordinates of the three optical marker points as vertices, and a two-dimensional parameterized grid is established in the spatial geometric curved surface sheet, with the first optical marker point as the starting point, the second optical marker point as the first edge direction, and the third optical marker point as the second edge direction; on the two-dimensional parameterized grid, the grid points with integer values of horizontal and vertical coordinates are selected as sampling points at an interval of 0.05 meters to obtain twenty-one sampling points; For each sampling point, the geometric straight line vector from the sampling point to each laser positioning base station is calculated as the theoretical line-of-sight vector; the real-time three-dimensional envelope range of the aircraft in space is determined according to the first real-time pose of the aircraft and the spatial three-dimensional coordinates of the three optical marker points, and the spatial range occupied by the load is determined based on the spatial positional relationship of the three optical marker points; It is judged whether each theoretical line-of-sight vector intersects with the real-time three-dimensional envelope range of the aircraft or the spatial range of the load, if the theoretical line-of-sight vector intersects with the spatial range of the aircraft or the load, it is determined that the corresponding theoretical line-of-sight vector is occluded; otherwise, it is determined that the corresponding theoretical line-of-sight vector is not occluded, to obtain the occlusion judgment result; For the unoccluded theoretical line-of-sight vector, the signal-to-noise ratio level of the received signal on the corresponding path is evaluated based on the transmission power of the laser positioning signal, the atmospheric attenuation characteristics of the signal on the transmission path, and the noise coefficient of the receiver, to obtain the signal-to-noise ratio evaluation result; The occlusion judgment results and the signal-to-noise ratio evaluation results of all sampling points for all laser positioning base stations are summarized as the occlusion judgment result of the theoretical line-of-sight vector and the evaluation result of the signal-to-noise ratio level.

5. The method according to claim 4, wherein, The step 2 includes: The coordinates of the optical marker point located at the center of the upper edge of the front end face of the load structure are selected from the spatial three-dimensional coordinates of the three optical marker points, and defined as the coordinates of the spatial reference origin; The difference between the spatial three-dimensional coordinates of the optical marker point located at the left top point of the lower edge of the rear end face and the coordinates of the spatial reference origin is calculated to obtain a first spatial vector; and the difference between the spatial three-dimensional coordinates of the optical marker point located at the right top point of the lower edge of the rear end face and the coordinates of the spatial reference origin is calculated to obtain a second spatial vector.

6. The method according to claim 5, wherein, The step 3 includes: unitize the first spatial vector as the first coordinate axis of the load local coordinate system, and calculate the dot product of the second spatial vector and the first coordinate axis; subtract the product of the dot product and the first coordinate axis from the second spatial vector to obtain an intermediate vector orthogonal to the first coordinate axis, and unitize the intermediate vector as the second coordinate axis of the load local coordinate system; calculate the cross product of the first coordinate axis and the second coordinate axis to obtain the third coordinate axis of the load local coordinate system; and define the load local coordinate system with the spatial reference origin, the first coordinate axis, the second coordinate axis and the third coordinate axis; calculate the initial position and the initial attitude of the load in the global coordinate system according to the rotation matrix and the translation vector between the load local coordinate system and the global coordinate system, and obtain the initial pose of the load; obtain the ranging data between the aircraft and the load measured by the onboard positioning and perception unit, fuse the ranging data with the first real-time pose, iteratively optimize and calibrate the initial pose of the load by the recursive state estimation algorithm, and obtain the optimized real-time pose of the load.

7. The method according to claim 6, wherein, obtain the ranging data between the aircraft and the load measured by the onboard positioning and perception unit, fuse the ranging data with the first real-time pose, iteratively optimize and calibrate the initial pose of the load by the recursive state estimation algorithm, and obtain the optimized real-time pose of the load, including: receive the initial pose of the load as the initial state vector of the recursive state estimation algorithm; receive the first real-time pose as the known input condition of the recursive state estimation algorithm; and receive the ranging data between the aircraft and the load measured by the onboard positioning and perception unit as the observation of the recursive state estimation algorithm; in the prediction step of the recursive state estimation algorithm, predict the load state vector and the state covariance matrix at the current time from the state vector at the previous time based on the motion of the aircraft described by the first real-time pose and the motion state of the load; in the update step of the recursive state estimation algorithm, calculate the residual between the predicted value and the actual measured value of the relative ranging observation at the current time, and calculate the Kalman gain based on the state covariance matrix and the observation noise covariance matrix; use the Kalman gain to correct and update the load state vector and the state covariance matrix obtained in the prediction step to obtain the corrected load state vector and the state covariance matrix at the current time, and output the position and attitude information contained in the corrected load state vector as the optimized real-time pose of the load.

8. The method according to claim 7, wherein, The step 4 includes: describe the motion state of each aircraft with the first real-time pose, describe the rigid body motion state of the load with the optimized real-time pose of the load, and establish a cooperative hoisting motion mechanics model containing aircraft translation, load rotation and sling constraint based on the sling length and the connection point geometry; based on the cooperative hoisting motion mechanics model, design a cooperative hoisting motion mechanics model predictive controller with the control objectives of suppressing load swing and tracking the preset formation trajectory, and in each control period, iteratively solve the optimization problem of minimizing the pose tracking error of the aircraft and the load in the future finite time domain according to the first real-time pose of the current aircraft and the optimized real-time pose of the load, calculate the thrust and moment command of each aircraft at the current time, and obtain the preliminary cooperative control command.

9. The method according to claim 8, wherein, The step 5 comprises: According to the load position coordinates in the load optimization real-time pose, the spatial distribution of the strong signal manifold and the weak signal homotopy group is updated, and the topology structure of the updated strong signal manifold and the weak signal homotopy group is obtained; Based on the topology structure of the updated strong signal manifold and the weak signal homotopy group, a control Lyapunov potential function is constructed, the input data of the control Lyapunov potential function is the coordinates of the load centroid in the global coordinate system, and the function value of the control Lyapunov potential function is monotonically increasing with the decrease of the shortest spatial distance between the load centroid coordinates and the weak signal homotopy group region; An equipotential surface with the control Lyapunov potential function value equal to a preset threshold is selected as a dividing surface, and the set of load centroid coordinates is divided into two regions, wherein the region with a potential function value lower than the preset threshold is defined as a central operation region, and the region with a potential function value greater than or equal to the preset threshold is defined as a boundary warning region; The region where the load centroid coordinates in the load optimization real-time pose is located is determined in real time, and when the load centroid coordinates are located in the boundary warning region, a dynamic attitude compensation instruction is calculated according to the positional relationship between the load centroid coordinates and the dividing surface, the dynamic attitude compensation instruction is vector superimposed with the preliminary cooperative control instruction to generate a final control instruction, and the final control instruction is sent to the corresponding aircraft.

10. A ground laser positioning mountainous area multi-machine hoisting cooperative control system, the system implements the method according to any one of claims 1 to 9, characterized in that, Comprise: The solving module is used for solving the laser positioning signal transmitted by the laser positioning base station, solving the first real-time pose of the aircraft, and solving the three-dimensional spatial coordinates of the three optical markers according to the image data of the three optical markers to define a spatial geometric surface sheet and decompose the surface sheet into a strong signal manifold and a weak signal homotopy group; The decomposition module is used for taking the optical marker of the center point of the upper edge of the front end face of the load structure as a spatial reference origin, connecting the two optical markers of the left top point and the right top point of the lower edge of the rear end face to form a first spatial vector and a second spatial vector; The construction module is used for constructing a load local coordinate system according to the spatial reference origin, the first and second spatial vectors, and solving the load initial pose; and the load initial pose is optimized and calibrated by fusing the ranging data and the first real-time pose to obtain the load optimization real-time pose; The calibration module is used for establishing a cooperative hoisting motion mechanics model according to the first real-time pose and the load optimization real-time pose to calculate the preliminary cooperative control instruction of each aircraft; The control module is used for determining the line-of-sight visibility of the load to each laser positioning base station according to the load optimization real-time pose; based on the topology structure of the strong signal manifold and the weak signal homotopy group, a control Lyapunov potential function is constructed, the equipotential surface of the control Lyapunov potential function divides the operation space of the load centroid into a central operation region and a boundary warning region; when the load centroid coordinates are located in the boundary warning region, a dynamic attitude compensation instruction is generated to correct the preliminary cooperative control instruction, and a final control instruction is output to each aircraft.

Images