A polar unmanned aerial vehicle navigation method and system based on image assisted positioning

Abstract

This invention provides a polar unmanned aerial vehicle (UAV) navigation method and system based on image-assisted localization, relating to the field of computer vision. The navigation method includes: acquiring the coarse three-dimensional position and coarse attitude quaternions of the polar UAV; acquiring real polar starry sky images and constructing a continuous optical residual strain field; calculating the divergence and curl of the optical residual strain field, identifying and proposing refraction distortion regions and auroral interference regions to obtain an initial starry sky region; extracting multidimensional features within the initial starry sky region, iteratively correcting the effective starry sky region boundary based on Mahalanobis distance to obtain candidate starry sky regions; extracting the actual motion optical flow field and theoretical optical flow field within the candidate starry sky regions, verifying them in conjunction with the three-dimensional translational velocity vector calculated by a mobile local area network, eliminating false starry sky regions, and obtaining the optimal starry sky region; and completing the final precise positioning within the optimal starry sky region.

Description

Technical Field

[0001] This invention relates to the field of computer vision, and more particularly to a polar unmanned aerial vehicle (UAV) navigation method and system based on image-assisted localization. Background Technology

[0002] The polar regions possess extreme natural environmental characteristics, including high latitude, strong geomagnetic interference, alternating periods of polar day and night, and harsh weather conditions. For unmanned aerial vehicles (UAVs) performing scientific expeditions, glacier mapping, environmental observation, and communication relay missions in the polar regions, a highly reliable navigation and positioning system is fundamental to ensuring their long-term safe operation. In traditional high-latitude navigation applications, Global Navigation Satellite Systems (GNSS) are often affected by ionospheric disturbances, auroral activity, and poor satellite geometry, leading to signal attenuation, positioning jumps, and even loss of lock-on. Furthermore, due to their proximity to the geomagnetic poles, conventional magnetic compasses are completely ineffective in the polar environment, failing to provide accurate heading references. Therefore, polar UAVs typically rely heavily on Inertial Navigation Systems (INS) for continuous positioning calculations. However, INS errors accumulate over time, necessitating continuous correction using external high-precision sensors during long-endurance missions.

[0003] To suppress error divergence in inertial navigation, existing technologies often incorporate terrestrial radio local area networks (LANs) as an auxiliary means. However, traditional LAN-assisted positioning heavily relies on fixed base station infrastructure with known prior locations. In polar scenarios, large-scale ice surface drift is frequent, and research teams are often in a dynamic relocation state, making it impossible to provide absolutely stationary fixed base stations. If mobile base stations mounted on snowmobiles or research vessels are forcibly introduced into the traditional solution architecture, the lack of joint estimation capability for dynamic network reference frames in conventional algorithms, treating them as isolated ranging nodes, will cause severe coordinate system drift.

[0004] To address the technical bottleneck of unobservable absolute heading in dynamic network architectures, stellar image-based assisted localization has become a preferred approach for obtaining absolute attitude references. However, the reliability of traditional stellar image localization methods faces severe challenges in the complex atmospheric and optical physical environment of the polar regions. The polar night sky is often accompanied by high-intensity auroral activity, and the low-frequency diffuse light of the aurora is easily misextracted as bright features by traditional image processing algorithms, leading to a large number of false star points being mixed into the matching process. Simultaneously, the extremely cold climate easily forms strong temperature inversion layers near the surface, causing severe low-altitude atmospheric refraction distortion, resulting in abnormal positional shifts of stars near the horizon, causing traditional blind star map matching algorithms based on static geometric topology to frequently fail.

[0005] In summary, existing navigation technologies lack a dynamic cooperative positioning architecture that can effectively cope with environments without fixed base stations. They also fail to solve the problems of adaptive feature extraction and multi-source kinematic cross-validation of traditional stellar image processing algorithms under conditions of strong optical noise, strong distortion, and limited computing power in polar environments. As a result, they are unable to meet the high reliability navigation requirements of polar UAVs in extreme environments. Summary of the Invention

[0006] This invention provides a polar unmanned aerial vehicle (UAV) navigation method based on image-assisted localization: Obtain the approximate 3D position and approximate attitude quaternions of the polar UAV; Real polar starry sky images are collected and matched with theoretical virtual starry sky images to construct a continuous optical residual strain field. Calculate the divergence and curl of the optical residual strain field, identify and propose the refractive distortion region and the aurora interference region, and obtain the initial starry sky region; Extract multidimensional features from the initial star region, calculate Mahalanobis distance, iteratively correct the effective star region boundary based on Mahalanobis distance, and obtain candidate sky regions; The actual and theoretical optical flow fields within the candidate sky regions are extracted and verified by combining the three-dimensional translational velocity vectors calculated by the mobile local area network. False star regions are eliminated to obtain the optimal star region. The final absolute attitude is obtained by performing star map calculation within the optimal starry sky region, and then fed back to the error state extended Kalman filter to complete the final accurate positioning.

[0007] Constructing a continuous optical residual strain field specifically includes: projecting visible stars onto the image pixel plane based on coarse 3D position and coarse attitude quaternions to obtain theoretical virtual pixel coordinates. , generate containing A theoretical virtual star map of 100 stars; local highlight extrema extraction is performed on real polar star images to obtain the actual centroid pixel coordinates of stars. And within the tolerance radius, the theoretical virtual pixel coordinates Perform matching and calculate its discrete pixel displacement residual vector. The discrete pixel displacement residual vector is diffused into a continuous two-dimensional optical residual strain field. Its coordinates at any pixel in the image The strain vector at that point is defined as: in, The number of successfully matched star pairs. The smoothing variance parameter of the kernel function. and Representing the optical strain field at shaft and Continuous offset components in the axial direction.

[0008] The divergence and curl of the optical residual strain field are calculated to identify and propose refractive distortion regions and auroral interference regions, thereby obtaining the initial starry sky region. Specifically, this includes calculating the divergence function of the optical residual strain field at each pixel. and curl function Define the zenith projection unit direction field. Calculate the characteristic quantity of refraction distortion : Calculate aurora characteristic quantities : in, The divergence coefficient; when Greater than the set refractive distortion threshold When, it is determined to be a region of refractive distortion; when Greater than the set chaos threshold At that time, it was determined to be an aurora interference zone; a topology partitioning mask matrix was generated based on the decision result. The areas that were not removed were taken as the initial starry sky area.

[0009] Multidimensional features include at least the proportion of low-frequency diffuse luminescence energy. Distance matching the initial topology The extraction process includes: dividing the initial starry sky region into... Local dynamic candidate subgrids For subgrid images Perform a two-dimensional discrete Fourier transform to define the proportion of low-frequency diffuse luminescence energy. for: in, Frequency domain coordinates; Extract the brightest Each local extremum pixel is converted into a unit observation vector. Extract its star-angle distance topological vector : Topological vector set based on theoretical candidate star library Calculate the initial topology matching distance .

[0010] Multidimensional features also include stellar sharpness feature parameters. With the final topological distance The extraction process includes: When detected and At that time, a two-dimensional Butterworth spatial high-pass filter is introduced for frequency domain reconstruction to obtain a high-frequency enhanced image. : in, and These represent the two-dimensional Fast Fourier Transform and its inverse transform, respectively. For the set cutoff frequency, The filter order; Calculate the zero-order moment of the target patch region in the high-frequency enhanced image. First-order moment and second-order central moment and The equivalent two-dimensional Gaussian variance is defined as... Calculate the stellar sharpness characteristic parameters Sub-pixel-level centroid coordinates are obtained based on the first and zeroth moments, and the fine-grained feature star-angle distance topological vector is recalculated. , and the target star catalog vector The difference is calculated to output the final high-precision topology distance. .

[0011] Calculating the Mahalanobis distance and iteratively correcting the effective star field boundary based on the Mahalanobis distance specifically includes: calculating the average observation elevation angle of the subgrid in inertial space. Construct the theoretical baseline curve for elevation angle-PSF ambiguity: Constructing joint observation vectors and the prior mean vector Based on polar physical constraints and coupling covariance matrix Calculate Mahalanobis distance : Iterative gradient correction of the star region boundary is performed using Mahalanobis distance: in, For ideal maximum sharpness, The polar atmospheric thickness scale angle. For the boundary learning rate, This is the Mahalanobis distance gradient.

[0012] Extracting the theoretical optical flow field within the candidate sky region specifically includes: obtaining the inertial angular velocity output by the airborne high-frequency gyroscope, and calculating the true rotation angular velocity of the camera coordinate system in inertial space by combining the camera's extrinsic parameters. For any image pixel in the candidate sky region Predict its theoretical optical flow vector : in, The calibrated equivalent focal length for the camera.

[0013] The verification was performed using the three-dimensional translational velocity vector calculated by the mobile local area network, and pseudo-star regions were eliminated. Specifically, this included: tracking and extracting the actual motion optical flow field within the candidate sky region. The residual strain field was calculated. ; Obtain the camera translation velocity vector mapped from the local area network Combining the theoretical spatial distribution Jacobian matrix of translation parallax Calculate the inner product correlation coefficient : like greater than the set motion threshold If the motion trajectory does not conform to the assumption of infinity, then the pseudo-starry sky region is determined and eliminated, and the optimal starry sky region is obtained.

[0014] Obtaining the coarse 3D position and coarse attitude quaternions of the polar UAV specifically includes: combining the 3D position, 3D velocity, and 3D acceleration of the mobile base station with the state of the UAV into a joint state vector; constructing an error state extended Kalman filter framework using the polar low-adhesion dynamics model and the cooperative ranging compensation model, and outputting the coarse 3D position and coarse attitude quaternions of the polar UAV as the prior state for joint estimation.

[0015] The present invention also provides a polar unmanned aerial vehicle (UAV) navigation system based on image-assisted positioning, the system comprising: Acquisition module: Acquires the approximate 3D position and approximate attitude quaternions of the polar UAV; Initial starry sky acquisition module: Acquire real polar starry sky images, combine them with theoretical virtual starry sky images for matching and subtraction, and construct a continuous optical residual strain field; Calculate the divergence and curl of the optical residual strain field, identify and propose the refractive distortion region and the aurora interference region, and obtain the initial starry sky region; Optimal star sky acquisition module: Extracts multi-dimensional features within the initial star sky region, calculates Mahalanobis distance, iteratively corrects the effective star sky region boundary based on Mahalanobis distance, and obtains candidate sky regions; The actual and theoretical optical flow fields within the candidate sky regions are extracted and verified by combining the three-dimensional translational velocity vectors calculated by the mobile local area network. False star regions are eliminated to obtain the optimal star region. The positioning module performs star map calculations within the optimal starry sky region to obtain the final absolute attitude, and feeds the final absolute attitude back into the error state extended Kalman filter to complete the final accurate positioning.

[0016] This application provides a polar UAV navigation method and system based on image-assisted positioning. By constructing an overall framework of ground-based mobile local area network collaborative positioning and star image-assisted positioning, the navigation reliability of UAVs under polar high-latitude strong magnetic fields and extreme weather conditions is improved. At the star image-assisted positioning level, to address the atmospheric refraction distortion caused by the intense aurora and strong near-surface temperature inversion layer in the polar regions, this application proposes an image feature extraction and topological region division mechanism based on optical residual strain field. Relying on the reference of coarse pose generation, the system uses computer vision technology to precisely compare real-world star images, constructing a continuous vector field reflecting the actual physical offset. Through spatial topological differential analysis of this vector field, this scheme can accurately separate the aurora chaotic interference region and the low-altitude refraction distortion region. This overcomes the technical defects of traditional static panoramic image matching algorithms, which are prone to image mismatch and attitude jumps under complex polar optical noise, and greatly improves the image matching and pattern recognition accuracy under extreme light pollution conditions on the server side.

[0017] To further filter out artifacts from low-altitude ice crystal reflections and the hidden interference from faint auroral coverage in the polar regions, this application constructs a multi-dimensional nonlinear mutual constraint evaluation mechanism based on the physical spatial characteristics of the polar regions. The system does not rely on isolated image brightness thresholds, but instead fuses features such as the proportion of low-frequency diffuse energy, point spread function ambiguity, observation elevation angle, and geometric topological matching distance. Combining the optical principle of low-elevation-angle blurring caused by polar inversion layers, image enhancement processing is performed through frequency domain reconstruction. The system calculates the Mahalanobis distance between the multi-dimensional feature distributions to cross-falsify the features of false stars with contradictory physical properties. Utilizing this distance for dynamic iteration and continuous inward squeezing of the usable boundary, the system achieves efficient screening of high-quality stellar image targets, ensuring accurate image classification and core feature extraction even in harsh environments.

[0018] Furthermore, this application introduces a kinematic closed-loop mutual verification mechanism, which rigorously maps the actual 3D translational velocity of the UAV decoded and output by the underlying sensor network to the image optical flow prediction model of the airborne camera, and cross-compares it with the continuous motion optical flow field of the actual tracked star image. Utilizing the kinematic residual strain between the actual ground physical velocity measurement and the parallax of the aerial image translation, pseudo-star clouds that do not conform to the laws of inertial space can be identified and eliminated. This yields optimal fine-grained sky region division, thereby outputting a highly reliable absolute attitude anchoring benchmark, fundamentally ensuring the absolute safety of long-endurance complex missions of polar UAVs. Attached Figure Description

[0019] Figure 1 A flowchart of image-assisted localization-based navigation for polar unmanned aerial vehicles; Figure 2This is a graph showing the relationship between the elevation angle and sharpness of polar stars. The black line represents the theoretical baseline, the blue dots represent actual observational data that conforms to the pattern, and the red cross in the upper left corner represents abnormal interference points. Detailed Implementation

[0020] This embodiment provides a polar UAV navigation method based on image-assisted positioning. In polar scientific research environments, mobile local area network base stations are typically mounted on actively moving vehicles such as snowmobiles or tracked vehicles. Due to the extremely low coefficient of adhesion on polar ice and snow surfaces, mobile base stations inevitably experience high-frequency track / tire slippage and nonlinear vibrations during active movement. Using conventional static node models or simple uniform velocity models will cause the filter to diverge rapidly. Furthermore, the extremely cold and dry atmospheric conditions of the polar regions alter the tropospheric propagation delay of radio frequency signals, and the vast, flat, unobstructed ice surface easily generates strong low-altitude grazing multipath effects on ranging radio waves.

[0021] To address this, this embodiment constructs an error-state extended Kalman filter architecture that considers active maneuvering under low polar adhesion and adaptive compensation for ice surface multipath. The specific implementation steps are as follows: To accurately track the dynamic coordinates of the active mobile base station, this method treats the three-dimensional position, three-dimensional velocity, and three-dimensional maneuvering acceleration of the mobile base station as estimable state variables, and combines them with the UAV state to form a high-dimensional joint state vector.

[0022] Define the system in Joint true state vector at time step for: in, The total number of active mobile base stations that constitute a local area network, ; For base station indexing, ; Let the true state vector of the UAV be defined as follows: ; , , , , These are the quaternions for the UAV's three-dimensional position vector, three-dimensional velocity vector, attitude, three-dimensional zero-bias vector of the onboard accelerometer, and three-dimensional zero-bias vector of the onboard gyroscope. For the first The real state vector of a mobile base station is defined as follows: ; , These are the three-dimensional position vector and the three-dimensional velocity vector of the mobile base station, respectively. It is a three-dimensional acceleration vector that couples the active driving force of the mobile base station with the adhesion force of the ice surface.

[0023] Synchronous definition of joint error state vector for: in, Let be the three-dimensional attitude error vector of the UAV. This is the error state vector corresponding to the physical quantity.

[0024] During intervals when no external collaborative observation data is available, the system performs high-frequency time updates to the joint state. (This is for the UAV error state.) The following is a recursive derivation based on the continuous-time error propagation equation for inertial navigation: in, This is the error state transition matrix for the UAV. The system noise driving matrix, White noise for airborne inertial sensors.

[0025] For mobile base stations, due to the frequent slippage and irregular acceleration of snowmobiles when actively driving on polar icy and snowy surfaces, this invention adds the sudden changes in slippage speed caused by the low adhesion of the polar environment to the ice surface noise in the speed integral term. The continuous-time dynamic equations for the mobile base stations are established as follows: The mechanical hysteresis physical delay for executing acceleration commands to polar tracked vehicles. In deep snow or ice environments in the polar regions, tracked vehicles have a much lower maneuver frequency than conventional wheeled vehicles due to their high mechanical damping. Let the maneuver frequency of the polar vehicle be... ,but In polar environments Values .

[0026] Using the state transition matrix For the joint error state covariance matrix Perform time step arrive Update: in, Let be the error state covariance matrix. For inclusion , and The system process noise covariance matrix.

[0027] When drones are Received at time When coordinating radio ranging signals from multiple mobile base stations, observation equations for the polar environment need to be established. The extreme cold of the polar regions leads to extremely dry and denser lower atmospheres, reducing the propagation speed of electromagnetic waves in such media. This invention introduces a compensation term to define the spatial observations for coordinating ranging. The equation is: in, For Euclidean L2 norm operators; To observe noise; This is the equivalent distance compensation amount caused by the propagation delay at the Earth pole.

[0028] in, This refers to the real-time absolute ambient temperature of the polar regions. For air pressure, The refractive index constant of the polar atmosphere. For drones relative to the first The observation elevation angle of each base station.

[0029] The above observation equations are linearized and expanded to extract the values ​​corresponding to the UAV position. and base station location The partial derivatives are used to construct the measurement Jacobian submatrix. : in, The values ​​for the corresponding velocity, acceleration, and attitude error terms are all assigned values. To include all visible base stations in the network. Stacking to form a global observation matrix .

[0030] Specifically, in order to complete the measurement update within the joint extended Kalman filter framework, the aforementioned partial derivatives need to be accurately mapped to the global high-dimensional state space. The joint error state vector of the system at the current moment is defined. The total dimension is Due to the error status of the drone Includes position, velocity, attitude, and two sets of zero biases, with dimensions of Dimension; The error states (position, velocity, acceleration) of each mobile base station are: Dimension; therefore, total dimension .

[0031] For the Scalar ranging equations for a mobile base station Its corresponding measurement Jacobian submatrix It is a size of The row vectors. Specifically, Only at the location corresponding to the drone (index (column) and corresponding to the first Location of each base station (index The block positions of the column are filled with non-zero partial derivative values, and the remaining block matrix parts corresponding to the UAV's velocity, attitude, and the velocities and accelerations of all base stations are assigned zero vectors. .

[0032] Therefore, a single observation matrix for: When polar drones are at the same time Receive all visible signals on the local area network. When multiple mobile base stations send concurrent cooperative ranging signals, all row vectors are... Rows are stacked according to the base station index order to form a structure of size. Global joint observation matrix : Considering the strong multipath effect of radio signals on the vast, flat ice surfaces of the polar regions, when drones fly at low altitudes or are far from base stations, the superposition of reflected signals from the ice surface and direct signals can lead to severe ranging jumps.

[0033] This embodiment constructs an adaptive dynamic measurement noise covariance matrix. . It is a diagonal matrix, and the first element on its diagonal is... variance of each element Defined as relating to the observation elevation angle Functions: in, The baseline distance measurement variance, This is the multipath reflection coefficient constant. The additional penalty term introduced in the formula... The purpose is to automatically increase the measurement noise variance and reduce the observation weight during unreliable periods when the system estimates that the mobile base station is undergoing drastic acceleration, deceleration, or slippage (i.e., the acceleration L2 norm increases and the base station antenna will vibrate violently due to vehicle pitch). This is the adjustment coefficient.

[0034] Calculate the Kalman gain matrix : Calculate the joint error state vector : in, This represents the actual ranging vector of the collaborative network. The predicted ranging vector is calculated by substituting the prior state of the system into the compensation model above.

[0035] Update covariance matrix : The calculated Each corresponds to compensation to the joint true state vector After the compensation is completed, the joint error state vector is cleared to zero for physical quantities such as position, velocity, acceleration and attitude.

[0036] Next, this embodiment provides a polar UAV navigation method based on image-assisted localization. It should be noted that this embodiment is combined with the previous embodiment, and the posterior state output in the previous embodiment is used as the prior coarse pose in this embodiment. When performing missions at night in polar regions, traditional star map matching algorithms typically search for bright features across the entire sky. However, intense polar auroral activity generates large areas of low-frequency bright regions, which are easily mistaken by algorithms for dense star regions. Simultaneously, the frigid air currents on the polar ice caps create strong temperature inversions, causing severe atmospheric refraction of starlight near the horizon and resulting in regular shifts in star positions. Directly performing star map matching on the entire image is not only computationally expensive but also highly susceptible to attitude calculation failures due to the introduction of auroral noise or refraction distortion points.

[0037] Therefore, this embodiment is based on the sky topology region division mechanism of optical residual strain field, and the specific implementation steps are as follows: First, obtain the image acquisition time of the drone at the current image acquisition moment. Prior rough three-dimensional position With prior rough pose quaternions , And satisfies the unitization constraint The coarse pose can be output by the local area network cooperative positioning in the previous embodiment, or provided by the UAV inertial navigation module.

[0038] Using a full-sky star catalog database pre-stored in the UAV's onboard computer, information on stars currently within the UAV's upper hemisphere field of view is extracted. Let the... The unit direction vector of a visible star in the Earth's inertial coordinate system (ECI) is: ,in For star index labels, , This represents the total number of visible stars.

[0039] Based on the prior coarse pose, the stellar direction vector is transformed into the camera coordinate system of the airborne astrophotography camera: in, For three-dimensional position With time The coordinate transformation matrix determined from the inertial frame to the navigation frame. For attitude quaternions The rotation matrix determined from the navigation system to the fuselage side. The intrinsic mounting extrinsic parameter matrix from the aircraft to the star camera.

[0040] Specifically, according to the rules of quaternion algebra, the rotation matrix for the transformation from the navigation system to the machine system. for: Furthermore, in high-latitude polar regions, geographic longitude lines exhibit strong geometric convergence with increasing latitude. In this case, even small perturbations to the position vector can cause a drastic rotation of the azimuth angle of the local reference coordinate system. This invention employs North-East-Earth as the local navigation system, transforming the prior coarse three-dimensional position... Convert to polar geographic coordinates: latitude With longitude .

[0041] Simultaneously, extract the image acquisition time. Calculate the Greenwich Mean Time (GMT) corresponding to that moment, denoted as . Therefore, the absolute inertial right ascension angle of the polar meridian where the UAV is currently located can be calculated. : Based on polar geographic coordinates With inertial right ascension angle Coordinate transformation matrix from Earth's inertial frame to the northeastern navigation frame for: In the matrix above, when the drone is operating in the polar region (latitude) )hour, and This means Matrix of longitude The sensitivity is amplified nonlinearly. Traditional star map projection algorithms often cause pixel-level global shifts in star projection positions in polar regions due to neglecting this nonlinear coupling. This embodiment, through rigorous analysis of the aforementioned analytical matrix, ensures that the theoretical virtual pixel coordinates can accurately absorb the rotation caused by the convergence of polar meridians, thereby guaranteeing the subsequent optical residual strain field. The study only isolates the physical refraction of the atmosphere and the interference of aurora diffusion, without incorporating the spurious mathematical residuals caused by the polar coordinate system.

[0042] Combined with the intrinsic parameter matrix of the starry sky camera ,Will Projected onto the image pixel plane, we obtain the first pixel. Theoretical virtual pixel coordinates of a star : This generates a set of coordinates. A theoretical virtual star map.

[0043] Control the airborne star camera to capture real polar starry sky images that include auroras and refraction interference. Local highlight extrema extraction was performed on the image to obtain the set of actual star centroid pixel coordinates. ,in This represents the total number of bright connected regions actually extracted, including real stars, aurora spots, and ice crystal reflection artifacts.

[0044] Based on the Euclidean distance in coordinate space, within the tolerance radius Inside, to and Perform nearest neighbor matching. For successfully matched point pairs, calculate their discrete pixel displacement residual vectors. : Subsequently, this invention introduces radial basis function interpolation to diffuse the discrete displacement residuals into a continuous two-dimensional optical residual strain field covering the entire image size. Define arbitrary pixel coordinates of the image. The strain vector at that point is: in, The number of successfully matched star pairs. represents the smoothing variance parameter of the kernel function. and These represent the optical strain field at the image pixel. shaft and Continuous offset components in the axial direction.

[0045] After constructing a continuous optical residual strain field, the divergence and curl of the vector field at each pixel are calculated.

[0046] divergence function Used to measure the degree of divergence or convergence of pixel displacement in this region: curl function Used to measure the degree of local rotation or shearing of pixel displacement in this region: based on and Generate region mask matrix .

[0047] Due to the influence of the polar ice inversion layer, light is severely bent when passing through the atmosphere. Furthermore, because the refractive index increases sharply as the elevation angle decreases, the actual imaging position of all low-elevation stars on the image will systematically shift and shrink upwards, that is, towards the zenith.

[0048] Therefore, the zenith projection unit direction field is defined. Calculate the characteristic quantity of refraction distortion : The second term represents the projected length of the residual vector towards the zenith. Polar refraction in the second term causes the distribution of star points to become denser, exhibiting a physical convergence phenomenon (i.e., divergence). Therefore, when Greater than the set refractive distortion threshold At that time, the region not only shifted but also compressed towards the zenith, indicating that the region was a polar refraction distortion zone.

[0049] The aurora is not a point light source, but rather a dynamic electromagnetic diffuse plasma in the upper atmosphere. The drift of the aurora causes the high-brightness pseudo-centroids extracted from local images to exhibit disordered, multi-directional random wandering, thereby inducing strong local shearing and chaotic divergence in the optical residual strain field.

[0050] Calculate aurora characteristic quantities : The deformation of aurora artifacts can cause anomalous vortices (high curl) and local bursts (high divergence) in the residual vector field. When Greater than the set chaos threshold At that time, the area was determined to be an aurora interference zone.

[0051] Based on the aforementioned physical mechanism, a topology partitioning mask matrix with the same resolution as the image is constructed. : In the marked Within the starry sky region, its optical residual vector It exhibits extremely low divergence and curl, or contains only uniform translations caused by prior coarse pose errors, and is approximately white noise.

[0052] After deconstructing the physical mechanisms of the panoramic starry sky, this invention obtains and Polluted image pixels and Images of the starry sky.

[0053] Next, this embodiment will specifically describe the reliability assessment and boundary correction method. In the previous embodiment, the aurora interference region and refractive distortion region were initially separated by the optical residual strain field, while the starry sky region, i.e., the mask matrix, was preserved. The area. However, in the real polar scientific research environment, extremely complex optical deceptions may still lurk within the reserve: for example, the strong reflection of drone-borne anti-collision strobe lights on the low-altitude ice surface at close range can form false star points with extremely high gradients; at the same time, some real stars covered by faint auroras, although their backgrounds show low-frequency diffuse reflection characteristics, still have the correct star catalog topological distances.

[0054] If traditional isolated threshold judgments or simple weighted summation are used for elimination, it is easy to mistakenly delete true stars and retain false stars. Therefore, this embodiment is based on the polar star observation credibility map, which uses the physical properties of stars to mutually support and falsify each other in space, and uses Mahalanobis distance to iteratively correct the boundary information.

[0055] Output starry sky area Further divided into Local dynamic candidate subgrids ( ).

[0056] For each subgrid Based on the prior coarse pose of the UAV, the average observation elevation angle of the subgrid in inertial space is calculated. At the same time, for Initial feature extraction is performed on the pixel image within the image, and the proportion of low-frequency diffuse reflection energy of the background intensity is calculated. This is used to characterize the basic noise level of suspected auroras, and the initial topological matching distance within the sub-grid is extracted through rapid pre-matching with the airborne local star catalog. , The smaller the value, the more the geometric topology of the representation matches the real starry sky.

[0057] At the same time, Pixel images within Initial feature extraction is performed, specifically calculating the proportion of low-frequency diffuse luminescence energy in auroras. and initial topology matching distance .

[0058] Polar auroras are essentially ionized light emitted from the upper atmosphere, exhibiting a diffuse characteristic of being concentrated in the low-frequency band in image spatial frequency representation. To measure this baseline noise level, subgrid images were analyzed. Perform a two-dimensional discrete Fourier transform to obtain its power spectral density matrix in the spatial frequency domain. : in, These are the spatial pixel coordinates within the subgrid. For frequency domain coordinates, and These represent the width and height of the grid pixels, respectively. The cutoff spatial frequency limit for diffuse aurora is... Define the proportion of low-frequency diffuse luminescence energy. This is the ratio of the integrated energy within the low-frequency cutoff circle to the total global energy. The ratio is between Between these values, the larger the value, the more deeply the background of the region is masked by low-frequency ionization emission of the aurora.

[0059] Calculate the initial topology matching distance : exist Extract the brightest Local extreme value pixels, preferred The star map triangle features are constructed, and their pixel coordinates are converted into unit observation vectors in the camera coordinate system. , Extracting the star-angle distance topological vector. Its elements are the spatial cosine angular distances between each vertex: Based on the current coarse pose retrieval of the airborne star catalog, a set of topological vectors for the theoretical candidate star library is extracted. , This represents the possible combined indexes in the star list. Initial topological matching distance. The minimum Euclidean distance between the two is: Let the star catalog combination corresponding to this minimum distance be the target star group. In polar optical noise, The smaller the value, the more accurately it indicates that even with aurora background noise, the geometric configuration of the bright spots still strictly conforms to the absolute relative positions of the stars in the real night sky. Frequency-domain adaptive reconstruction of aurora diffuse luminescence and topological constraints: This step executes the first layer of the constraint mechanism: the mandatory constraint on the diffuse properties and topological structure of the aurora. The aurora region exhibits low-frequency diffuse characteristics optically, but the high-frequency point light source information of real stars may be obscured within it.

[0060] Pair of subgrids The following logical judgment is made: The system detects that the background sharpness in the region is extremely low, i.e., the proportion of low-frequency diffuse energy meets the following condition. The area appears to be covered by a faint aurora, and the bright spot combinations extracted from it conform to the star catalog topology, meaning the initial topological matching distance satisfies... If the subgrid is not discarded, it is instead identified as a high-value aurora diffuse region containing real stars, and a two-dimensional Butterworth space high-pass filter is introduced into the subgrid image. Frequency domain reconstruction was performed to filter out low-frequency auroras and restore high-frequency stellar features. The reconstructed high-frequency enhanced image. for: in, and These represent the two-dimensional Fast Fourier Transform and its inverse transform, respectively. For frequency domain coordinates, The cutoff frequency is dynamically set based on the frequency of aurora fluctuations. This represents the filter order.

[0061] After high-frequency enhancement or if the condition is not triggered and the original image is preserved, The final sub-pixel centroid extraction is performed within the region, and the reciprocal of the Gaussian variance of the point spread function of the bright spot in that region is fitted as the final stellar sharpness feature parameter. And output the final topological distance. .

[0062] Specifically, stellar sharpness feature parameters are calculated based on the image center moment. : For high-frequency enhanced images The extracted first Target patch area Calculate its zeroth moment and first moment , : This yields sub-pixel-level centroid coordinates unaffected by low-frequency aurora gradient stretching. Further calculations were performed to determine the second-order central moment, which reflects the energy diffusion range of the light spot. and : The PSF of the light spot is defined as the equivalent two-dimensional Gaussian variance. If the patch is a near-range reflection artifact, its light spot will exhibit a large, discrete distribution. It will rise sharply. Define the final stellar sharpness feature parameter for this subgrid. for The reciprocal of the mean of the Gaussian variances of the PSF of each star point: Output the final high-precision topology distance : After filtering out aurora background gradient interference and calculating sub-pixel centroids with extremely high precision Then, it is converted back into a unit observation vector. Recalculate the topological vector of the characteristic star angular distance. and the locked target star table vector Calculate the difference and output the final topology matching distance. : Real Star Group The convergence will be reduced to sub-arcsecond precision; while polar artifact noise that coincidentally forms a similar initial topology cannot pass this sub-pixel-level rigorous comparison and will be eliminated in the subsequent covariance evaluation.

[0063] Construction of the covariance matrix coupled with the physical laws of the near-surface polar inversion layer: This step implements the second layer of the constraint mechanism: the physical constraint on the observation elevation angle (B) and stellar sharpness (A). In polar environments, the extreme cold at the surface leads to a strong near-surface temperature inversion layer. The closer to the horizon (low elevation angle region), the greater the atmospheric refractive index gradient, and the more severely the star's point spread function (PSF) will be stretched and blurred.

[0064] Based on this natural physical law, a theoretical benchmark curve for elevation angle-PSF ambiguity is constructed. : in, This represents the ideal maximum sharpness when there is no distortion in the zenith direction. This is the polar atmospheric thickness scaling angle. In the theoretical baseline curve of elevation angle-PSF ambiguity, low elevation angles correspond to low sharpness.

[0065] This invention constructs a joint observation vector from the extracted multidimensional features. At the same time, define the prior mean vector. .

[0066] In the process of mutually falsifying the observed elevation angle (B) and stellar sharpness, a polar physical constraint coupling covariance matrix is ​​constructed. : In this covariance matrix, the core lies in the cross-correlation coefficient. If in a low elevation angle area ( ), and extracted an exceptionally sharp bright spot with a very high gradient (i.e. This severely contradicts the natural phenomenon of polar inversions causing obscuration. (Through off-diagonal elements) The constraint mapping will be falsified by the system at the underlying mathematical logic level, and it will be determined to be a reflection artifact of the drone anti-collision strobe light.

[0067] Based on the joint observation vector and physical coupling covariance matrix constructed above, for each subgrid... Calculate the Mahalanobis distance of its feature point distribution. : The smaller the Mahalanobis distance, the more the optical characteristics of the region conform to the objective laws verified by the current UAV attitude and the multiple physical environments of the polar region. Regarding the aforementioned low-elevation and exceptionally sharp reflection artifacts, the calculated... It will tend to infinity, and thus will be automatically assigned an extremely low credibility weight.

[0068] Finally, the Mahalanobis distance is used to define the effective star region boundary of the output. Iterative correction based on gradient: in, For the boundary learning rate, This represents the Mahalanobis distance gradient of the local neighborhood of a pixel.

[0069] Based on the constraints of polar physics, the Mahalanobis distance continuously squeezes and corrects the boundary inward, peeling away those contradictory areas of the starry sky, such as topologically correct but with incorrect sharpness, or low elevation angle but too clear, layer by layer, until it converges to obtain the final starry sky solution region that conforms to the natural laws of the polar regions.

[0070] In the harsh polar environment, some high-altitude drifting clouds can form a localized bright topological structure that is very similar to that of stars under certain lighting conditions. However, real stars are at infinity, and their image motion (optical flow) on the camera's focal plane is only affected by the drone's attitude rotation, not by the drone's translational motion. Meanwhile, the polar ice crystal clouds disguised as stars are at a finite altitude (with a finite depth), and their image motion includes parallax optical flow caused by the drone's translational speed.

[0071] Candidate sky region sequences and coarse attitude calculations were extracted. After initial segmentation, three high-confidence candidate regions were selected in the sky, denoted as Region 1 (…). ), Zone 2 ( ) and Zone 3 ( ).

[0072] Extracting continuous time periods A sequence of starry sky images within the area. Selected area 1 ( Using this as an initial reference benchmark, preliminary Wahba attitude calculations are performed based on the coordinates of the star points within it to obtain... Coarse attitude rotation matrix of the UAV from navigation system to machine system .

[0073] Within the synchronized time window of image acquisition, a mobile local area network base station deployed on the polar ice surface is used to perform joint Doppler frequency shift and carrier phase velocity measurements on the UAV. The three-dimensional translational velocity vector of the UAV relative to the polar ice surface navigation system is then calculated. : in, These represent the actual ground-to-ground flight speeds of the drone in the north, east, and ground directions, respectively.

[0074] Combined with the coarse attitude calculated from region 1 and the fixed body-to-camera extrinsic parameter matrix The ground velocity vector is mapped onto the 3D camera coordinate system of the star camera to obtain the camera translation velocity vector. : Obtain the inertial angular velocity output by the airborne high-frequency gyroscope By combining the camera's extrinsic parameters, the true rotational angular velocity of the camera coordinate system in inertial space is calculated. .

[0075] Since the real star is at infinity, the drone's translational speed It does not produce any parallax. Therefore, for zone 2 ( ) and Zone 3 ( Any image pixel in ) Predict its theoretical optical flow vector : in, The calibrated equivalent focal length for the camera.

[0076] Using the Lucas-Kanade (LK) continuous pyramid optical flow algorithm, the actual motion trajectories of bright spots in regions 2 and 3 of the image sequence were tracked, and their actual motion optical flow fields were extracted. By subtracting the actual optical flow from the theoretical optical flow, the residual strain field of the local area network velocity-starlight optical flow is calculated. : If the lurking creature in Zone 2 is at a limited altitude The images of polar high-altitude ice crystal clouds are affected by the translational motion of the drone, resulting in a translational parallax optical flow field. According to projection geometry, the theoretical spatial distribution of translational parallax is represented by the Jacobian matrix. Speed ​​of camera coordinates and local area network mapping Decide: Calculate the residual strain field Parallax structure with local area network The inner product correlation coefficient between them : Drone speed calculated by local area network This is the true value given by the ground. If region 2 represents the real starry sky, its residual... Caused solely by camera white noise, and having no directional correlation with the drone's translational speed, at this time... If region 2 is a star disguised as a polar high-altitude ice crystal cloud, its residual... This must be caused by the parallax created when the drone moves at a horizontal speed across the clouds, resulting in its difference from... They are highly linearly aligned in the mathematical direction, at this point .

[0077] Set motion threshold .when At that time, pseudo-star regions whose motion trajectory does not conform to the assumption of infinity are eliminated.

[0078] Through the aforementioned cross-validation, all pseudo-star regions affected by parallax interference from ice crystal clouds are eliminated, yielding the optimal star region. Within this confirmed optimal star region, the highest precision sub-pixel star map calculation is performed to obtain the 3D attitude. This final absolute attitude is then used as the observation input and fed back to the joint error state extended Kalman filter to achieve global convergence and final precise positioning of the inertial navigation error and local area network coordinate system drift.

[0079] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), Rambus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0080] In this specification, the same or similar parts between the various embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the descriptions of the embodiments described later are relatively simple, and relevant parts can be referred to the descriptions of the foregoing embodiments.

[0081] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A polar unmanned aerial vehicle (UAV) navigation method based on image-assisted localization, characterized in that, include: Obtain the approximate 3D position and approximate attitude quaternions of the polar UAV; Real polar starry sky images are collected and matched with theoretical virtual starry sky images to construct a continuous optical residual strain field. Calculate the divergence and curl of the optical residual strain field, identify and propose the refractive distortion region and the aurora interference region, and obtain the initial starry sky region; Extract multidimensional features from the initial star region, calculate Mahalanobis distance, iteratively correct the effective star region boundary based on Mahalanobis distance, and obtain candidate sky regions; The actual and theoretical optical flow fields within the candidate sky regions are extracted and verified by combining the three-dimensional translational velocity vectors calculated by the mobile local area network. False star regions are eliminated to obtain the optimal star region. The final absolute attitude is obtained by performing star map calculation within the optimal starry sky region, and then fed back to the error state extended Kalman filter to complete the final accurate positioning.

2. The polar UAV navigation method based on image-assisted positioning according to claim 1, characterized in that, Constructing a continuous optical residual strain field specifically includes: projecting visible stars onto the image pixel plane based on coarse 3D position and coarse attitude quaternions to obtain theoretical virtual pixel coordinates. , generate containing A theoretical virtual star map of 100 stars; local highlight extrema extraction is performed on real polar star images to obtain the actual centroid pixel coordinates of stars. And within the tolerance radius, the theoretical virtual pixel coordinates Perform matching and calculate its discrete pixel displacement residual vector. The discrete pixel displacement residual vector is diffused into a continuous two-dimensional optical residual strain field. Its coordinates at any pixel in the image The strain vector at that point is defined as: ； in, The number of successfully matched star pairs. The smoothing variance parameter of the kernel function. and Representing the optical strain field at shaft and Continuous offset components in the axial direction.

3. The polar UAV navigation method based on image-assisted positioning according to claim 2, characterized in that, The divergence and curl of the optical residual strain field are calculated to identify and propose refractive distortion regions and auroral interference regions, thereby obtaining the initial starry sky region. Specifically, this includes calculating the divergence function of the optical residual strain field at each pixel. and curl function Define the zenith projection unit direction field. Calculate the characteristic quantity of refractive distortion : ； Calculate aurora characteristic quantities : ； in, The divergence coefficient; when Greater than the set refractive distortion threshold When, it is determined to be a region of refractive distortion; when Greater than the set chaos threshold At that time, it was determined to be an aurora interference zone; a topology partitioning mask matrix was generated based on the decision result. The areas that were not removed were taken as the initial starry sky area.

4. The polar UAV navigation method based on image-assisted positioning according to claim 1, characterized in that, Multidimensional features include at least the proportion of low-frequency diffuse luminescence energy. Distance matching the initial topology ; The extraction process includes: dividing the initial starry sky region into... Local dynamic candidate subgrids For subgrid images Perform a two-dimensional discrete Fourier transform to define the proportion of low-frequency diffuse luminescence energy. for: ； in, Frequency domain coordinates; Extract the brightest Each local extremum pixel is converted into a unit observation vector. Extract its star-angle distance topological vector : ； ； Topological vector set based on theoretical candidate star library Calculate the initial topology matching distance .

5. A polar UAV navigation method based on image-assisted positioning according to claim 4, characterized in that, Multidimensional features also include stellar sharpness feature parameters. With the final topological distance ; Its extraction process includes: When detected and At that time, a two-dimensional Butterworth spatial high-pass filter is introduced for frequency domain reconstruction to obtain a high-frequency enhanced image. : ； in, and These represent the two-dimensional Fast Fourier Transform and its inverse transform, respectively. For the set cutoff frequency, The filter order; Calculate the zero-order moment of the target patch region in the high-frequency enhanced image. First-order moment and second-order central moment and The equivalent two-dimensional Gaussian variance is defined as... Calculate the stellar sharpness characteristic parameters Sub-pixel-level centroid coordinates are obtained based on the first and zeroth moments, and the fine-grained feature star-angle distance topological vector is recalculated. , and the target star catalog vector The difference is calculated to output the final high-precision topology distance. .

6. A polar UAV navigation method based on image-assisted positioning according to claim 5, characterized in that, Calculating the Mahalanobis distance and iteratively correcting the effective star field boundary based on the Mahalanobis distance specifically includes: calculating the average observation elevation angle of the subgrid in inertial space. Construct the theoretical baseline curve for elevation angle-PSF ambiguity: ; Constructing joint observation vectors and the prior mean vector Based on polar physics-constrained coupling covariance matrix Calculate Mahalanobis distance : ； Iterative gradient correction of the star region boundary is performed using Mahalanobis distance: ； in, For ideal maximum sharpness, The polar atmospheric thickness scale angle. For the boundary learning rate, This is the Mahalanobis distance gradient.

7. A polar UAV navigation method based on image-assisted positioning according to claim 1, characterized in that, Extracting the theoretical optical flow field within the candidate sky region specifically includes: obtaining the inertial angular velocity output by the airborne high-frequency gyroscope, and calculating the true rotation angular velocity of the camera coordinate system in inertial space by combining the camera's extrinsic parameters. For any image pixel in the candidate sky region Predict its theoretical optical flow vector : ； in, The calibrated equivalent focal length for the camera.

8. A polar UAV navigation method based on image-assisted positioning according to claim 7, characterized in that, The verification was performed using the three-dimensional translational velocity vector calculated by the mobile local area network, and pseudo-star regions were eliminated. Specifically, this included: tracking and extracting the actual motion optical flow field within the candidate sky region. The residual strain field was calculated. ; Obtain the camera translation velocity vector mapped from the local area network Combining the theoretical spatial distribution Jacobian matrix of translation parallax Calculate the inner product correlation coefficient : ； like greater than the set motion threshold If the motion trajectory does not conform to the assumption of infinity, then the pseudo-starry sky region is determined and eliminated, and the optimal starry sky region is obtained.

9. A polar UAV navigation method based on image-assisted positioning according to claim 1, characterized in that, Obtaining the coarse 3D position and coarse attitude quaternions of the polar UAV specifically includes: combining the 3D position, 3D velocity, and 3D acceleration of the mobile base station with the state of the UAV into a joint state vector; constructing an error state extended Kalman filter framework using the polar low-adhesion dynamics model and the cooperative ranging compensation model, and outputting the coarse 3D position and coarse attitude quaternions of the polar UAV as the prior state for joint estimation.

10. A polar unmanned aerial vehicle (UAV) navigation system based on image-assisted positioning, characterized in that, The system includes: Acquisition module: Acquires the approximate 3D position and approximate attitude quaternions of the polar UAV; Initial starry sky acquisition module: Acquire real polar starry sky images, combine them with theoretical virtual starry sky images for matching and subtraction, and construct a continuous optical residual strain field; Calculate the divergence and curl of the optical residual strain field, identify and propose the refractive distortion region and the aurora interference region, and obtain the initial starry sky region; Optimal star sky acquisition module: Extracts multi-dimensional features within the initial star sky region, calculates Mahalanobis distance, iteratively corrects the effective star sky region boundary based on Mahalanobis distance, and obtains candidate sky regions; The actual and theoretical optical flow fields within the candidate sky regions are extracted and verified by combining the three-dimensional translational velocity vectors calculated by the mobile local area network. False star regions are eliminated to obtain the optimal star region. The positioning module performs star map calculations within the optimal starry sky region to obtain the final absolute attitude, and feeds the final absolute attitude back into the error state extended Kalman filter to complete the final accurate positioning.

Images