A navigation method suitable for field geological surveys

Abstract

This invention discloses a navigation method suitable for field geological surveys, specifically relating to the fields of geological exploration and navigation positioning technology. It involves collecting three-component geomagnetic field data to calculate rotational angular velocity and extracting travel cycles to calculate triaxial cross-correlation coefficients. When both rotational and correlation anomalies occur simultaneously, a geomagnetic anomaly zone is identified, and the mechanical and electronic compasses are forcibly shut down. The accurate magnetic azimuth and satellite coordinates before entering the anomaly zone are extracted as initial references. Inertial dead reckoning is initiated, and the heading change or absolute heading angle is obtained. Periodic drift corrections are performed on the dead reckoning to correct the navigation heading. This invention accurately identifies geomagnetic anomaly zones using dual anomaly detection based on geomagnetic field vector rotational angular velocity and triaxial correlation coefficients, forcibly shutting down compass output to avoid misleading information. Furthermore, by integrating inertial dead reckoning with visual or polarized light heading correction, a corrected navigation heading independent of the geomagnetic field is output, ensuring the accuracy and safety of field geological surveys.

Description

Technical Field

[0001] This invention relates to the field of geological exploration and navigation positioning technology, and more specifically, to a navigation method suitable for field geological surveys. Background Technology

[0002] In field geological surveys, investigators frequently need to traverse areas with strongly magnetic rock masses, such as magnetite deposits, basalt outcrops, and ophiolite tectonic zones. These areas contain abundant ferromagnetic minerals, causing significant distortions in the local geomagnetic field. The total geomagnetic field strength and magnetic declination can fluctuate dramatically within a range of tens or even several meters. Both mechanical and electronic compasses rely on measuring the direction of the surrounding geomagnetic field to determine geographic north. When the geomagnetic field is distorted, the azimuth given by the compass is no longer true geographic north, but a distorted and incorrect direction. Furthermore, this magnetic anomaly is usually imperceptible to the human senses. Investigators often only discover the problem after wandering off the planned route for an extended period, by which time they have already wasted considerable time and energy, and may even find themselves in dangerous situations due to repeated detours.

[0003] In existing field navigation technologies, global navigation satellite systems (GPS) are prone to signal loss or multipath effects in obstructed environments such as canyons and dense forests, failing to provide continuous and reliable positioning. While inertial navigation systems do not rely on external magnetic fields or satellite signals, they suffer from long-term integral drift, making it difficult to maintain pointing accuracy over long distances. Positioning methods based on geomagnetic fingerprint matching typically require pre-collection of geomagnetic reference maps for the entire area, which is difficult to implement in unknown or large-scale geological survey areas. Therefore, this invention proposes a navigation method suitable for field geological surveys to address the aforementioned problems. Summary of the Invention

[0004] To achieve the above objectives, the present invention provides the following technical solution: A navigation method suitable for field geological surveys includes the following steps: The three components of the geomagnetic field data output by the triaxial magnetometer are continuously collected. The spatial angle between the geomagnetic field vectors at two adjacent sampling times is divided by the sampling time interval to obtain the angular velocity of the geomagnetic field vector rotation. A complete walking cycle is selected as the phase analysis window. The cross-correlation coefficients between the output sequences of any two axes of the triaxial magnetometer within the phase analysis window are calculated to obtain a set of triaxial correlation coefficients containing three cross-correlation coefficients. Rotational anomalies are determined by using the rotational angular velocity of the geomagnetic field vector, and correlation anomalies are determined by using the triaxial correlation coefficient set. When both rotational and correlation anomalies are found in the judgment results, the current area is determined to be a geomagnetic anomaly zone. In response to the determination that the current area is a geomagnetic anomaly zone, the azimuth output of the mechanical compass and electronic compass is forcibly blocked, and an alarm is issued; Extract the accurate magnetic azimuth and the last valid satellite positioning coordinates recorded at the last moment before the geomagnetic anomaly zone was identified. Using the accurate magnetic azimuth as the initial heading reference, start the gyroscope and accelerometer in the inertial measurement unit to perform dead reckoning and obtain continuous relative heading angles. The camera unit is activated to acquire continuous surface images. The optical flow features between adjacent image frames are extracted using the visual odometry method to calculate the course change. Alternatively, the sky polarization sensor is activated to acquire the atmospheric polarization mode and calculate the absolute course angle. The course change or the absolute course angle is used to perform periodic drift correction on the continuous relative course angle, and a corrected navigation course that does not depend on the geomagnetic field is output.

[0005] In a preferred embodiment, performing rotation anomaly detection and related anomaly detection refers to: The rotational angular velocity of the geomagnetic field vector is compared with the first historical baseline, which is the moving average of the rotational angular velocity of the geomagnetic field vector over a continuous period before entering the current region. When the rotational angular velocity of the geomagnetic field vector exceeds a predetermined multiple of the first historical baseline, it is marked as a rotational anomaly. The triaxial correlation coefficient set is compared with the second historical baseline, which is the moving average of the triaxial correlation coefficient set over a continuous period before entering the current region. When the minimum value in the triaxial correlation coefficient set is lower than a predetermined proportion of the corresponding minimum value in the second historical baseline, it is marked as a correlation anomaly.

[0006] In a preferred embodiment, the method for extracting the complete walking cycle is as follows: The vertical acceleration value output by the accelerometer in the inertial measurement unit is continuously monitored. When the vertical acceleration value crosses the zero point twice in a row and the time interval between the two crossings is within the preset normal walking frequency range, the time interval between the two zero point crossings is determined as a complete walking cycle. When step frequency cannot be detected, a fixed duration is used as a complete walking cycle. The cross-correlation coefficient is calculated as follows: for any two axis output sequences, multiply the corresponding data points of the two sequences within the same phase analysis window, sum them, and then divide by the product of the square root of each sequence to obtain the cross-correlation coefficient.

[0007] In a preferred embodiment, dead reckoning is specifically implemented as follows: The accurate magnetic azimuth recorded at the last moment before entering the geomagnetic anomaly zone is used as the initial heading angle, and the last effective satellite positioning coordinates are used as the initial position coordinates. By integrating the three-axis angular velocity output by the gyroscope over time, the change in heading angle relative to the initial heading angle at each moment can be obtained. Gait detection is performed using the triaxial acceleration output from the accelerometer. When the vertical acceleration value crosses the zero point twice consecutively, it is determined as a complete step. The step length is estimated based on the composite amplitude of the acceleration, and the step length is proportional to the fourth root of the acceleration amplitude. Multiply the step size by the cosine and sine of the current heading angle and add them to the position coordinates to obtain continuous position calculation results; After each predetermined number of steps is calculated, the angular velocity drift of the integral is corrected once using the change in heading calculated by the visual odometry method or the sky polarization light sensor. The correction amount is the average difference between the gyroscope integral heading and the visual heading within the most recent predetermined number of steps.

[0008] In a preferred embodiment, the specific steps for extracting optical flow features between adjacent image frames using the visual odometry method are as follows: The continuous surface images acquired by the camera unit are converted into grayscale images, and fast corner features are extracted from each grayscale image; The pyramid optical flow method is used to track the pixel displacement of corner features between two adjacent frames; Median filtering is applied to the displacement of all successfully tracked corner pixels to remove outliers whose displacement vector magnitude exceeds a predetermined multiple of the median standard deviation. The remaining effective pixel displacements are converted into unit pixel movement vectors, and then the heading change of the carrier in the horizontal plane is deduced by using the intrinsic parameter matrix of the camera unit and the assumed height of the carrier. When the number of effective corner points is lower than the predetermined number, the visual odometry method is deemed to be ineffective, and the system automatically switches to the sky polarization light sensor for heading calculation.

[0009] In a preferred embodiment, the specific steps for the sky polarization sensor to acquire the atmospheric polarization mode and calculate the absolute heading angle are as follows: The atmospheric polarization degree and polarization azimuth angle in three different directions in the sky are measured using a time-division or pixel-division method. The three directions are selected from the zenith direction and the direction of the sun at predetermined angles on both sides. According to the Rayleigh scattering model, there is a geometric constraint relationship between the polarization azimuth angles measured in the three directions and the position of the sun. By solving the system of constraint equations, the azimuth angle of the sun in the horizontal coordinate system can be deduced. The solar azimuth angle is compared with the theoretical solar azimuth angle provided by the internal clock of the device. The difference between the two is the absolute heading angle between the current orientation of the device and geographic north. When a situation arises where it is impossible to measure three valid directions simultaneously due to obstruction, the most recently successfully calculated absolute heading angle is automatically saved and used as the benchmark for subsequent dead reckoning until a valid measurement is obtained again.

[0010] In a preferred embodiment, the real-time calculation of the geomagnetic field vector rotation angular velocity and the determination of rotation anomalies employ the following steps: A three-component geomagnetic field data sequence was established, with each sampling point containing the northward component, the eastward component, and the vertical component; For two geomagnetic field vectors at two adjacent sampling times, calculate their spatial angle. The spatial angle is equal to the dot product of the two vectors divided by the product of the magnitudes of the two vectors and then taking the inverse cosine. Divide the spatial angle by the sampling time interval to obtain the instantaneous rotational angular velocity; Calculate the average instantaneous rotational angular velocity within the sliding window, with the sliding window length set to a preset value; Calculate the coefficient of variation of the rotational angular velocity within the sliding window. The coefficient of variation is equal to the standard deviation of the angular velocity within the window divided by the average value. Calculate the average value of the coefficient of variation over a continuous period before entering the current region, and use it as the historical baseline coefficient of variation; When the coefficient of variation of the historical baseline is less than the preset variation threshold, it is marked as a rotation anomaly when the angular velocity of the geomagnetic field vector rotation exceeds a predetermined multiple of the first historical baseline. When the coefficient of variation of the historical baseline is greater than or equal to the preset variation threshold, and the angular velocity of the geomagnetic field vector rotation exceeds a predetermined multiple of the first historical baseline and the current coefficient of variation exceeds a predetermined multiple of the coefficient of variation of the historical baseline, it is marked as a rotation anomaly. When a predetermined number of consecutive sampling times are all marked as rotational anomalies, the final rotational anomaly marker is output.

[0011] In a preferred embodiment, the correlation anomaly determination of the triaxial correlation coefficient set adopts the following steps: The first step is to denote the northward component sequence, eastward component sequence, and vertical component sequence of the geomagnetic field data collected within a complete walking cycle as the first sequence, the second sequence, and the third sequence, respectively. The length of each sequence is equal to the number of sampling points within that walking cycle. The second step is to calculate the cross-correlation coefficient between the first and second sequences. The cross-correlation coefficient is equal to the product of the covariance of the two sequences and their respective standard deviations. Similarly, calculate the cross-correlation coefficients between the first and third sequences, and between the second and third sequences, to obtain the set of triaxial correlation coefficients. The third step is to calculate the statistical characteristics of the correlation coefficients of the number of complete walking cycles before entering the current area, and record the minimum value of the triaxial correlation coefficient set in each cycle to form a historical minimum value sequence. The fourth step is to calculate the moving median and moving interquartile range of the historical minimum value sequence; Fifth, calculate the minimum value in the set of triaxial correlation coefficients for the current walking cycle; Step 6: When the current minimum value is less than the sliding median minus 1.5 times the sliding interquartile range and the current minimum value is less than the preset lower limit of the correlation coefficient, it is determined that a related abnormality has occurred in the current walking cycle. When the above conditions are met for a number of consecutive preset walking cycles, the final relevant anomaly flag is output.

[0012] The technical effects and advantages of this invention are as follows: This invention calculates the angular velocity of the geomagnetic field vector by continuously collecting three-component geomagnetic field data and calculates the triaxial cross-correlation coefficient by capturing a complete walking cycle. When rotational anomalies and correlation anomalies occur simultaneously, a geomagnetic anomaly zone is identified. This dual-determination mechanism utilizes two independent physical characteristics: abrupt changes in the rate of change of the geomagnetic field direction and a breakdown in the phase-locked relationship of the triaxial output waveform. Both must occur simultaneously to trigger the determination, effectively eliminating false judgments caused by normal walking posture changes or single sensor spike noise. This allows investigators to obtain accurate identification in the early stages of geomagnetic interference, avoiding delays in response due to missed or false reports.

[0013] This invention forcibly blocks the azimuth output of both mechanical and electronic compasses after identifying a geomagnetic anomaly zone and issues an alarm, fundamentally cutting off the presentation of erroneous directional information. Unlike existing practices that only issue audible alerts or display warning text, the forced blocking prevents compass readings from being displayed on the navigation interface, eliminating the risk that investigators might subconsciously rely on incorrect compass guidance in emergency or fatigued states. This completely avoids repeated circling, route deviations, and the resulting safety threats such as excessive physical exertion, hypothermia, and falls caused by geomagnetic anomalies, significantly improving the safety of field geological survey operations.

[0014] This invention uses the last accurate magnetic azimuth and the last effective satellite positioning coordinates before entering a geomagnetic anomaly zone as the initial heading reference. It activates the gyroscopes and accelerometers in the inertial measurement unit to perform dead reckoning, obtaining continuous relative heading angles. Simultaneously, it activates the camera unit to calculate the heading change using visual odometry, or activates the sky polarization sensor to calculate the absolute heading angle. The heading change or absolute heading angle is used to periodically correct the dead reckoning drift, ultimately outputting a corrected navigation heading independent of the geomagnetic field. This scheme can continuously provide reliable heading even in environments where satellite signals are blocked by canyons or dense forests, and effectively suppresses time drift in inertial navigation. It allows surveyors to obtain accurate directional guidance even when completely isolated from geomagnetic field information, achieving autonomous recovery capability for field geological survey navigation in magnetic anomaly environments. Attached Figure Description

[0015] To facilitate understanding by those skilled in the art, the present invention will be further described below with reference to the accompanying drawings; Figure 1This is a schematic diagram of a navigation method applicable to field geological surveys according to the present invention. Detailed Implementation

[0016] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0017] Reference Figure 1 The following examples were obtained: Example 1: A navigation method suitable for field geological surveys, comprising the following steps: The three components of the geomagnetic field data output by the triaxial magnetometer are continuously collected. The spatial angle between the geomagnetic field vectors at two adjacent sampling times is divided by the sampling time interval to obtain the angular velocity of the geomagnetic field vector rotation. A complete walking cycle is selected as the phase analysis window. The cross-correlation coefficients between the output sequences of any two axes of the triaxial magnetometer within the phase analysis window are calculated to obtain a set of triaxial correlation coefficients containing three cross-correlation coefficients. Rotational anomalies are determined by using the rotational angular velocity of the geomagnetic field vector, and correlation anomalies are determined by using the triaxial correlation coefficient set. When both rotational and correlation anomalies are found in the judgment results, the current area is determined to be a geomagnetic anomaly zone. In response to the determination that the current area is a geomagnetic anomaly zone, the azimuth output of the mechanical compass and electronic compass is forcibly blocked, and an alarm is issued; Extract the accurate magnetic azimuth and the last valid satellite positioning coordinates recorded at the last moment before the geomagnetic anomaly zone was identified. Using the accurate magnetic azimuth as the initial heading reference, start the gyroscope and accelerometer in the inertial measurement unit to perform dead reckoning and obtain continuous relative heading angles. The camera unit is activated to acquire continuous surface images. The optical flow features between adjacent image frames are extracted using the visual odometry method to calculate the course change. Alternatively, the sky polarization sensor is activated to acquire the atmospheric polarization mode and calculate the absolute course angle. The course change or the absolute course angle is used to perform periodic drift correction on the continuous relative course angle, and a corrected navigation course that does not depend on the geomagnetic field is output.

[0018] When the navigation device determines that the current area is a geomagnetic anomaly zone, it first performs a forced shutdown operation: at the software level, it cuts off the transmission channel of azimuth data from the mechanical and electronic compasses to the display interface and navigation calculation module. Specifically, the compass driver service in the device's operating system receives a shutdown command from the anomaly detection module, which marks the compass data stream as invalid. The compass dial or digital reading area on the display interface, which was originally used to indicate direction, is replaced with a red lock icon, and the text prompt "Compass Locked" is overlaid. No azimuth data is updated. At the same time, all algorithms in the navigation calculation module that rely on magnetic azimuth are suspended, and the inertial dead reckoning module automatically takes over the direction output.

[0019] Regarding alarm issuance, the device simultaneously activates visual and auditory alarms. The visual alarm uses a full-screen flashing method: the screen background alternates between red and black, flashing twice per second for fifteen seconds, then briefly flashing once every thirty seconds until the anomaly zone detection is complete. The auditory alarm plays a high-priority voice prompt through the device's built-in speaker, stating, "Currently entering a strong magnetic anomaly zone, compass direction locked, please follow the inertial visual navigation direction." This voice prompt is played twice consecutively at the beginning of the blockade, and then repeated every sixty seconds until the compass blockade is lifted. If the device is equipped with a vibration motor, it will also simultaneously issue a tactile alarm with a two-second long vibration followed by a one-second interval, alerting investigators wearing the device on their waist or wrist. All alarm methods are designed to be non-ignorable by a single button press; the device must automatically detect that the area has been cleared of the anomaly zone and that the recovery conditions are met before stopping the alarm and lifting the blockade. The above alarm parameters can be adjusted in the device configuration interface according to the actual usage scenario, such as adjusting the flashing frequency and voice repetition interval, to adapt to the perception needs of different field environments.

[0020] Performing rotation anomaly detection and related anomaly detection refers to: During field exploration, the rotational angular velocity of the geomagnetic field vector reflects the degree of drastic change in the geomagnetic field direction with spatial location. The first historical baseline is the sliding average of all geomagnetic field vector rotational angular velocities over a continuous period before entering the current area, representing the background level of magnetic field direction change under normal conditions in that area. When the real-time measured rotational angular velocity exceeds three times the first historical baseline, it is marked as a rotational anomaly. This is because during normal field exploration, the change in geomagnetic field direction is gradual, and the rotational angular velocity is usually stable at 0.3 to 0.6 radians per second. However, when entering the boundary of a magnetic anomaly, this value can suddenly rise to over 2.0 radians per second. Three times the baseline can effectively distinguish between the slow changes caused by normal terrain undulations and the abrupt changes caused by ore body boundaries.

[0021] The second historical baseline is the moving average of all triaxial correlation coefficient sets over a continuous period prior to entering the current region. Each triaxial correlation coefficient set contains three cross-correlation coefficients, corresponding to the waveform similarity between pairs of magnetometer outputs in the north, east, and vertical directions, respectively. During normal walking, the periodic swaying of the human body results in a stable phase-locked relationship in the triaxial outputs, with cross-correlation coefficients typically between 0.85 and 0.95, and a minimum value not lower than 0.80. When the magnetic field is disturbed by a local magnetic source, the waveform correlation of the three axes is disrupted, and the minimum cross-correlation coefficient may plummet to below 0.30. A predetermined ratio of half is set, meaning that when the minimum value in the current triaxial correlation coefficient set is less than half of the corresponding minimum value in the second historical baseline, it is marked as a correlation anomaly. This is because a 50% drop far exceeds the fluctuation range caused by minor gait changes during normal walking, avoiding misjudgment while ensuring that no true magnetic anomalies are missed. For example: Suppose that within the first 30 seconds of entering a basalt outcrop area, the sliding average value of the rotational angular velocity of the geomagnetic field vector is 0.5 radians per second, and the real-time measured rotational angular velocity is 1.6 radians per second, which exceeds three times the baseline (1.5), triggering a rotational anomaly.

[0022] Meanwhile, the moving average of the minimum value of the triaxial correlation coefficient set during this period is 0.90, while the minimum value measured during the current walking cycle is 0.35, which is less than half (0.45), triggering a correlation anomaly. With both anomalies occurring simultaneously, the system determines that it has entered a geomagnetic anomaly zone. If only one condition is met, such as a high rotational angular velocity but a normal correlation coefficient, it may be due to a brief phenomenon of rapid turning by personnel or terrain undulations, and will not trigger the final judgment, thus avoiding a false alarm.

[0023] The method for extracting the complete walking cycle is as follows: The vertical acceleration value output by the accelerometer in the inertial measurement unit is continuously monitored. This vertical acceleration value changes regularly with the rise and fall of the feet during human walking. When the human body pushes off the ground, the vertical acceleration is positive, and when the foot lands and decelerates, the vertical acceleration is negative. Therefore, a complete step is always accompanied by the vertical acceleration value passing through zero from positive to negative or passing through zero from negative to positive.

[0024] The criterion of crossing the zero point twice consecutively is used because a single zero-point crossing might be a false zero-point crossing caused by noise interference or slight body swaying, while two consecutive zero-point crossings truly correspond to a complete movement of lifting a leg and landing. The time interval between the two crossings needs to be within a preset normal cadence range. This preset normal cadence range is determined based on statistical cadence data from field geological surveyors walking under load, typically ranging from 0.8 seconds to 1.5 seconds, corresponding to 40 to 75 steps per minute. When two consecutive zero-point crossings are detected simultaneously and the time interval is within this range, the time interval between the two zero-point crossings is considered a complete walking cycle. For example, if the vertical acceleration crosses the zero point from top to bottom at time 1.2 seconds and then from bottom to top at time 2.3 seconds, the time interval is 1.1 seconds, falling within the preset range. Then, the data between time 1.2 seconds and 2.3 seconds constitutes a complete walking cycle. In cases where gait frequency cannot be detected, such as when a person is standing still or walking at a very slow speed, resulting in an acceleration amplitude that is too small to generate a stable zero crossing point, a fixed duration is used as a complete walking cycle. This fixed duration is set based on the average gait frequency cycle of the human body, usually 1.2 seconds, to ensure that a phase analysis window can still be obtained when there is no gait information.

[0025] The cross-correlation coefficient is calculated as follows: For any two axis output sequences, such as the northward component sequence and the eastward component sequence, multiply the corresponding data points of the two sequences within the same phase analysis window and sum them. This sum reflects the synchronicity of the two sets of data at the same time position. Then divide this sum by the product of the square roots of their respective squares, where the square roots represent the energy magnitude of each sequence. After dividing, the normalized correlation coefficient is obtained. The cross-correlation coefficient calculated in this way ranges from -1 to +1, where +1 indicates perfect positive correlation, -1 indicates perfect negative correlation, and zero indicates no correlation. During normal walking, the three-axis output waveforms within the same phase analysis window exhibit a high degree of consistency due to the fixed pattern of human body swaying, and the cross-correlation coefficient is usually above 0.85. For example, if the northward component sequence is [0.5, 0.8, 0.3, -0.1, -0.4] and the eastward component sequence is [0.4, 0.7, 0.2, -0.2, -0.5] during a complete walking cycle, the cross-correlation coefficient calculated according to the above formula is approximately 0.92, indicating that the two axes have a good synchronization relationship.

[0026] The specific implementation method of dead reckoning is as follows: The initial heading angle is the accurate magnetic azimuth recorded at the last moment before entering the geomagnetic anomaly zone, and the initial position coordinates are the last valid satellite positioning coordinates. This is because mechanical and electronic compasses are forcibly blocked after entering the geomagnetic anomaly zone and can no longer provide reliable directional information. Therefore, it is necessary to rely on the last reliable magnetic azimuth and satellite positioning coordinates before entering the anomaly zone as the starting point for subsequent dead reckoning. For example, suppose that at the moment of determining entry into the geomagnetic anomaly zone, the last recorded accurate magnetic azimuth is 30 degrees east of north, and the last valid satellite positioning coordinates are 118.2 degrees east longitude and 39.8 degrees north latitude. Then, the navigation equipment will use 30 degrees as the initial heading angle for dead reckoning and these latitude and longitude coordinates as the initial position coordinates.

[0027] By integrating the three-axis angular velocities output by the gyroscope over time, the change in heading angle relative to the initial heading angle at each instant is obtained. The gyroscope measures the angular velocity of the carrier around its three axes, in radians per second. Continuously summing these angular velocity values ​​over time yields the total rotation angle from the initial moment to the current moment. For example, if the gyroscope's heading angular velocity is ten degrees per second within 0.1 seconds, then the heading angle change after 0.1 seconds is one degree. The longer the integration time, the larger the accumulated value. To prevent excessive accumulation of integration drift, it is necessary to periodically correct using heading changes from other sensors.

[0028] Gait detection utilizes the triaxial acceleration output from an accelerometer. A complete step is defined as two consecutive zero-point crossings of the vertical acceleration value. The vertical acceleration value exhibits a regular alternation of positive and negative values ​​during walking: positive when the foot leaves the ground and negative when the foot lands and decelerates. A complete step necessarily involves two zero-point crossings (once from positive to negative, and once from negative to positive). Two consecutive zero-point crossings ensure that a complete closed-loop movement of leg lift and landing is detected, rather than a single zero-point crossing caused by body swaying. For example, if the vertical acceleration crosses the zero point from top to bottom in 1.2 seconds and from bottom to top in 2.3 seconds, with an interval of 1.1 seconds, then a step is determined to have been completed between 1.2 and 2.3 seconds.

[0029] Step length is estimated based on the composite amplitude of acceleration, which is proportional to the fourth root of the acceleration amplitude. The composite amplitude of acceleration is equal to the square root of the sum of the squares of the three-axis accelerations, reflecting the force exerted by the human body when pushing off the ground during walking. Experimental data shows that when walking with a load on rugged terrain, the composite amplitude of acceleration is approximately linearly related to the fourth root of the amplitude when it ranges from 1.2 to 1.8 times the gravitational acceleration. For example, with a proportionality coefficient of 0.45, when the composite amplitude of acceleration is 1.5 times the gravitational acceleration, the fourth root is approximately 1.1, and the estimated step length is approximately 0.5 meters. When the composite amplitude of acceleration increases to 1.8 times the gravitational acceleration, the fourth root is approximately 1.16, and the estimated step length is approximately 0.52 meters. This pattern conforms to the biomechanical characteristics of human walking, and step length estimation can be achieved without additional sensors.

[0030] The step size is multiplied by the cosine and sine of the current heading angle and then added to the position coordinates to obtain a continuous position calculation result. The current heading angle is equal to the initial heading angle plus the heading angle change obtained by gyroscope integration, and then further corrected. For example, if the current heading angle is 30 degrees and the estimated step size is 0.5 meters, then the position increment in the north direction is 0.5 multiplied by the cosine of 30 degrees, approximately 0.43 meters, and the position increment in the east direction is 0.5 multiplied by the sine of 30 degrees, approximately 0.25 meters. These two increments are added to the previously calculated position coordinates to obtain the new position at the current moment. This accumulation is performed after each step, thus continuously updating the position.

[0031] After each predetermined number of steps, the angular velocity drift of the integral is zeroed out using the heading change calculated by visual odometry or a sky polarized light sensor. The correction is the average difference between the gyroscope integral heading and the visual heading within the most recent predetermined number of steps. The predetermined number of steps is typically ten, as ten steps correspond to a walking distance of approximately five to eight meters, within which the visual odometry and polarized light sensor can provide reliable relative heading changes. For example, if the total heading change obtained by the gyroscope integration within the most recent ten steps is fifteen degrees, while the total heading change calculated by visual odometry is twelve degrees, the difference is three degrees, so the correction is 0.3 degrees per step. This correction is evenly distributed across the heading angle of each subsequent step to eliminate the gyroscope's zero-point drift error. The corrected heading angle is then used for the next position calculation, ensuring that position accuracy does not deteriorate due to integral drift during long-term walking.

[0032] The specific steps for extracting optical flow features between adjacent image frames using the visual odometry method are as follows: The continuous surface images captured by the camera unit are converted into grayscale images. Color images contain a large amount of data across three channels (red, green, and blue). After conversion to grayscale, each pixel has only one brightness value, reducing the computational load to one-third of the original. At the same time, the light and dark variations in the surface texture are completely preserved. For example, the rock texture on the ground is composed of different grayscale gradients, which remain clearly discernible in the grayscale image.

[0033] For each grayscale image, fast corner features are extracted. Fast corner feature extraction is a computationally efficient corner detection method that determines whether a point is a corner by comparing the grayscale values ​​of sixteen pixels around the candidate pixel's circumference. A grayscale difference threshold of twenty is set; when the grayscale values ​​of nine consecutive pixels on the circumference differ from the grayscale value of the center pixel by more than twenty, the center pixel is marked as a corner. Corners typically correspond to features with obvious directional changes in surface images, such as the edges of gravel, the tips of vegetation leaves, and soil particles. These features are easily and stably tracked between adjacent frames. For example, in an image taken on a gravel road, fast corner feature detection can extract two to three hundred evenly distributed corners.

[0034] The pyramid optical flow method is used to track the pixel displacement of corner features between two adjacent image frames. The pyramid optical flow method is a sparse optical flow tracking algorithm in existing technology. Its core operation is to construct an image pyramid by downsampling the original image layer by layer. For example, the pyramid has three layers: the first layer is the original image, the second layer reduces the original image's width and height by half, and the third layer reduces it by half again. A coarse optical flow displacement is calculated in the top layer, the image with the lowest resolution, and then the result is passed to the next layer for refinement, calculating the final pixel displacement layer by layer down to the original image layer. This process can track larger pixel displacements. For example, when an investigator walks at a speed of 0.5 meters per second, the pixel movement of the same ground feature between two adjacent image frames may reach 20 to 30 pixels. Single-layer optical flow methods are prone to losing tracking, while the three-layer pyramid structure can stably track within this range. For each fast corner feature extracted from the previous frame, the pixel point whose grayscale distribution best matches it is searched in each layer of the pyramid in the current frame. The lateral and vertical movement of this corner point on the image plane are calculated, i.e., the pixel displacement.

[0035] Median filtering is applied to the pixel displacements of all successfully tracked corner points to remove outliers whose displacement vector magnitudes exceed a predetermined multiple of the median standard deviation. Median filtering involves arranging the magnitudes of all successfully tracked corner point displacement vectors in ascending order, taking the median value as the median, and then calculating the median standard deviation as the median of the absolute deviations of all magnitudes from the median. The predetermined multiple is set to three times the median standard deviation. When the magnitude of a displacement vector is greater than the median plus three times the median standard deviation or less than the median minus three times the median standard deviation, it is identified as an outlier and removed. For example, assuming that among one hundred successfully tracked corner points, the median displacement vector magnitude is fifteen pixels and the median standard deviation is three pixels, then corner points with a magnitude less than six pixels or greater than twenty-four pixels are removed. These outliers are usually caused by swaying leaves, changes in lighting, or slight camera unit vibration. The remaining corner points after removal reflect the true movement trend of the ground surface in the horizontal plane.

[0036] The remaining effective pixel displacements are converted into unit pixel movement vectors. Then, using the camera unit's intrinsic parameter matrix and the assumed height of the carrier, the carrier's heading change in the horizontal plane is deduced. The unit pixel movement vector is a direction vector of length 1 obtained by removing the pixel value of each effective corner point from its magnitude, representing the direction of movement of that corner point. The unit pixel movement vectors of all effective corner points are averaged to obtain the average direction of the overall movement. The camera unit's intrinsic parameter matrix includes the focal length and principal point coordinates. These parameters are calibrated at the factory and are used to transform pixel coordinates on the image plane to spatial orientation in the camera coordinate system. The assumed height of the carrier refers to the estimated height of the camera unit above the ground. In geological surveys, the camera unit is usually fixed at chest level or waist level, with a height of approximately 0.8 meters to 1.2 meters above the ground. The midpoint value of 1.0 meters is used for calculation.

[0037] Based on the pinhole camera model, given the pixel displacement, intrinsic parameter matrix, and height, the rotation angle of the carrier in the horizontal plane relative to the last image acquisition time can be deduced through geometric relationships; that is, the change in heading. For example, if the average direction of the overall pixel displacement is five degrees to the right, the calculated change in heading might correspond to a 3.5-degree rightward rotation of the carrier.

[0038] When the number of valid corner points falls below a predetermined number, the visual odometry method is deemed ineffective, and the system automatically switches to the sky polarization sensor for heading calculation. The predetermined number can be set to ten, because with fewer than ten valid corner points, the vector averaging result is easily affected by individual residual outliers and becomes unreliable. For example, when walking on smooth rock surfaces, snow, or uniform sand, the grayscale image lacks sufficient texture features, and fast corner detection can only extract five or six corner points. Furthermore, some of these corner points may fail to be tracked between adjacent frames, leaving fewer than ten valid corner points. In this case, the visual odometry method is deemed ineffective, and the navigation device automatically switches to the sky polarization sensor to obtain the absolute heading angle.

[0039] The specific steps for the sky polarization sensor to acquire the atmospheric polarization pattern and calculate the absolute heading angle are as follows: When using a time-division multiplexing method, the sky polarization sensor sequentially blocks the light paths in two directions, collecting atmospheric polarization degree and polarization azimuth angles in the zenith direction, the direction to the left of the sun at a predetermined angle, and the direction to the right of the sun at a predetermined angle. When using a pixel-based method, the sensor simultaneously receives light from three directions through different pixel areas, synchronously acquiring data from all three directions. The three directions are selected from the zenith direction and the directions to the left and right of the sun at predetermined angles. This predetermined angle is 45 degrees because the atmospheric polarization information is richest at 45 degrees and less affected by direct sunlight. For example, if the azimuth angle of the sun in the horizontal coordinate system is 180 degrees (due south), then the two lateral directions are taken as 135 degrees and 225 degrees respectively.

[0040] The sky-polarized light sensor uses a time-division multiplexing method to sequentially measure the degree of atmospheric polarization and polarization azimuth angle in the zenith direction, 45 degrees to the left of the sun, and 45 degrees to the right of the sun. The Rayleigh scattering model states that a single scattering of sunlight by atmospheric molecules produces a specific polarization pattern, in which the polarization azimuth angle of the scattered light satisfies a definite geometric relationship with the sun's direction, the observation direction, and the scattering angle. Specifically, for any observation direction in the sky, its polarization azimuth angle is equal to the tangent direction of the great circle containing that direction and the sun's direction. This relationship can be expressed by the formula that the tangent of the polarization azimuth angle equals the sine of the difference between the sun's azimuth angle and the observation azimuth angle multiplied by a certain function. However, in actual calculations, the system of equations is usually established directly using the measurements from the three directions.

[0041] Given the zenith angle (the angle between each direction and the zenith) and azimuth angle (relative to geographic north) of three observation directions, the polarization azimuth angles measured in the three directions are P1, P2, and P3, respectively. According to the Rayleigh scattering model, there is a nonlinear functional relationship between the polarization azimuth angle Pi of each direction and the solar azimuth angle S and altitude angle H: tan(Pi) = sin(S-Ai) / (cos(H)*tan(θi)-sin(H)*cos(S-Ai)), where Ai is the azimuth angle of the i-th observation direction, and θi is the zenith angle of that direction. Substituting the measured values ​​of the three observation directions, we obtain a system of three equations containing two unknowns (solar azimuth angle S and solar altitude angle H). By solving the system iteratively using the least squares method or by using analytical methods (e.g., first eliminating H through the two equations, and then solving for S), we can deduce the solar azimuth angle S in the horizontal coordinate system.

[0042] For example: Suppose the device's current location is 118 degrees east longitude and 40 degrees north latitude, and the device's internal clock shows local time as 14:30 on October 15, 2024 (Beijing time). At this time, the theoretical azimuth of the sun is calculated to be 198 degrees (southwest). The sky polarization sensor measures the zenith direction (zenith angle 0 degrees, azimuth angle arbitrary), the direction 45 degrees to the left of the sun (taken relative to the current device orientation, but to avoid cyclic dependence, the sensor first roughly determines the sun's direction area during actual measurement), and the direction 45 degrees to the right of the sun. The measured polarization azimuth angles in the three directions are as follows: the polarization azimuth angle in the zenith direction is 0 degrees (because the polarization direction in the zenith direction is always perpendicular to the great circle of the sun direction. If the sun azimuth angle is 198 degrees, then the polarization direction in the zenith direction should be 198-90=108 degrees. In this example, we assume that the measured value is exactly 108 degrees, but in order to illustrate the calculation process, we set a set of typical measured values). We assume that the sun altitude angle is 30 degrees and the sun azimuth angle is 180 degrees (due south).

[0043] The three observation directions are set as follows: Direction 1 (zenith) with a zenith angle of 0 degrees and an arbitrary azimuth angle, its theoretical polarization azimuth angle is 90 degrees (because the polarization direction of the zenith is perpendicular to the direction of the sun); Direction 2 (45 degrees to the left of the sun) with a zenith angle of 45 degrees and an azimuth angle of 135 degrees, the measured polarization azimuth angle is -45 degrees; Direction 3 (45 degrees to the right of the sun) with a zenith angle of 45 degrees and an azimuth angle of 225 degrees, the measured polarization azimuth angle is 45 degrees. Substituting these three sets of data (zenith angle, azimuth angle, polarization azimuth angle) into the Rayleigh scattering equation, the solar azimuth angle is calculated to be 180 degrees and the altitude angle to be 30 degrees. Then, the calculated solar azimuth angle of 180 degrees is compared with the theoretical solar azimuth angle calculated by the internal clock of the equipment based on time and geographical location (assuming the theoretical value is 178 degrees at this time). The difference of 2 degrees is the deviation of the current orientation of the equipment from geographical north (if the equipment is facing north, the two values ​​should be equal). Through this calculation, the absolute heading angle can be obtained.

[0044] When the sky is obscured by trees, mountains, or clouds, preventing the sensor from simultaneously acquiring stable polarization data in all three valid directions, the device automatically saves the most recently successfully calculated absolute heading angle and uses it as the baseline for subsequent dead reckoning. The saved absolute heading angle is used continuously during the obscuration period, but the device checks sky conditions every second. Once valid measurements in all three directions are restored, the absolute heading angle is immediately recalculated and updated to eliminate accumulated errors caused by solar motion. For example, after traveling for 30 minutes in dense forest and exiting the obscured area, the sun has moved approximately 7.5 degrees; the device will then recalculate and obtain a new absolute heading angle, updating the baseline.

[0045] By comparing the deduced solar azimuth with the theoretical solar azimuth provided by the device's internal clock (which contains precise year, month, day, hour, minute, and second information), and combining this with pre-stored latitude and longitude coordinates of the current location, the theoretical solar azimuth in the horizontal coordinate system at that moment can be calculated. The difference between the deduced and theoretical azimuth is the angle of deviation of the device's current orientation relative to geographic north. For example, if the device deduces a solar azimuth of 190 degrees, while the internal clock calculates a theoretical solar azimuth of 180 degrees, the difference is +10 degrees, indicating that the device's current orientation is 10 degrees east of geographic north. This difference is the absolute heading angle.

[0046] When obstruction prevents simultaneous measurement of three valid directions—such as dense tree canopy, mountainsides on either side of a canyon, or clouds obscuring the zenith or a lateral direction—the system automatically saves the most recently successfully calculated absolute heading angle and uses it as the baseline for subsequent dead reckoning until a valid measurement is obtained again. The saved absolute heading angle remains valid for no more than thirty minutes, as the sun's position changes over time, and prolonged periods without updates can lead to error accumulation. For example, when navigating through dense forest, if the sensors cannot simultaneously collect valid data in three directions, the system locks the absolute heading angle of 30 degrees calculated the last time the system exited the forest. Afterward, it maintains the heading solely through gyroscope integration and visual odometry until it exits the obstructed area, at which point the absolute heading angle is remeasured and updated.

[0047] The real-time calculation of the rotational angular velocity of the geomagnetic field vector and the identification of rotational anomalies are carried out using the following steps: A three-component geomagnetic field data sequence was established. Each sampling point included a northward component, an eastward component, and a vertical component. The three components of the geomagnetic field vector changed continuously as the device moved in space, forming a time series. The sampling interval was between 0.02 seconds and 0.1 seconds. This frequency can capture the rapid changes in the direction of the geomagnetic field without generating too much data.

[0048] For two geomagnetic field vectors at two adjacent sampling times, their spatial angle is calculated. The spatial angle is equal to the dot product of the two vectors divided by the product of their magnitudes, and then the inverse cosine is taken. This spatial angle reflects the actual angle that the geomagnetic field direction has rotated in three-dimensional space from the previous moment to the current moment, and is not affected by the magnitude of the magnetic field strength. For example, when walking in a region with a uniform geomagnetic field, the magnetic field direction is almost the same at different times, and the angle is close to zero degrees; when entering the boundary of a magnetite body, the magnetic field direction may deflect by twenty degrees in 0.1 seconds, and the spatial angle increases significantly.

[0049] Dividing the spatial angle by the sampling time interval yields the instantaneous rotational angular velocity, measured in radians per second, which indicates how quickly the direction of the geomagnetic field changes over time. During normal outdoor walking, due to changes in ground slope and equipment vibration, the instantaneous rotational angular velocity is typically in the range of 0.3 to 0.6 radians per second.

[0050] The average instantaneous rotational angular velocity is calculated within a sliding window. The sliding window length is a preset value, ranging from 0.6 seconds to 3 seconds, corresponding to 30 to 150 sampling points. This window length can smooth out random noise from individual sampling points while preserving the key characteristics of abrupt changes in magnetic field anomalies. For example, when the sliding window has 30 sampling points (0.6 seconds), the average value represents the general level of the rotational angular velocity within the most recent 0.6 seconds.

[0051] The coefficient of variation (COP) of the rotational angular velocity within the sliding window is calculated. The COP is equal to the standard deviation of the angular velocity within the window divided by the mean. The COP is a dimensionless statistic that reflects the degree of fluctuation in the rotational angular velocity within the window. When the magnetic field is stable, the angular velocity fluctuation is small, and the COP is usually below 0.2. When the magnetic field is strongly disturbed, the angular velocity fluctuates greatly, and the COP may exceed 0.5.

[0052] The average value of the coefficient of variation over a continuous period prior to entering the current region is calculated as the historical baseline coefficient of variation. This continuous period consists of data from 30 to 60 seconds prior to entering the current region. The historical baseline coefficient of variation represents the stability of the rotational angular velocity of the region under normal conditions. For example, when traversing a typical sedimentary rock area, the historical baseline coefficient of variation is approximately 0.15; while in a basalt area, even without significant magnetic anomalies, the historical baseline coefficient of variation may reach 0.25.

[0053] When the coefficient of variation of the historical baseline is less than a preset variation threshold, a rotational anomaly is identified when the angular velocity of the geomagnetic field vector exceeds a predetermined multiple of the first historical baseline. The preset variation threshold is set to 0.3, because when the coefficient of variation of the historical baseline is below 0.3, it indicates that the background fluctuation of the geomagnetic field before entering the region is relatively small, and an anomaly can be reliably detected simply by the angular velocity of rotation exceeding three times the first historical baseline. For example, if the coefficient of variation of the historical baseline is 0.2, the first historical baseline is 0.5 radians per second, and the real-time angular velocity increases to 1.6 radians per second, exceeding three times the baseline, a rotational anomaly can be identified.

[0054] When the historical baseline coefficient of variation is greater than or equal to a preset variation threshold, a rotational anomaly is identified when both the geomagnetic field vector rotational angular velocity exceed a predetermined multiple of the first historical baseline and the current coefficient of variation exceeds a predetermined multiple of the historical baseline coefficient of variation. This condition is added because when the geomagnetic field background before entering the region already exhibits significant fluctuations (coefficient of variation higher than 0.3), relying solely on a rotational angular velocity exceeding three times the baseline may lead to frequent false alarms due to background fluctuations. For example, if the historical baseline coefficient of variation is 0.35, the first historical baseline is 0.6 radians per second, the real-time rotational angular velocity is 1.9 radians per second (exceeding three times the baseline), but the current coefficient of variation is only 0.2 (not exceeding a predetermined multiple of the historical baseline coefficient of variation, e.g., 1.5 times, or 0.525), this indicates that the increase in angular velocity is merely the normal upper limit of background fluctuations, not a true magnetic anomaly jump, and therefore, a rotational anomaly is not identified. Only when the current coefficient of variation also simultaneously exceeds 1.5 times the historical baseline coefficient of variation is it confirmed that the magnetic field direction change is indeed in a state of drastic and abnormal fluctuation, thereby reducing the false alarm rate in complex geological backgrounds. The predetermined multiple is 1.5 times the current coefficient of variation. This value has been tested in the field and ensures that the false alarm rate is controlled to below 5% while maintaining the detection rate.

[0055] When a predetermined number of consecutive sampling moments are marked as rotational anomalies, the final rotational anomaly marker is output. The predetermined number of consecutive sampling moments is five, corresponding to a time window of 0.1 to 0.5 seconds. A single marker may be sporadic noise; if the condition is met five times consecutively, it proves that the magnetic field is indeed in a continuous abnormal state, avoiding erroneous outputs due to transient interference. For example, if five consecutive sampling points within 0.2 seconds meet the above branch judgment condition, the system will output the final rotational anomaly marker, which will then be used in conjunction with related anomalies to determine the geomagnetic anomaly region.

[0056] The following steps are used to determine the correlation anomalies in the triaxial correlation coefficient set: The first step involves denoting the northward, eastward, and vertical component sequences of the geomagnetic field data collected within a complete walking cycle as the first, second, and third sequences, respectively. The length of each sequence is equal to the number of sampling points within that walking cycle. The duration of a complete walking cycle is determined by the step frequency, typically between 0.8 and 1.5 seconds. Therefore, the number of sampling points in each sequence equals the sampling frequency multiplied by the cycle duration. For example, if the sampling frequency is 50 Hz and a complete walking cycle is 1.2 seconds, then each sequence contains 60 sampling points. Naming the three sequences separately facilitates differentiation when calculating cross-correlation coefficients later.

[0057] The second step is to calculate the cross-correlation coefficient between the first and second sequences. The cross-correlation coefficient is equal to the product of the covariance of the two sequences divided by their respective standard deviations. Similarly, calculate the cross-correlation coefficients between the first and third sequences, and between the second and third sequences, to obtain the set of three-axis correlation coefficients. The cross-correlation coefficient is a value between -1 and +1, where +1 indicates a perfect positive correlation, -1 indicates a perfect negative correlation, and zero indicates no correlation. The covariance is calculated by subtracting the mean of each sequence from the data at each corresponding position, multiplying the deviations at corresponding positions, summing the results, and finally dividing by the sequence length minus one. The standard deviation is the square root of the sum of the squared deviations of each sequence from its mean, divided by the length minus one.

[0058] For example, suppose the first sequence (northward component) is [0.5, 0.8, 0.3, -0.1, -0.4] and the second sequence (eastward component) is [0.4, 0.7, 0.2, -0.2, -0.5]. The calculated covariance is approximately 0.86, and the two standard deviations are approximately 0.47 and 0.46, respectively. The cross-correlation coefficient is approximately 0.86 divided by (0.47 multiplied by 0.46), which equals 0.99, indicating that the two sequences are highly correlated. After calculating the cross-correlation coefficients of the three pairs of combinations, we obtain a set of three values, such as {0.99, 0.97, 0.96}.

[0059] The third step involves calculating the statistical characteristics of the correlation coefficients for a predetermined number of complete walking cycles before entering the current area. The minimum value in the triaxial correlation coefficient set within each cycle is recorded, forming a historical minimum value sequence. The predetermined number of cycles is twenty complete walking cycles, corresponding to a time window of approximately twenty to thirty seconds. Each cycle will inevitably contain a minimum value in the triaxial correlation coefficient set; for example, the minimum value in the first cycle might be 0.94, and in the second cycle, it might be 0.92. Collecting twenty such values ​​constitutes the historical minimum value sequence. [0.94,0.92,0.95,0.91,0.93,0.90,0.92,0.94,0.91,0.89,0.93,0.92,0.90,0.91,0.88,0.92,0.93,0.90,0.91,0.89], this sequence reflects the normal fluctuation range of the minimum value of the triaxial correlation coefficient of the magnetic field before entering the current region.

[0060] The fourth step is to calculate the moving median and moving interquartile range of the historical minimum value sequence. The moving median is the middle value after arranging the values ​​in the sequence in ascending order. For twenty values, the median is the average of the tenth and eleventh values. The interquartile range is the difference between the upper quartile (75th percentile) and the lower quartile (25th percentile). For example, after sorting the above sequence, we get: [0.88,0.89,0.89,0.90,0.90,0.90,0.91,0.91,0.91,0.92,0.92,0.92,0.93,0.93,0.94,0.94,0.95,0.95], the lower quartile is 0.91 (the tenth quartile), the upper quartile is 0.94 (the fifteenth quartile), the interquartile range is 0.03, and the median is 0.915. This median and interquartile range are used to generate a dynamic anomaly detection threshold.

[0061] The fifth step is to calculate the minimum value in the set of triaxial correlation coefficients for the current walking cycle. Assuming the triaxial cross-correlation coefficients measured in the current cycle are 0.45, 0.48, and 0.42, the minimum value is 0.42. This minimum value represents whether the phase lock between the three axes of the current magnetic field has broken down.

[0062] Step 6: When the current minimum value is less than the moving median minus 1.5 times the moving interquartile range and less than the preset lower limit of the correlation coefficient, a correlation anomaly is determined to have occurred in the current walking cycle. The calculation result of the moving median minus 1.5 times the moving interquartile range is 0.915 minus 1.5 multiplied by 0.03 equals 0.87. The current minimum value of 0.42 is obviously less than 0.87. At the same time, the preset lower limit of the correlation coefficient is 0.30, and 0.42 is less than 0.30, so the condition is not met. Therefore, no correlation anomaly is triggered in this example. If the current minimum value drops to 0.25, then both conditions of being less than 0.87 and less than 0.30 are met, and a correlation anomaly is determined. 1.5 times the moving interquartile range is a commonly used parameter in statistical outlier detection. It can effectively identify a sudden drop in correlation caused by magnetic anomalies while ensuring that more than 99% of normal samples fall within the threshold.

[0063] The preset lower limit of the correlation coefficient of 0.30 is determined based on a large amount of field measurement data. When the cross-correlation coefficients of all three axes are below 0.30, it indicates that the outputs of the three axes are completely out of sync, which must be caused by strong magnetic interference, forming a double insurance with the statistical threshold. When the above conditions are met for a preset number of consecutive walking cycles, the final correlation anomaly flag is output. The preset number of consecutive cycles is three walking cycles, because a single correlation anomaly may be an accidental gait anomaly or a brief disturbance. Only when three consecutive cycles are abnormal can it be confirmed that the magnetic field has entered a state of continuous instability, eliminating false alarms. For example, if the minimum values ​​of the three consecutive cycles are 0.25, 0.22, and 0.28, respectively, all less than 0.87 and 0.30, then a correlation anomaly flag is output.

[0064] All parameters in this method, such as predetermined multiples, predetermined ratios, predetermined times, predetermined thresholds, predetermined quantities, and predetermined lower limits of correlation coefficients, are not fixed values. Instead, they can be adjusted through experimental calibration before equipment calibration or field use, based on the actual application scenario, sensor model, geological background noise level, and required detection sensitivity and false alarm rate. The specific selection criteria and reasonable value ranges for these parameters are as follows: The rotational angular velocity exceeding the predetermined multiple of the first historical baseline is selected based on the statistical distribution of the random fluctuation amplitude of rotational angular velocity during normal walking and the abrupt change amplitude of magnetic anomalies. Those skilled in the art calculate the standard deviation by collecting rotational angular velocity data under conditions without magnetic anomalies in the target area, typically choosing a multiple between 2 and 5 times. Too low a multiple may lead to false alarms due to normal walking swaying, while too high a multiple may miss weak magnetic anomalies. A multiple of 3 times is recommended based on the empirical rule that 99.7% of the samples in a Gaussian distribution fall within the range of plus or minus three standard deviations of the mean, but in a noisy environment, it can be adjusted to 4 or 5 times.

[0065] The predetermined proportion of the minimum value in the triaxial correlation coefficient set that is lower than the minimum value corresponding to the second historical baseline is selected based on the diurnal fluctuation range of the correlation coefficient during normal walking. The confidence interval is calculated by collecting the minimum correlation coefficient values ​​from at least 200 normal walking cycles in the target area beforehand, and is typically set to 50% to 70% as the alarm threshold. A predetermined proportion of one-half (i.e., lower than half of the historical baseline) is the most commonly used low-sensitivity threshold, but it can be adjusted to two-thirds or three-fifths in areas with complex geological backgrounds.

[0066] The upper and lower limits of the preset normal cadence range are determined based on physiological data of human walking. The lower limit is usually 0.6 to 0.8 seconds, corresponding to the upper limit of cadence for brisk walking, and the upper limit is usually 1.5 to 1.8 seconds, corresponding to the lower limit of cadence for slow, weight-bearing walking. Specific values ​​can be fine-tuned on-site based on the average height and load weight of the participants. The fixed duration is based on the median of the average human cadence cycle and is usually set to any value between 1.0 and 1.5 seconds. When cadence cannot be detected, this fixed duration is used as a temporary walking cycle, with an allowable error of up to 20%.

[0067] The preset sliding window length, used to calculate the average rotational angular velocity and coefficient of variation, is determined based on the sensor sampling frequency and the desired time resolution. At a sampling frequency of 50 Hz, a window length between 30 and 150 sampling points corresponds to 0.6 to 3.0 seconds. A window that is too short generates significant noise, while a window that is too long results in lag; a recommended value is 50 sampling points (1 second). The predetermined number of steps (the interval between steps to correct for gyroscope integral drift) is determined based on the minimum distance from which a visual odometry or polarized light sensor provides a reliable change in heading. Typically, 6 to 15 steps are used, equivalent to a walking distance of 3 to 10 meters. Too few steps result in frequent corrections but may introduce visual measurement noise, while too many steps lead to excessive drift accumulation; 10 steps are recommended.

[0068] The median standard deviation (MAD) multiple is used to remove anomalous pixel displacement vectors in optical flow tracking. It is determined based on the difference in distribution between normal and anomalous displacements in the field image. Typically, it is 2.5 to 4 times the MAD. A value below 2.5 times may mistakenly delete normal corner points, while a value above 4 times may miss anomalous points; 3 times is recommended. The predetermined number (the lower limit of the effective corner point count when visual odometry fails) is determined based on the number of samples required for vector averaging in optical flow tracking. Typically, it is 8 to 15. Fewer than 8 samples result in unstable statistics, while more than 15 samples lead to high computational cost and diminishing returns; 10 samples are recommended.

[0069] The predetermined angle (the angles on either side of the sun in the three measurement directions of sky-polarized light) is determined based on the analysis of the sensitivity of polarization information using the Rayleigh scattering model. It is typically set between 30 and 60 degrees. If the angle is too small, the two lateral directions are too close to the sun, resulting in low polarization; if the angle is too large, the scattering angle increases and the signal weakens. 45 degrees is recommended. The preset variation threshold (the threshold for distinguishing between stable and unstable magnetic field backgrounds) is determined based on the statistical distribution of the coefficient of variation from a large amount of field magnetic field data. The coefficient of variation is typically less than 0.2 in areas without magnetic anomalies, between 0.2 and 0.4 in areas with slight magnetic backgrounds, and greater than 0.5 in areas with strong magnetic anomalies. A preset variation threshold of 0.3 can effectively distinguish between stable and fluctuating backgrounds.

[0070] The current coefficient of variation exceeds the historical baseline coefficient of variation by a predetermined multiple to increase the stringency of the criterion in areas of background fluctuation. Typically, this is 1.2 to 2.0 times; exceeding 2.0 times is too stringent and may result in missed detections, while below 1.2 times is equivalent to no additional conditions; 1.5 times is recommended. The predetermined number of consecutive sampling times (the number of consecutive anomalies required to output the final rotating anomaly marker) is determined based on the sampling frequency and expected response time. At a sampling frequency of 50 Hz, 3 to 10 consecutive anomalies corresponding to 0.06 to 0.2 seconds can effectively filter out isolated noise spikes; 5 times is recommended. The preset number of periods (the number of walks required to establish the historical minimum value sequence) is determined based on statistical stability requirements. Typically, 15 to 30 periods are used; too few periods result in large statistical fluctuations, while too many periods lead to a sluggish response to background changes; 20 periods are recommended.

[0071] The 1.5-fold sliding interquartile range (box plot outlier criterion) is a statistically accepted standard for detecting mild outliers, typically fixed at 1.5 times, but can be increased to 2 or 3 times to reduce sensitivity depending on actual needs. The preset lower limit of the correlation coefficient (absolute low correlation threshold) is determined based on the measured range of cross-correlation coefficients when the triaxial magnetometer completely loses phase lock. Field measurements show that the cross-correlation coefficient rarely falls below 0.4 during normal walking; when it falls below 0.3, it is almost certainly caused by magnetic interference. Therefore, the preset lower limit of the correlation coefficient can be any value between 0.2 and 0.4, with 0.3 recommended. The preset number of consecutive cycles (the number of consecutive abnormality cycles required to output the final correlation anomaly marker) is determined based on the human walking cycle and the spatial continuity of the magnetic anomaly. Typically, 2 to 5 consecutive cycles are used; fewer than 2 is insufficient to confirm persistence, and more than 5 may miss short-term magnetic anomalies; 3 cycles are recommended. All parameters in this method are adjustable design parameters. The specific values ​​are conventional optimizations made by those skilled in the art based on actual working conditions without changing the core discrimination logic of this invention.

[0072] The above algorithms or formulas are all dimensionless and numerical calculations, and the results are obtained by software simulation based on a large amount of collected data to obtain the most recent real-world results. The preset parameters are set by those skilled in the art according to the actual situation.

[0073] It should be understood that in the various embodiments of this application, the order of the above-mentioned processes does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0074] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0075] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the devices and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0076] The above are merely specific embodiments 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 scope of the technology 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 navigation method suitable for field geological survey, characterized by, Includes the following steps: The three components of the geomagnetic field data output by the triaxial magnetometer are continuously collected. The spatial angle between the geomagnetic field vectors at two adjacent sampling times is divided by the sampling time interval to obtain the angular velocity of the geomagnetic field vector rotation. A complete walking cycle is selected as the phase analysis window. The cross-correlation coefficients between the output sequences of any two axes of the triaxial magnetometer within the phase analysis window are calculated to obtain a set of triaxial correlation coefficients containing three cross-correlation coefficients. Rotational anomalies are determined by using the rotational angular velocity of the geomagnetic field vector, and correlation anomalies are determined by using the triaxial correlation coefficient set. When both rotational and correlation anomalies are found in the judgment results, the current area is determined to be a geomagnetic anomaly zone. In response to the determination that the current area is a geomagnetic anomaly zone, the azimuth output of the mechanical compass and electronic compass is forcibly blocked, and an alarm is issued; Extract the accurate magnetic azimuth and the last valid satellite positioning coordinates recorded at the last moment before the geomagnetic anomaly zone was identified. Using the accurate magnetic azimuth as the initial heading reference, start the gyroscope and accelerometer in the inertial measurement unit to perform dead reckoning and obtain continuous relative heading angles. By acquiring continuous surface images, the optical flow features between adjacent image frames are extracted using the visual odometry method to calculate the heading change, or the atmospheric polarization mode is obtained by activating the sky polarization sensor to calculate the absolute heading angle. The heading change or absolute heading angle is used to perform periodic drift correction on the continuous relative heading angle, and a corrected navigation heading that does not depend on the geomagnetic field is output. Among them, performing rotation anomaly detection and related anomaly detection refers to: The rotational angular velocity of the geomagnetic field vector is compared with the first historical baseline, which is the moving average of the rotational angular velocity of the geomagnetic field vector over a continuous period before entering the current region. When the rotational angular velocity of the geomagnetic field vector exceeds a predetermined multiple of the first historical baseline, it is marked as a rotational anomaly. The triaxial correlation coefficient set is compared with the second historical baseline, which is the moving average of the triaxial correlation coefficient set over a continuous period before entering the current region. When the minimum value in the triaxial correlation coefficient set is lower than a predetermined proportion of the corresponding minimum value in the second historical baseline, it is marked as a correlation anomaly.

2. The navigation method for field geological survey according to claim 1, wherein, The method for extracting the complete walking cycle is as follows: The vertical acceleration value output by the accelerometer in the inertial measurement unit is continuously monitored. When the vertical acceleration value crosses the zero point twice in a row and the time interval between the two crossings is within the preset normal walking frequency range, the time interval between the two zero point crossings is determined as a complete walking cycle. When step frequency cannot be detected, a fixed duration is used as a complete walking cycle. The cross-correlation coefficient is calculated as follows: for any two axis output sequences, multiply the corresponding data points of the two sequences within the same phase analysis window, sum them, and then divide by the product of the square root of each sequence to obtain the cross-correlation coefficient.

3. The navigation method for field geological survey according to claim 2, wherein, The specific implementation method of dead reckoning is as follows: The accurate magnetic azimuth recorded at the last moment before entering the geomagnetic anomaly zone is used as the initial heading angle, and the last effective satellite positioning coordinates are used as the initial position coordinates. By integrating the three-axis angular velocity output by the gyroscope over time, the change in heading angle relative to the initial heading angle at each moment can be obtained. Gait detection is performed using the triaxial acceleration output from the accelerometer. When the vertical acceleration value crosses the zero point twice consecutively, it is determined as a complete step. The step length is estimated based on the composite amplitude of the acceleration, and the step length is proportional to the fourth root of the acceleration amplitude. Multiply the step size by the cosine and sine of the current heading angle and add them to the position coordinates to obtain continuous position calculation results; After each predetermined number of steps is calculated, the angular velocity drift of the integral is corrected once using the change in heading calculated by the visual odometry method or the sky polarization light sensor. The correction amount is the average difference between the gyroscope integral heading and the visual heading within the most recent predetermined number of steps.

4. A navigation method suitable for field geological surveys according to claim 3, characterized in that, The specific steps for extracting optical flow features between adjacent image frames using the visual odometry method are as follows: The continuous surface images captured by the camera unit are converted into grayscale images, and corner features are quickly extracted from each grayscale image; The pyramid optical flow method is used to track the pixel displacement of corner features between two adjacent frames; Median filtering is applied to the displacement of all successfully tracked corner pixels to remove outliers whose displacement vector magnitude exceeds a predetermined multiple of the median standard deviation. The remaining effective pixel displacements are converted into unit pixel movement vectors, and then the heading change of the carrier in the horizontal plane is deduced by using the intrinsic parameter matrix of the camera unit and the assumed height of the carrier. When the number of effective corner points is lower than the predetermined number, the visual odometry method is deemed to be ineffective, and the system automatically switches to the sky polarization light sensor for heading calculation.

5. A navigation method suitable for field geological surveys according to claim 4, characterized in that, The specific steps for the sky polarization sensor to acquire the atmospheric polarization pattern and calculate the absolute heading angle are as follows: The atmospheric polarization degree and polarization azimuth angle in three different directions in the sky are measured using a time-division or pixel-division method. The three directions are selected from the zenith direction and the direction of the sun at predetermined angles on both sides. According to the Rayleigh scattering model, there is a geometric constraint relationship between the polarization azimuth angles measured in the three directions and the position of the sun. The azimuth angle of the sun in the horizontal coordinate system can be deduced by solving the constraint equations. The solar azimuth angle is compared with the theoretical solar azimuth angle provided by the internal clock of the device. The difference between the two is the absolute heading angle between the current orientation of the device and geographic north. When a situation arises where it is impossible to measure three valid directions simultaneously due to obstruction, the most recently successfully calculated absolute heading angle is automatically saved and used as the benchmark for subsequent dead reckoning until a valid measurement is obtained again.

6. A navigation method suitable for field geological surveys according to claim 5, characterized in that, The real-time calculation of the rotational angular velocity of the geomagnetic field vector and the identification of rotational anomalies are carried out using the following steps: A three-component geomagnetic field data sequence was established, with each sampling point containing the northward component, the eastward component, and the vertical component; For two geomagnetic field vectors at two adjacent sampling times, calculate their spatial angle. The spatial angle is equal to the dot product of the two vectors divided by the product of the magnitudes of the two vectors and then taking the inverse cosine. Divide the spatial angle by the sampling time interval to obtain the instantaneous rotational angular velocity; Calculate the average instantaneous rotational angular velocity within the sliding window, with the sliding window length set to a preset value; Calculate the coefficient of variation of the rotational angular velocity within the sliding window. The coefficient of variation is equal to the standard deviation of the angular velocity within the window divided by the average value. Calculate the average value of the coefficient of variation over a continuous period before entering the current region, and use it as the historical baseline coefficient of variation; When the coefficient of variation of the historical baseline is less than the preset variation threshold, it is marked as a rotation anomaly when the angular velocity of the geomagnetic field vector rotation exceeds a predetermined multiple of the first historical baseline. When the coefficient of variation of the historical baseline is greater than or equal to the preset variation threshold, and the angular velocity of the geomagnetic field vector rotation exceeds a predetermined multiple of the first historical baseline and the current coefficient of variation exceeds a predetermined multiple of the coefficient of variation of the historical baseline, it is marked as a rotation anomaly. When a predetermined number of consecutive sampling times are all marked as rotational anomalies, the final rotational anomaly marker is output.

7. A navigation method suitable for field geological surveys according to claim 6, characterized in that, The following steps are used to determine the correlation anomalies in the triaxial correlation coefficient set: The first step is to denote the northward component sequence, eastward component sequence, and vertical component sequence of the geomagnetic field data collected within a complete walking cycle as the first sequence, the second sequence, and the third sequence, respectively. The length of each sequence is equal to the number of sampling points within that walking cycle. The second step is to calculate the cross-correlation coefficient between the first and second sequences. The cross-correlation coefficient is equal to the product of the covariance of the two sequences and their respective standard deviations. Similarly, calculate the cross-correlation coefficients between the first and third sequences, and between the second and third sequences, to obtain the set of triaxial correlation coefficients. The third step is to calculate the statistical characteristics of the correlation coefficients of the number of complete walking cycles before entering the current area, and record the minimum value of the triaxial correlation coefficient set in each cycle to form a historical minimum value sequence. The fourth step is to calculate the moving median and moving interquartile range of the historical minimum value sequence; Fifth, calculate the minimum value in the set of triaxial correlation coefficients for the current walking cycle; Step 6: When the current minimum value is less than the sliding median minus 1.5 times the sliding interquartile range and the current minimum value is less than the preset lower limit of the correlation coefficient, it is determined that a related abnormality has occurred in the current walking cycle. When the above conditions are met for a number of consecutive preset walking cycles, the final relevant anomaly flag is output.

Images