A heading angle disambiguation method based on gradient change of polarization angle

Abstract

The application discloses a heading angle disambiguation method based on polarization angle gradient change, extracts the asymmetric information of polarization angle change by analyzing the gradient characteristics of polarization angle space distribution, and realizes the rapid disambiguation of the heading angle by combining the gradient direction consistency test in the dynamic window. The application uses the polarization degree gradient vectors at different positions on the same polarization angle to construct local constraints, directly calculates the unique heading angle through the geometric relationship between the gradient direction and the heading angle, and simultaneously introduces the short-time angular velocity information of the inertial measurement unit to dynamically compensate the gradient change, so that the robustness and real-time performance are significantly improved.

Description

Technical Field

[0001] This invention relates to the field of heading angle measurement technology, and in particular to a heading angle deambiguation method based on polarization angle gradient variation. Background Technology

[0002] Traditional geomagnetic induction-based heading angle measurement methods are susceptible to interference from environmental magnetic fields (such as buildings and metal structures), and their accuracy significantly decreases under complex weather conditions (such as cloudy skies and haze). Therefore, heading angle sensing technology based on sky polarization modes has gradually become an important research direction in the field of autonomous navigation. Polarized navigation, by analyzing the atmospheric polarization angle distribution characteristics, can provide stable heading information without accumulated errors. Polarization sensors have advantages such as small size, low power consumption, and passive sensing, complementing the characteristics of Inertial Navigation Systems (INS): polarization navigation can provide INS with a long-term stable heading reference, suppressing heading angle drift; while INS can provide motion state compensation for polarization navigation in short-term dynamic scenarios, improving anti-interference capabilities.

[0003] Existing methods for calculating polarization heading angles mainly fall into two categories: single-point polarization angle calculation and multi-point polarization mode matching. Single-point polarization angle calculation measures the local polarization angle using a single polarization sensor and then calculates the heading angle using a solar position model or a sky polarization distribution model. Multi-point polarization mode matching, on the other hand, obtains the polarization angle distribution using a multi-directional polarization sensor array, compares it with a theoretical or empirical polarization mode library, and selects the optimal matching result to determine the heading angle. However, these methods have the following limitations in practical applications:

[0004] 1. Ambiguity in heading angle: The calculation of single-point polarization angle depends on the linear mapping relationship between polarization angle and heading angle. However, due to polarization symmetry, the heading angle has an inherent ambiguity of 180 degrees, making it impossible to distinguish between positive and negative directions.

[0005] 2. Poor environmental adaptability: Existing methods rely on high-precision sky polarization models or preset model libraries, but in non-ideal weather (such as local clouds, aerosol scattering) or dynamic occlusion scenarios, polarization patterns are distorted, leading to matching failure or model inaccuracy.

[0006] 3. High computational complexity: Multi-point pattern matching requires traversing a large-scale polarization pattern library, making it difficult to guarantee real-time performance; the heading angle optimization algorithm based on nonlinear models has a slow iterative convergence speed, making it difficult to meet the needs of high-speed motion scenarios.

[0007] 4. Strong external dependence: Existing methods require real-time acquisition of prior information on the sun's position or sky polarization mode, and cannot work independently in environments without satellite signals or in dynamic environments (such as indoors or tunnels). Summary of the Invention

[0008] To address the limitations and shortcomings of existing technologies, this invention provides a heading angle deambiguation method based on polarization angle gradient variation, comprising:

[0009] The polarization camera is calibrated to obtain the focal length and principal point parameter values ​​of the polarization camera;

[0010] The polarization camera is used to acquire sky images, and the degree of polarization and polarization angle information of each point are calculated.

[0011] A sky polarization pattern distribution image is obtained based on the polarization degree information and polarization angle information;

[0012] Based on the parameter values ​​of the polarization camera, obtain the polarization camera coordinate system;

[0013] Calculate the polarization angle in the polarization camera coordinate system according to Stokes' theorem;

[0014] The polarization angle in the polarization camera coordinate system is transformed into the polarization angle in the meridional coordinate system according to the coordinate system redrawing method;

[0015] Obtain the polarization angle image after coordinate system redrawing, and get the polarization angle on the meridian;

[0016] Extract the spatial gradient distribution of the polarization degree information and construct a polarization degree gradient vector field;

[0017] A dynamic window is divided in the meridian coordinate system, and the consistency of the polarization gradient direction within the dynamic window is checked to select a preset gradient region.

[0018] Eliminate the 180° polarization angle ambiguity based on gradient direction asymmetry to generate a deambiguous polarization angle distribution image;

[0019] Short-time angular velocity data from an inertial measurement unit are used to compensate for gradient direction drift within the dynamic window.

[0020] Based on the deambigued polarization angle distribution image, the solar azimuth angle is calculated;

[0021] Calculate navigation baselines based on astronomical yearbooks;

[0022] A unique heading angle is output by the geometric constraint relationship between the gradient vector and the heading angle.

[0023] Optional, also includes:

[0024] The polarization state of light and the intensity of light in different polarization directions are described using four Stokes parameters, expressed in matrix form as follows:

[0025] ;

[0026] Where S0 represents the total incident light intensity, S1 represents the light intensity difference between the x and y components, S2 represents the light intensity difference between the +45° and -45° polarization components, and S3 represents the light intensity difference between the left-handed and right-handed circular polarization components.

[0027] Optionally, under natural light conditions, S3=0, and the direction of sunlight is... At that time, the expression for the light intensity after transmission is as follows:

[0028] ,

[0029] When the polarizer angle is taken as 0°, 45°, 90° and 135° respectively, the expression is as follows:

[0030] ,

[0031] Where S0, S1, and S2 are Stokes parameters.

[0032] Optionally, the expression for the polarization angle in the polarization camera coordinate system is as follows:

[0033] (8),

[0034] Where DoLP is the degree of linear polarization, DoP is the degree of polarization, Aop is the polarization angle, and S1 and S2 represent the light intensities of two sets of linearly polarized light that are perpendicular to each other in the polarization direction. .

[0035] Optionally, the polarization angle is the angle between the direction of electric vector vibration and the reference coordinate system of the polarization camera.

[0036] Optionally, the step of transforming the polarization angle in the polarization camera coordinate system to the polarization angle in the meridional coordinate system according to the coordinate system redrawing method includes:

[0037] Obtain the polarization scattering azimuth angle of observation point M in the polarization camera coordinate system. The polarization scattering azimuth angle of observation point N in the polarization camera coordinate system is obtained. ;

[0038] The polarization angles of observation points M and N in the meridian coordinate system are obtained through geometric transformation, as expressed below:

[0039] (9),

[0040] in, Let M be the polarization angle of the observation point in the meridian coordinate system. Let N be the polarization angle of the observation point in the meridian coordinate system. Let M be the polarization angle of the observation point in the polarization camera coordinate system. Let N be the polarization angle of the observation point in the coordinate system of the polarization camera.

[0041] Optional, also includes:

[0042] The polarization angle image is binarized based on the numerical range of the polarization angle along the meridian, as shown in the following expression:

[0043] (10),

[0044] Where Aop is the polarization angle, and Aop_Gray is the degree of polarization of the polarized grayscale image after binarization, with a value of 0 or 1.

[0045] Optional, also includes:

[0046] The least squares method was used to fit a straight line to the binarized polarization angle image;

[0047] The solar azimuth angle is obtained when the sum of the squares of the residuals at all observation points is minimized.

[0048] Optional, also includes:

[0049] Define a circular area to the left of the meridian. Define a circular area to the right of the meridian. This is used to statistically analyze local polarization characteristics, calculating the sum of polarization degrees within the circular region to the left of the meridian in the polarization degree image, and calculating the sum of polarization degrees within the circular region to the right of the meridian in the polarization degree image. The expression is as follows:

[0050] (12),

[0051] Where Dop is the degree of polarization, i is the pixel index, and the summation range is limited to the circular region. and the circular region All pixels inside, The circular region The sum of internal polarization degrees, The circular region The sum of internal polarization degrees;

[0052] Compare the sum of the degrees of polarization in the circular region to the left of the meridian with the sum of the degrees of polarization in the circular region to the right of the meridian;

[0053] The sun's orientation points from the region with a larger sum of polarization degrees to the region with a smaller sum of polarization degrees.

[0054] Optionally, the step of calculating the navigation baseline based on the astronomical almanac includes:

[0055] The angle difference between the solar azimuth and north is obtained using the following expression:

[0056] (13),

[0057] in, The azimuth of the sun. The solar altitude angle, Indicates the geographical latitude of the observation point. Indicates the solar declination angle. The solar hour angle (c) is determined by the geographical latitude (L) of the observation point and the solar declination angle. and solar hour angle An intermediate variable c, formed by combining the elements, is used to reflect the adjustment effect of the geometric relationship between the Earth's rotation, the location of the observation point, and the position of the sun on the calculation of the solar azimuth angle.

[0058] The present invention has the following beneficial effects:

[0059] This invention provides a heading angle deambiguation method based on polarization angle gradient variation. By analyzing the gradient characteristics of the spatial distribution of polarization angles, it extracts the asymmetric information of polarization angle variation and combines it with gradient direction consistency checks within a dynamic window to achieve rapid heading angle deambiguity. This invention utilizes polarization degree gradient vectors at different positions on the same polarization angle to construct local constraints, directly calculating a unique heading angle through the geometric relationship between the gradient direction and the heading angle. Simultaneously, it introduces short-time angular velocity information from the inertial measurement unit to dynamically compensate for gradient changes, significantly improving robustness and real-time performance. Attached Figure Description

[0060] Figure 1 This is a schematic diagram of the celestial coordinate system provided in Embodiment 1 of the present invention.

[0061] Figure 2 This is a schematic diagram of the sky polarization light distribution pattern provided in Embodiment 1 of the present invention.

[0062] Figure 3 This is a schematic diagram of the polarization degree and polarization angle distribution at azimuth angle 0° provided in Embodiment 1 of the present invention.

[0063] Figure 4 This is a schematic diagram of the polarization degree and polarization angle distribution at an azimuth angle of 30° provided in Embodiment 1 of the present invention.

[0064] Figure 5 This is a schematic diagram of the polarization degree and polarization angle distribution at an azimuth angle of 60° provided in Embodiment 1 of the present invention.

[0065] Figure 6 This is a schematic diagram of the polarization degree and polarization angle distribution at an azimuth angle of 90° provided in Embodiment 1 of the present invention.

[0066] Figure 7 This is a schematic diagram of the degree of polarization and polarization angle distribution at an elevation angle of 0° provided in Embodiment 1 of the present invention.

[0067] Figure 8 This is a schematic diagram of the degree of polarization and polarization angle distribution at an elevation angle of 30° provided in Embodiment 1 of the present invention.

[0068] Figure 9 This is a schematic diagram of the polarization degree and polarization angle distribution at a height angle of 60° provided in Embodiment 1 of the present invention.

[0069] Figure 10 This is a schematic diagram of the degree of polarization and polarization angle distribution at a height angle of 90° provided in Embodiment 1 of the present invention.

[0070] Figure 11 The original image and polarization angle image are provided in Embodiment 1 of the present invention.

[0071] Figure 12 This is a schematic diagram of the polarization angle image coordinate system transformation provided in Embodiment 1 of the present invention.

[0072] Figure 13 The image shows the polarization angle after coordinate redrawing, as provided in Embodiment 1 of the present invention.

[0073] Figure 14 This is a schematic diagram of the azimuth angle calculation process provided in Embodiment 1 of the present invention.

[0074] Figure 15 This is a schematic diagram of polarization ambiguity provided in Embodiment 1 of the present invention.

[0075] Figure 16 This is a schematic diagram of the deambiguity process provided in Embodiment 1 of the present invention.

[0076] Figure 17 This is a flowchart of the heading angle calculation provided in Embodiment 1 of the present invention. Detailed Implementation

[0077] To enable those skilled in the art to better understand the technical solution of the present invention, the heading angle deambiguation method based on polarization angle gradient change provided by the present invention will be described in detail below with reference to the accompanying drawings.

[0078] Example 1

[0079] This embodiment provides a heading angle deambiguation method based on polarization angle gradient variation, including:

[0080] The polarization camera is calibrated to obtain the focal length and principal point parameter values ​​of the polarization camera;

[0081] Use a polarization camera to acquire sky images;

[0082] Calculate the degree of polarization and the angle of polarization of each point in the sky image;

[0083] A sky polarization pattern distribution image is obtained based on the polarization degree information and the polarization angle information;

[0084] Based on the parameter values ​​of the polarization camera, the coordinate system of the polarization camera is obtained;

[0085] Calculate the polarization angle in the polarization camera coordinate system according to Stokes' theorem;

[0086] The polarization angle in the polarization camera coordinate system is transformed into the polarization angle in the meridional coordinate system according to the coordinate system redrawing method;

[0087] Obtain the polarization angle image after coordinate system redrawing;

[0088] The polarization angle along the meridian is obtained from the polarization angle image.

[0089] Extract the spatial gradient distribution of the polarization degree information and construct a polarization degree gradient vector field;

[0090] In the meridian coordinate system, a dynamic window is divided, and the consistency of the polarization gradient direction within the window is checked to select the effective gradient region.

[0091] Eliminate the 180° polarization angle ambiguity based on gradient direction asymmetry to generate a deambiguous polarization angle distribution image;

[0092] Short-time angular velocity data from an inertial measurement unit (IMU) is introduced to compensate for gradient direction drift within a dynamic window;

[0093] The solar azimuth angle is calculated based on the unambiguous polarization angle distribution.

[0094] Calculate navigation baselines based on astronomical yearbooks;

[0095] A unique heading angle is output by the geometric constraint relationship between the gradient vector and the heading angle.

[0096] Optionally, the parameters of the polarization camera include focal length and principal point.

[0097] Optionally, polarization cameras can be used to obtain images of polarized light at different polarization angles. To comprehensively describe the polarization state of the light, it is necessary to analyze the intensity of the light in different polarization directions. This information is typically represented by four Stokes parameters, written in matrix form as follows:

[0098] ;

[0099] Where S0 represents the total incident light intensity, S1 represents the intensity difference between the x and y components, S2 represents the intensity difference between the +45° and -45° polarization components, and S3 represents the intensity difference between the left-hand and right-hand circular polarization components. It's worth noting that these parameters are usually expressed as average light intensity values ​​over a period of time. Under natural light conditions, S3 = 0, therefore the direction of sunlight is... The intensity of light after transmission can be obtained using the following formula:

[0100] ,

[0101] When the polarizer angle is 0°, 45°, 90°, and 135°, the Stokes component... , , It can be obtained using the following formula:

[0102] ,

[0103] The degree of polarization and polarization angle of the incident light can be obtained using the following formulas:

[0104] Optionally, the expression for the polarization angle in the polarization camera coordinate system is as follows:

[0105] (8),

[0106] Where DoLP is the degree of linear polarization, DoP is the degree of polarization, Aop is the polarization angle, and S1 and S2 represent the light intensities of linearly polarized light in two mutually perpendicular polarization directions, respectively. .

[0107] Optionally, the polarization angle is the angle between the direction of electric vector vibration and the reference coordinate system of the polarization camera.

[0108] Optionally, the step of transforming the polarization angle in the polarization camera coordinate system to the polarization angle in the meridional coordinate system according to the coordinate system redrawing method includes:

[0109] Obtain the polarization scattering azimuth angle of observation point M in the polarization camera coordinate system. The polarization scattering azimuth angle of observation point N in the polarization camera coordinate system is obtained. ;

[0110] The polarization angles of observation points M and N in the meridian coordinate system are obtained through geometric transformation, as expressed below:

[0111] (9),

[0112] in, Let M be the polarization angle of the observation point in the meridian coordinate system. Let N be the polarization angle of the observation point in the meridian coordinate system. Let M be the polarization angle of the observation point in the polarization camera coordinate system. Let N be the polarization angle of the observation point in the coordinate system of the polarization camera.

[0113] Optional, also includes:

[0114] The polarization angle image is binarized based on the numerical range of the polarization angle along the meridian, as shown in the following expression:

[0115] (10),

[0116] Where Aop is the polarization angle, and Aop_Gray is the degree of polarization of the polarized grayscale image after binarization, with a value of only 0 or 1.

[0117] Optional, also includes:

[0118] The least squares method was used to fit a straight line to the binarized polarization angle image;

[0119] The solar azimuth angle is obtained when the sum of the squares of the residuals at all observation points is minimized.

[0120] Optional, also includes:

[0121] Define circular regions on both sides of the meridian. and This is used to statistically analyze local polarization characteristics, calculating the sum of polarization degrees within the circular region to the left of the meridian in the polarization degree image, and calculating the sum of polarization degrees within the circular region to the right of the meridian in the polarization degree image. The expression is as follows:

[0122] (12),

[0123] Where Dop is the degree of polarization, i is the pixel index, and the summation range is limited to the left and right circular regions. and All pixels inside, and That is, the sum of the degrees of polarization within the two circular ranges on the left and right;

[0124] Compare the sum of the degrees of polarization in the circular region to the left of the meridian with the sum of the degrees of polarization in the circular region to the right of the meridian;

[0125] The sun's orientation points from the region with a larger sum of polarization degrees to the region with a smaller sum of polarization degrees.

[0126] Optionally, the step of calculating the navigation baseline based on the astronomical almanac includes:

[0127] The angle difference between the solar azimuth and north is obtained using the following expression:

[0128] (13),

[0129] in, The azimuth of the sun. The solar altitude angle, Indicates the geographical latitude of the observation point. Indicates the solar declination angle. The solar hour angle (c) is determined by the geographical latitude (L) of the observation point and the solar declination angle. and solar hour angle The intermediate variables formed by the combination reflect the adjustment effect of the geometric relationship between the Earth's rotation, the observation point's position, and the Sun's position on the azimuth calculation.

[0130] To address the aforementioned issues, this embodiment improves upon the existing meridian fitting method for heading calculation by proposing a fast heading angle deambiguation method based on polarization degree gradient changes. By analyzing the gradient characteristics of the spatial distribution of polarization angles, the asymmetric information of polarization angle changes is extracted, and combined with gradient direction consistency checks within a dynamic window, fast heading angle deambiguity is achieved. Specifically, local constraints are constructed using polarization degree gradient vectors at different positions on the same polarization angle. The unique heading angle is directly calculated through the geometric relationship between the gradient direction and the heading angle. Simultaneously, short-time angular velocity information from the Inertial Measurement Unit (IMU) is introduced to dynamically compensate for gradient changes, significantly improving robustness and real-time performance in complex environments.

[0131] 1. Atmospheric polarization model analysis and polarization mode simulation

[0132] First, set the observer's location as the origin of the coordinate system, with the east and north directions as the coordinate axes, and establish a navigation coordinate system.

[0133] Figure 1 This is a schematic diagram of the celestial coordinate system provided in Embodiment 1 of the present invention. Figure 1 As shown, S represents the intersection of the sphere with the direction of the sun, and P represents the intersection of a certain observation direction with the sphere. The angle between the line connecting the observer's position and the sun's position and the horizontal plane is defined as the solar altitude angle. The angle between the line connecting the observer's position and the sun's position, projected onto the horizontal plane, and the positive X-axis is defined as the solar azimuth angle. Similarly, the definitions of the altitude angle and azimuth angle of the observation point can be obtained.

[0134] Next, based on the Rayleigh scattering model, the degree of polarization and polarization angle at any observation point P are determined. According to equation (1), the degree of polarization at the observation point needs to be obtained. This can be achieved by solving the cosine formula of the line connecting the observer's location and the sun's location, and the line connecting the observer and the observation point.

[0135] (1),

[0136] The scattering angle can be obtained through calculation. :

[0137] (2),

[0138] The degree of polarization at any observation point P can be obtained:

[0139] (3),

[0140] According to the Rayleigh scattering model, The direction of incident light scattering at point A is perpendicular to the plane defined by the observer, the observation point, and the sun. Therefore, the polarization angle of the incident light... Defined as Point scattering direction and the direction of scattering through N on the sphere The angle formed by the tangent of the circular arc. The coordinate system is defined by the origin of the observation point's coordinates, with the perpendicular and parallel directions on the sphere as the coordinate circumference. Navigation Horizontal Coordinate System to the coordinate system of the observation point The transition matrix is ​​defined as:

[0141] (4),

[0142] In the original coordinate system middle, The direction of point light vector vibration can be expressed as:

[0143] (5),

[0144] It can be obtained in coordinate system The vector of the vibration direction of the point light vector:

[0145] (6),

[0146] Based on this, we can obtain Polarization azimuth at point It can be represented as:

[0147] (7),

[0148] Figure 2This is a schematic diagram of the sky polarization light distribution pattern provided in Embodiment 1 of the present invention. Figure 2 As shown, an atmospheric polarization model can be obtained by analyzing the distribution of polarization degree and polarization angle at sky observation points. In the figure, the thickness of the dashed line represents the magnitude of polarization degree, and the tangent of the dashed line represents the polarization angle and direction.

[0149] The atmospheric polarization pattern appears as a series of concentric circles with the sun at the center. As shown in the figure above, the dashed line is thickest when the observation point is perpendicular to the sun's direction, representing the greatest degree of polarization. Furthermore, the degree of polarization gradually decreases as the angle between the observation point and the sun's direction decreases.

[0150] Furthermore, the figure also illustrates the axial symmetry of the polarization distribution pattern. The axis of symmetry is the solar meridian, which lies on the arc of the sphere connecting the sun's position and the zenith. Along the meridian, the polarization azimuth angle is 90 degrees at all locations. Moreover, the degree of polarization and the polarization angle are symmetrically distributed in the region surrounding the meridian. Throughout the day, as the sun rises in the east, reaches its southern high, and then sets in the west, its altitude angle gradually increases and then decreases, while the azimuth angle also rotates over time. For the zenith region, the degree of polarization first increases, then decreases, and then gradually increases again, while the polarization angle continuously rotates.

[0151] To more intuitively observe the distribution of polarization degree and polarization angle of sky polarized light under different solar altitude and azimuth angles, and to facilitate the analysis and simulation of heading angle calculation algorithms, a three-dimensional simulation of the atmospheric polarization pattern was performed using Matlab software, and pseudo-color was used to represent the distribution of polarization degree and polarization angle. With a uniform solar altitude angle of 90°, the polarization pattern distribution was observed at solar azimuth angles of 0°, 30°, 60°, and 90°, yielding the following results: Figure 3 , Figure 4 , Figure 5 and Figure 6 , Figure 3-6 The left side shows the distribution of polarization degree, and the right side shows the distribution of polarization angle.

[0152] The 3D simulation model vividly demonstrates the overall layout of the sky's polarization pattern. The distribution of polarization degree presents a clear concentric circle shape, with the polarization degree reaching its maximum at a position 90° to the sun, while the polarization degree on both sides of this position is symmetrically distributed. As the viewing angle deviates from this position, the polarization degree gradually decreases. In the distribution of polarization angle, the solar meridian exists, with the figures on both sides being symmetrical in shape and anti-symmetrical in color. When focusing on the zenith region, a regular periodic change in the polarization angle around the zenith can be observed, with a period of 180°.

[0153] By observing these charts, it can be seen that changes in the solar azimuth angle do not alter the overall shape of the sky polarization pattern; rather, they cause the symmetrical points of the polarization pattern to rotate with the rotation of the solar azimuth angle. With the solar azimuth angle set to 0°, the polarization pattern distributions were observed at solar altitude angles of 0°, 30°, 60°, and 90°, yielding the following results: Figure 7 , Figure 8 , Figure 9 and Figure 10 , Figure 7-10 The left side shows the distribution of polarization degree, and the right side shows the distribution of polarization angle.

[0154] The Sun's position in the sky directly affects the distribution of polarization patterns in scattered light from the atmosphere, particularly the shape of the polarization angle. As the Sun's position rises, the points of symmetry in the distribution also rise accordingly. Although the shape of the polarization angle remains constant, the positions of these points of symmetry rotate continuously with changes in the Sun's azimuth angle.

[0155] However, even though the shape of the polarization angle distribution pattern is constantly changing, there is always a symmetrical line passing through the direction of the sun, namely the solar meridian. When the solar altitude angle reaches 90°, an extreme case occurs, where the polarization angle has only two distributions: ±90°, located on both sides of the solar meridian.

[0156] 2. A heading acquisition method based on an improved fitting meridian method

[0157] The distribution of atmospheric polarization patterns is closely related to the solar elevation and azimuth angles. By observing and analyzing real-time atmospheric polarization patterns, important information about the solar elevation and azimuth angles can be inferred. Under ideal atmospheric conditions based on Rayleigh scattering theory, the solar meridian is a straight line in the sky polarization distribution model. Using Stokes' theorem... The polarization angle of the image obtained by the polarization camera is calculated to obtain... Figure 11 . Figure 11 The original image is shown on the left, and the polarization angle image is shown on the right.

[0158] In the sky polarization mode, the axis of symmetry is the location of the solar meridian, and the polarization angle around the meridian is ±90°. However, the polarization angle image is not centrally symmetrical. The reason can be found by analyzing the polarization angle calculation formula (8), because the polarization angle calculated in the formula is the angle between the electric vector vibration direction and the camera's reference coordinate system.

[0159] (8),

[0160] Where DoLP is the degree of linear polarization, DoP is the degree of polarization, Aop is the polarization angle, and S1 and S2 represent the light intensities of two sets of linearly polarized light that are perpendicular to each other in the polarization direction. .

[0161] Therefore, there is no central symmetry, making meridian extraction difficult and unable to accurately extract the solar azimuth angle. Therefore, this embodiment proposes a coordinate system redrawing method for polarization angle image coordinate system transformation, such as... Figure 12 As shown.

[0162] First, for ease of analysis, a polarization observation plane coordinate system Qxy is established, with the origin at the camera's optical center O, where OS is the solar meridian. M is an observation point in the sky, and N is another observation point symmetrical about the meridian. According to the atmospheric scattering model, the polarization scattering direction of the incident light at any observation point in the sky must be perpendicular to the line connecting that point and the sun's position; that is, SM is perpendicular to ML, and SN is perpendicular to NL. Using the horizontal axis of the coordinate system as the polarization reference value, the polarization angle between the two points can be obtained as follows: and Obviously, these two angles are not equal. Therefore, there is no polarization angle symmetry about the solar meridian in the polarization observation plane coordinate system, and the solar meridian cannot be obtained by fitting the initial polarization angle image obtained by the polarization camera.

[0163] Therefore, in this embodiment, the line connecting the point projected by the sun into the coordinate system to the origin and the direction perpendicular to this line are used as coordinate axes to establish the solar meridian coordinate system Ost. The direction of the solar meridian is defined as the horizontal axis. From the geometric relationship in the figure, it can be seen that the polarization angle between the two points is... and The two angles are equal in magnitude and remain symmetrical about the meridian regardless of the rotation of the observation plane, conforming to the atmospheric polarization model. However, directly solving for these two points is difficult; geometric analysis can be used to obtain... , Furthermore, through geometric transformation, it can be obtained that and The calculation formula is as follows:

[0164] (9),

[0165] In the formula, and Let M and N be the polarization scattering azimuth angles in the initial coordinate system of the observation points, which can be obtained from the pixel coordinates of the two points. and The observation point M and N The polarization angle in the initial observation plane.

[0166] The polarization angle in the camera coordinate system is calculated using Stokes' theorem, and then the polarization angle in the meridional coordinate system is obtained through coordinate transformation. The image is then redrawn as follows: Figure 13 As shown, the meridian features are obvious in the image. The polarization angle along the meridian is ±90°, and the magnitude of the polarization angle is centrally symmetrical about the meridian. Redrawing the image by rotating the camera 90° three times consecutively shows that the meridian rotation angle is close to the camera rotation angle. Therefore, a preliminary relationship between the meridian azimuth and the camera rotation angle can be observed.

[0167] After obtaining the polarization angle image, it is binarized according to the numerical range of the polarization angle along the meridian, as shown in the following formula:

[0168] (10),

[0169] To obtain the solar azimuth angle from a binarized meridian image, a least-squares method is needed for line fitting. Least-squares fitting is a mathematical optimization technique used to find the best-fit curve for a set of data. It achieves this by minimizing the sum of squared errors, where the error refers to the vertical distance between each data point and the curve. The goal of least-squares fitting is to determine the model parameters that minimize the difference between the model's predictions and the actual observations (i.e., the sum of squared residuals).

[0170] Suppose we have a set of observation data points I hope to achieve this through a function. Fit these data, where A vector representing the model parameters. For each data point Its residual r i Defined as observation value Compared with the predicted values ​​calculated by the model The difference between the two, i.e. .

[0171] After calculating the residuals for each observation point, construct the overall error function:

[0172] (11),

[0173] The goal of the least squares method is to minimize the sum of squares of the residuals at all observation points. This is achieved by considering each parameter in S... Taking the partial derivatives and setting them to zero yields a set of linear equations. Solving these equations provides the optimal estimates of the model parameters. Where k is the slope, its value in the camera coordinate system can be obtained by performing an arctangent transformation on it. , Figure 14 To solve the azimuth angle The process.

[0174] 3. Heading Angle Deambiguation Method Based on Polarization Degree Variation

[0175] The heading angle calculated using the above method ranges from 0° to 180°. However, in practical applications, the heading angle of the carrier varies from 0° to 360°. Therefore, when the carrier's rotation angle exceeds 180°, the heading angle calculated by the above method will jump from 180° to 0°, which is called the heading angle ambiguity problem. This causes significant interference to polarized light navigation applications. This embodiment proposes a heading angle deambiguation technique based on polarization gradient variation, achieving full observability of the carrier's heading and ensuring the feasibility of practical navigation applications.

[0176] First, analysis of atmospheric polarization patterns reveals that the polarization angle is symmetrical about the solar meridian, which represents the sun's location. The degree of polarization varies along the meridian; the closer to the sun, the closer the degree of polarization is to 0, and the farther away, the closer it is to 1. Therefore, this embodiment determines the sun's absolute location by comparing the magnitudes of polarization along the meridian, thus eliminating ambiguity. Figure 15 It can be seen that the angle of polarization (AOP) images of the camera at the initial position and the image rotated 180° are almost identical, meaning that the two cases will calculate the same solar azimuth angle, but the actual solar azimuth angle differs by 180°. By observing the degree of polarization (DOP) images, it can be seen that the DOP values ​​of the two are completely opposite.

[0177] Figure 16 The process of deambiguity is demonstrated by calculating the sum of the degrees of polarization in the circular regions to the left and right of the meridian in the Dop diagram using the following formula:

[0178] (12),

[0179] Comparing the magnitudes of the two sums, the sun's azimuth is the direction from the region with the larger sum of polarization degrees to the region with the smaller sum. The yaw angle calculation result is from... Become Thus, the problem of ambiguity regarding the sun's azimuth angle has been resolved.

[0180] However, information about the sun's azimuth alone can only determine the relative position between the sun and the observer. In navigation practice, to implement accurate navigation, in addition to determining the sun's azimuth mentioned above, it is necessary to further use tools such as astronomical almanacs to calculate the navigation baseline, which is the angle between the sun's meridian and a fixed direction (such as north). Figure 17This is a flowchart of the heading angle calculation provided in Embodiment 1 of the present invention. By combining these two pieces of information, the specific location of the designated direction (e.g., north) at the observation point can be clearly specified, thereby achieving the navigation goal. The use of astronomical almanacs involves a detailed set of calculation methods. In short, the formula for calculating the angle difference between the solar azimuth and north is (13):

[0181] (13),

[0182] in, The azimuth of the sun. The solar altitude angle, Indicates the geographical latitude of the observation point. Indicates the solar declination angle. The solar hour angle (c) is determined by the geographical latitude (L) of the observation point and the solar declination angle. and solar hour angle An intermediate variable c, formed by combining the elements, is used to reflect the adjustment effect of the geometric relationship between the Earth's rotation, the location of the observation point, and the position of the sun on the calculation of the solar azimuth angle.

[0183] (1) Solve the problem of ambiguous heading angle and improve navigation accuracy.

[0184] By proposing a heading angle deambiguation technique based on polarization degree gradient variation, this embodiment effectively solves the inherent 180-degree heading angle ambiguity problem in traditional single-point polarization angle calculation, extending the original 180° observability to a 360° all-sky observability, ensuring the continuity and accuracy of heading angle calculation. Compared with existing technical solutions, this embodiment directly calculates a unique heading angle through the asymmetric characteristics of polarization angle gradient, greatly improving the accuracy of the navigation system. Especially in practical applications where the heading angle variation exceeds 180°, it can continuously and stably calculate the correct heading angle information. The proposed heading angle deambiguation method based on polarization degree variation effectively resolved the ambiguity during experiments.

[0185] (2) Accurate determination of solar azimuth based on polarized light mode

[0186] This embodiment proposes a solar azimuth angle determination method based on polarization degree by deeply analyzing the solar polarization mode and its propagation characteristics in the atmosphere, combined with the variation law of solar meridian polarization degree. Compared with existing schemes, the improved scheme more accurately determines the absolute azimuth of the sun, eliminates errors in traditional methods, improves the calculation accuracy of heading angle, and ensures higher stability and reliability.

[0187] (3) Improve the robustness and real-time performance of the system

[0188] Compared to existing polarization mode library traversal optimization methods, this embodiment improves computational efficiency by more than 3 times, enabling millisecond-level real-time heading angle output to meet the navigation requirements of high-speed moving vehicles. At the same time, it does not rely on external solar position information, significantly enhancing system independence and deployment flexibility.

[0189] Meanwhile, by incorporating polarization analysis and improving the solar meridian fitting method, this embodiment demonstrates stronger robustness in varying lighting environments. In multiple outdoor experiments, compared to existing technologies' dependence on ideal weather and static scenarios, this embodiment maintains heading angle calculation accuracy even under conditions of localized cloud cover and dynamic interference (such as rapid UAV maneuvers), exhibiting a 50% improvement in anti-interference capability and adaptability to diverse environments such as urban canyons and underground spaces. Therefore, this embodiment's ability to rapidly calculate accurate heading angles ensures the real-time performance and reliability of the navigation system, meeting the dual requirements of real-time performance and accuracy in practical navigation applications.

[0190] To address the ambiguity issue of heading angles, besides the deambiguation method based on polarization angle gradient changes proposed in this embodiment, there are theoretically several alternative technical paths that can achieve similar goals. However, these solutions differ significantly in practicality, cost, or performance. The following lists possible alternatives and analyzes their advantages and disadvantages:

[0191] (1) Fusion of multispectral polarization sensors

[0192] By utilizing the differences in polarization angles across different spectral bands (such as visible light and near-infrared), the 180-degree ambiguity can be eliminated through the asymmetry of multispectral polarization distribution. However, this requires the integration of a multispectral polarization sensor, significantly increasing hardware costs, and the multispectral data fusion algorithm is complex, making real-time performance difficult to guarantee.

[0193] (2) Loosely coupled assistance of inertial navigation system (INS)

[0194] This method uses inertial navigation system (INS) trajectory estimation to provide short-term heading change trends, and combines polarization angle measurement results for direction filtering (e.g., eliminating reverse headings by angular velocity integration). This method directly utilizes existing INS hardware resources, resulting in low system modification costs. However, the accumulated error of the INS increases over time, and the auxiliary effect becomes ineffective after long-term operation. Furthermore, it cannot independently resolve ambiguity issues in static scenarios (such as when the carrier is stationary).

[0195] (3) Geomagnetic-assisted correction

[0196] This method combines an absolute heading reference provided by a geomagnetic sensor, eliminating ambiguity by using the difference between magnetic heading and polarization heading. The logic is simple, and it works even in static scenarios, but its drawback is that it is susceptible to environmental interference (such as from metal structures and electromagnetic devices).

[0197] Traditional heading angle calculations using the sun meridian method inherently suffer from 180° ambiguity, especially when the vehicle undergoes large-angle rotations (exceeding 180°), where abrupt changes in the heading angle calculation can occur, affecting the stability and accuracy of the navigation system. To address this critical issue, this embodiment innovatively proposes a heading angle disambiguation method based on polarization gradient variation. By introducing a polarization gradient algorithm, it effectively corrects the heading angle calculation, ensuring the continuity of the solution results and significantly improving the system's accuracy and reliability.

[0198] (1) Combining atmospheric polarization distribution characteristics with the solar meridian method

[0199] This embodiment fully utilizes the distribution patterns of atmospheric polarization modes and combines them with the solar meridian method to improve the accuracy of heading angle calculation. First, through in-depth modeling and analysis of atmospheric polarization modes, the mechanism by which changes in the sun's position affect the distribution of polarized light is revealed. Then, a heading angle calculation method based on the symmetry of the solar meridian is proposed. Using this method, the solar azimuth angle can be calculated with high precision, thus providing solid theoretical support for subsequent heading angle disambiguation and accurate calculation. This approach not only improves the accuracy of heading angle calculation but also lays the foundation for the application of polarized light navigation systems in complex environments.

[0200] (2) Polarization angle redrawing and meridian extraction techniques

[0201] To achieve accurate heading angle calculation, this embodiment employs the Stokes parameter calculation method to process polarized light images, redraws the polarization angle, and combines coordinate system transformation to obtain a more stable and accurate meridian feature map. By acquiring image information in real time using a polarized light camera, this embodiment further extracts and redraws the polarization angle, making the calculated meridian information clearer and thus enhancing the robustness of the heading angle calculation. Furthermore, this method can effectively establish the correlation between the trend of polarization angle variation along the meridian direction and the solar azimuth angle, providing reliable data support for accurate heading angle calculation. The application of this technology enables the polarized light navigation system to maintain high accuracy under different lighting conditions and possess strong environmental adaptability.

[0202] (3) Heading angle deambiguation technique based on polarization degree gradient

[0203] To address the 180° ambiguity issue in the traditional solar meridian method for calculating heading angles, this embodiment innovatively improves upon the fitting meridian method by proposing a rapid heading angle deambiguation method based on polarization degree gradient changes. This method analyzes the polarization degree changes within the meridian region and rapidly compares the magnitudes of the polarization degree gradients, effectively eliminating the ambiguity caused by the symmetry of the solar azimuth angle. Compared to traditional deambiguation methods, this embodiment's solution requires no additional reference information and can achieve continuous heading angle calculation across the entire 360° omnidirectional range. It not only eliminates the heading angle jump problem but also significantly improves the accuracy and stability of heading angle calculations, enabling polarized light navigation systems to operate reliably in a wider range of application scenarios.

[0204] In summary, this embodiment achieves a breakthrough improvement over the traditional solar meridian method by deeply exploring the distribution characteristics of atmospheric polarized light, optimizing polarization angle redrawing and meridian extraction techniques, and innovatively applying a heading angle deambiguation method based on polarization degree gradient. This provides a more accurate and robust solution for autonomous polarized light navigation.

[0205] It is understood that the above embodiments are merely exemplary implementations used to illustrate the principles of the present invention, and the present invention is not limited thereto. For those skilled in the art, various modifications and improvements can be made without departing from the spirit and essence of the present invention, and these modifications and improvements are also considered to be within the scope of protection of the present invention.

Claims

1. A heading angle deambiguation method based on polarization angle gradient variation, characterized in that, include: The polarization camera is calibrated to obtain the focal length and principal point parameter values ​​of the polarization camera; The polarization camera is used to acquire sky images, and the degree of polarization and polarization angle information of each point are calculated. A sky polarization pattern distribution image is obtained based on the polarization degree information and polarization angle information; Based on the parameter values ​​of the polarization camera, the coordinate system of the polarization camera is obtained; Calculate the polarization angle in the polarization camera coordinate system according to Stokes' theorem; The polarization angle in the polarization camera coordinate system is transformed into the polarization angle in the meridional coordinate system according to the coordinate system redrawing method; Obtain the polarization angle image after coordinate system redrawing, and get the polarization angle on the meridian; Extract the spatial gradient distribution of the polarization degree information and construct a polarization degree gradient vector field; A dynamic window is divided in the meridian coordinate system, and the consistency of the polarization gradient direction within the dynamic window is checked to select a preset gradient region. Eliminate the 180° polarization angle ambiguity based on gradient direction asymmetry to generate a deambiguous polarization angle distribution image; Short-time angular velocity data from an inertial measurement unit are used to compensate for gradient direction drift within the dynamic window. Based on the deambigued polarization angle distribution image, the solar azimuth angle is calculated; Calculate navigation baselines based on astronomical yearbooks; A unique heading angle is output by the geometric constraint relationship between the gradient vector and the heading angle.

2. The heading angle deambiguation method based on polarization angle gradient change according to claim 1, characterized in that, Also includes: The polarization state of light and the intensity of light in different polarization directions are described using four Stokes parameters, expressed in matrix form as follows: ； Where S0 represents the total incident light intensity, S1 represents the light intensity difference between the x and y components, S2 represents the light intensity difference between the +45° and -45° polarization components, and S3 represents the light intensity difference between the left-handed and right-handed circular polarization components.

3. The heading angle deambiguation method based on polarization angle gradient change according to claim 2, characterized in that, Under natural light conditions, S3=0, and the direction of sunlight is... At that time, the expression for the light intensity after transmission is as follows: ， When the polarizer angle is taken as 0°, 45°, 90° and 135° respectively, the expression is as follows: ， Where S0, S1, and S2 are Stokes parameters.

4. The heading angle deambiguation method based on polarization angle gradient change according to claim 3, characterized in that, The expression for the polarization angle in the polarization camera coordinate system is as follows: （8）， Where DoLP is the degree of linear polarization, Dop is the degree of polarization, Aop is the polarization angle, and S1 and S2 represent the light intensities of two sets of linearly polarized light that are perpendicular to each other in the polarization direction.

5. The heading angle deambiguation method based on polarization angle gradient change according to claim 4, characterized in that, The polarization angle is the angle between the direction of electric vector vibration and the reference coordinate system of the polarization camera.

6. The heading angle deambiguation method based on polarization angle gradient change according to claim 5, characterized in that, The step of transforming the polarization angle in the polarization camera coordinate system to the polarization angle in the meridional coordinate system according to the coordinate system redrawing method includes: Obtain the polarization scattering azimuth angle of observation point M in the polarization camera coordinate system. The polarization scattering azimuth angle of observation point N in the polarization camera coordinate system is obtained. ; The polarization angles of observation points M and N in the meridian coordinate system are obtained through geometric transformation, as expressed below: （9）， in, Let M be the polarization angle of the observation point in the meridian coordinate system. Let N be the polarization angle of the observation point in the meridian coordinate system. Let M be the polarization angle of the observation point in the polarization camera coordinate system. Let N be the polarization angle of the observation point in the coordinate system of the polarization camera.

7. The heading angle deambiguation method based on polarization angle gradient change according to claim 6, characterized in that, Also includes: The polarization angle image is binarized based on the numerical range of the polarization angle along the meridian, as shown in the following expression: （10）， Where Aop is the polarization angle, and Aop_Gray is the degree of polarization of the polarized grayscale image after binarization, with a value of 0 or 1.

8. The heading angle deambiguation method based on polarization angle gradient change according to claim 7, characterized in that, Also includes: The least squares method was used to perform linear fitting on the binarized polarization angle image; The solar azimuth angle is obtained when the sum of the squares of the residuals at all observation points is minimized.

9. The heading angle deambiguation method based on polarization angle gradient change according to claim 8, characterized in that, Also includes: Define a circular area to the left of the meridian. Define a circular area to the right of the meridian. This is used to statistically analyze local polarization characteristics, calculating the sum of polarization degrees within the circular region to the left of the meridian in the polarization degree image, and calculating the sum of polarization degrees within the circular region to the right of the meridian in the polarization degree image. The expression is as follows: （12）， Where Dop is the degree of polarization, i is the pixel index, and the summation range is limited to the circular region. and the circular region All pixels inside, The circular region The sum of internal polarization degrees, The circular region The sum of internal polarization degrees; Compare the sum of the degrees of polarization in the circular region to the left of the meridian with the sum of the degrees of polarization in the circular region to the right of the meridian; The sun's orientation points from the region with a larger sum of polarization degrees to the region with a smaller sum of polarization degrees.

10. The heading angle deambiguation method based on polarization angle gradient change according to claim 9, characterized in that, The steps for calculating the navigation baseline based on the astronomical almanac include: The expression for obtaining the solar azimuth angle is as follows: （13）， in, The azimuth of the sun. The solar altitude angle, Indicates the geographical latitude of the observation point. Indicates the solar declination angle. The solar hour angle (c) is determined by the geographical latitude (L) of the observation point and the solar declination angle. and solar hour angle An intermediate variable c, formed by combining the elements, is used to reflect the adjustment effect of the geometric relationship between the Earth's rotation, the location of the observation point, and the position of the sun on the calculation of the solar azimuth angle.

Images