Surveying and mapping positioning method based on multi-sensor fusion

Abstract

The application belongs to the technical field of surveying and mapping, and particularly relates to a surveying and mapping positioning method based on multi-sensor fusion, which comprises the following steps: acquiring multi-sensor observation data, solving a floating-point ambiguity vector and a variance covariance matrix; extracting a sick eigenvalue through eigenvalue decomposition, constructing a regularization transformation operator by using a remaining stable eigenvector, and projecting data to a stable subspace; performing LDLT decomposition on the subspace matrix and selecting an optimal integer multiplier through discrete trial to complete decorrelation; setting a boundary based on a chi-square distribution, and preferentially performing node search along a dimension with a minimum sum of covariance absolute values to generate a candidate set; projecting the candidate set back to the subspace to extract an optimal solution and solve a three-dimensional positioning coordinate. The application effectively avoids blind traversal and greatly improves the success rate of fixing integer ambiguity and the solving efficiency under complex environment through isolation of sick eigenvalues, optimal multiplier decorrelation and preferential dimension reduction search.

Description

Technical Field

[0001] This invention relates to the field of surveying and mapping technology. More specifically, this invention relates to a surveying and positioning method based on multi-sensor fusion. Background Technology

[0002] In the field of high-precision mapping and positioning within global navigation satellite systems, carrier phase observation is a crucial method for achieving centimeter-level and even millimeter-level positioning due to its high accuracy. The carrier phase observation equation contains unknown integer ambiguity parameters, and their fixed velocity and accuracy directly affect the reliability and timeliness of high-precision mapping and positioning. In complex mapping environments, receivers are susceptible to factors such as signal obstruction, multipath effects, atmospheric delay residues, and poor observation geometry, causing the normal matrix of the observation equation to exhibit varying degrees of ill-conditioned characteristics. This ill-conditioned structure reduces the accuracy of floating-point ambiguity resolution and creates strong statistical correlations in the dimensions of the floating-point ambiguity variance-covariance matrix. This leads to the formation of a flattened, stretched ellipsoid in the ambiguity search domain, increasing the complexity and time consumption of integer searches, and may even misclassify erroneous solutions as optimal solutions, affecting the stability and operational efficiency of the mapping and positioning system.

[0003] To reduce ambiguity correlation and improve search efficiency, ambiguity decorrelation algorithms are typically used to reconstruct and circularize the search space, accelerating high-precision positioning coordinate calculation. However, when dealing with high-dimensional covariance matrices with ill-conditioned characteristics, existing ambiguity resolution and decorrelation methods often struggle to balance decorrelation effectiveness with numerical stability of the transformation process. Due to a lack of effective identification and isolation of ill-conditioned components caused by adverse observation environments, conventional decorrelation transformations easily amplify potential observation noise and errors, disrupting the correct ambiguity search boundary and leading to candidate set distortion. Furthermore, existing search processes often employ fixed traversal logic, failing to fully consider the independence differences of parameters across dimensions, easily wasting computational resources and reducing the success rate of ambiguity fixation. Therefore, a ambiguity decorrelation method that can isolate ill-conditioned interference, optimize the search path, and maintain numerical stability is urgently needed to improve the accuracy and efficiency of high-precision mapping and positioning. Summary of the Invention

[0004] To address the technical problem that existing technologies exhibit ill-conditioned characteristics of the covariance matrix under complex observation environments, leading to strong correlations and excessive stretching of the ambiguity search space, resulting in large search computations and a high risk of misjudging the optimal solution, thereby reducing the stability and efficiency of the mapping and positioning system, this invention proposes a mapping and positioning method based on multi-sensor fusion.

[0005] This invention provides a mapping and positioning method based on multi-sensor fusion, comprising: S1, acquiring multi-sensor observation data and calculating the floating-point ambiguity vector and variance-covariance matrix; S2, performing eigenvalue decomposition on the variance-covariance matrix and extracting ill-conditioned eigenvalues ​​greater than a judgment threshold; S3, constructing a regularization transformation operator from the eigenvectors corresponding to the remaining stable eigenvalues, and projecting the floating-point ambiguity vector and the variance-covariance matrix to a stable subspace using the regularization transformation operator; S4, performing LDLT decomposition on the variance-covariance matrix in the stable subspace to obtain a lower triangular matrix, and traversing the off-diagonal elements of the lower triangular matrix and substituting them into a preset set of candidate multipliers. In step S5, discrete trial calculations are performed, and the integer that minimizes the sum of squares of local off-diagonal elements is selected as the optimal integer multiplier to complete decorrelation; S6, the chi-square distribution quantile value is determined according to the dimension of the stable subspace as the boundary square radius of the initial ellipsoid search space, the sum of the absolute values ​​of the covariance between each unfixed dimension and other unfixed dimensions is calculated, and node search is preferentially performed along the dimension with the smallest sum of absolute values ​​of the covariance to generate full-dimensional integer candidate solutions to form an integer ambiguity candidate set; S7, the integer ambiguity candidate set is projected back to the stable subspace through the regularization transformation operator to calculate the quadratic norm, the optimal solution corresponding to the minimum quadratic norm is extracted, and the three-dimensional positioning coordinates are calculated using the optimal solution.

[0006] This invention isolates ill-conditioned features through eigenvalue decomposition and constructs a regularized transformation operator. This not only eliminates the destruction of the search boundary by high-variance directions and improves the numerical stability of the transformation, but also selects the optimal multiplier through discrete trial calculations to complete the decorrelation operation, greatly reducing the correlation between parameters. At the same time, combined with dimensionality reduction evaluation based on the sum of absolute covariance values, it prioritizes locking the dimension with the strongest independence, avoiding blind traversal of the grid search space. Thus, even in complex observation environments such as severe multipath propagation and occlusion, it can still quickly and accurately search for and lock the optimal solution, significantly improving the success rate of integer ambiguity fixation and the efficiency of three-dimensional coordinate calculation in high-precision multi-source fusion mapping and positioning systems.

[0007] Preferably, the step of acquiring multi-sensor observation data and calculating the floating-point ambiguity vector and variance covariance matrix includes: extracting carrier phase data, pseudorange data, and auxiliary observation data from the multi-sensor observation data to construct a joint observation model; calling the Kalman filter algorithm to perform state parameter estimation based on the joint observation model; extracting ambiguity parameter blocks from the state vector to form the floating-point ambiguity vector; and extracting the row and column submatrices corresponding to the ambiguity parameters in the filter covariance matrix as the variance covariance matrix.

[0008] Preferably, the step of extracting ill-conditioned feature values ​​that are greater than a judgment threshold includes: arranging all feature values ​​obtained by eigenvalue decomposition of the variance-covariance matrix in descending order of numerical value to obtain a descending feature value sequence; calculating the ratio between adjacent feature values ​​in the descending feature value sequence to obtain an adjacent ratio sequence; detecting the mutation point position in the adjacent ratio sequence where the ratio is greater than a preset mutation threshold, and setting the feature value in the descending feature value sequence that is before the mutation point position as the judgment threshold.

[0009] This invention dynamically determines the judgment threshold by calculating the adjacent ratio of the descending eigenvalue sequence and detecting abrupt change points. It can accurately identify eigenvalues ​​that cause spectral discontinuities due to the influence of weak observation direction, effectively isolate deep pathological components that cause excessive stretching of the search space, and ensure the absolute stability of the subsequent projection space.

[0010] Preferably, constructing a regularization transformation operator from the eigenvectors corresponding to the remaining stable eigenvalues ​​includes: extracting the eigenvectors corresponding to the remaining stable eigenvalues ​​and concatenating them in order to construct a basic projection matrix; using a constant diagonal matrix of a preset dimension as a weight matrix; and multiplying the weight matrix with the basic projection matrix to generate the regularization transformation operator.

[0011] Preferably, the step of traversing the off-diagonal elements of the lower triangular matrix and substituting them into a preset set of candidate multipliers for discrete trial calculations, selecting the integer that minimizes the sum of squares of the local off-diagonal elements as the optimal integer multiplier, includes: calculating the nearest integer of the off-diagonal elements to be reduced; constructing the set of candidate multipliers containing zero multipliers and non-zero multipliers with the nearest integer as the center and a set of candidate radii; substituting each candidate multiplier into the column element update formula of the lower triangular matrix for trial calculations, calculating the sum of squares of all off-diagonal elements in the column after each trial calculation, and extracting the candidate multiplier that minimizes the sum of squares as the optimal integer multiplier.

[0012] This invention constructs a set of candidate multipliers containing the intervals around the nearest integer and performs discrete optimization by minimizing the sum of squares of local off-diagonal elements. This can minimize the absolute magnitude of the column elements of the lower triangular matrix, effectively reduce the coupling between parameters of different dimensions, and thus improve the decorrelation effect of the transformed search space.

[0013] Preferably, the step of prioritizing node search along the dimension with the smallest sum of absolute covariance values ​​includes: determining the dimension with the smallest sum of absolute covariance values ​​as the priority fixed target, calculating the nearest integer of the floating-point fuzziness estimate corresponding to the priority fixed target; using the nearest integer as the center, alternately selecting integer candidate values ​​on both sides of the center in order from near to far from the center to generate a search sequence, and performing node search based on the search sequence.

[0014] This invention, by determining the most independent priority fixed target and then using a zigzag search sequence that alternates between near and far integers as the center, can quickly hit the core region with the highest probability density, greatly reducing the number of invalid branches and iteration time when searching for the optimal integer solution.

[0015] Preferably, the step of projecting the integer ambiguity candidate set back to the stable subspace using the regularization transform operator to calculate the quadratic norm and extracting the optimal solution corresponding to the minimum quadratic norm includes: projecting the full-dimensional integer candidate solutions onto the stable subspace using the regularization transform operator to obtain candidate vectors; calculating the difference vector between the candidate vectors and the floating-point ambiguity vectors in the stable subspace; calculating the squared Mahalanobis distance of the difference vectors using the inverse matrix of the variance-covariance matrix in the stable subspace as the quadratic norm; extracting the candidate solution corresponding to the minimum quadratic norm as the optimal solution, and extracting the candidate solution corresponding to the second smallest quadratic norm as the suboptimal solution; and confirming that the extracted optimal solution passes the reliability test when the ratio obtained by dividing the quadratic norm of the suboptimal solution by the quadratic norm of the optimal solution is greater than a preset threshold.

[0016] Preferably, the step of extracting ill-conditioned feature values ​​that are greater than the judgment threshold includes: using the K-means clustering algorithm to divide all feature values ​​obtained by eigenvalue decomposition of the variance-covariance matrix into principal component clusters and ill-conditioned component clusters; and extracting the minimum feature value within the ill-conditioned component cluster as the judgment threshold.

[0017] Preferably, traversing the off-diagonal elements of the lower triangular matrix includes: scanning the off-diagonal elements of the lower triangular matrix from the last row upwards, and from left to right within each row; performing discrete trial calculations and selecting the optimal integer multiplier only for off-diagonal elements whose absolute values ​​are greater than a preset value, until the absolute values ​​of all off-diagonal elements of the lower triangular matrix are not greater than the preset value.

[0018] Preferably, the step of calculating the three-dimensional positioning coordinates using the optimal solution includes: determining the integer ambiguity solution by combining the confirmed optimal solution with the generalized inverse mapping of the regularization transform operator back to the original ambiguity domain; substituting the integer ambiguity solution as a fixed ambiguity constraint value back into the state update equation set to reconstruct the measurement matrix; and performing forward filtering and smoothing processing to calculate the three-dimensional positioning coordinates.

[0019] This invention uses the generalized inverse mapping and verified integer ambiguity solution as an absolute constraint to perform smoothing filtering on the state update equation, effectively eliminating the error accumulation caused by floating-point truncation, so that the accuracy of the calculated three-dimensional positioning coordinates can be stably converged to the centimeter or even millimeter level.

[0020] The beneficial effects of this invention are as follows:

[0021] This invention identifies and isolates ill-conditioned structural components by performing eigenvalue decomposition on the variance-covariance matrix and setting an adaptive threshold based on the eigenvalue distribution. It also improves numerical computation stability by projecting ambiguity onto a stable subspace using a regularization transformation operator. During decorrelation, it reduces the correlation of ambiguity parameters and improves the search space shape by optimizing the integer multiplier selection strategy to handle off-diagonal elements.

[0022] Furthermore, during the ellipsoid search phase, the dimension with the strongest independence is evaluated and prioritized for fixing in real time, suppressing blind expansion of search branches, shortening search time, and improving the fixing success rate. By combining quadratic norm ratio testing, generalized inverse mapping reference, and consistency verification of the original ambiguity domain, reliable transformation of the optimal solution is ensured, enabling rapid and accurate fixing of integer ambiguities even in complex observation environments, and obtaining high-precision and high-reliability mapping and positioning coordinates. Attached Figure Description

[0023] Figure 1 This is a flowchart of the mapping and positioning method based on multi-sensor fusion in this invention; Figure 2 This is a schematic diagram of integer multiplier optimization in this invention; Figure 3 This is a schematic diagram of the non-fixed-dimensional evaluation in this invention. Detailed Implementation

[0024] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0025] This invention discloses a mapping and positioning method based on multi-sensor fusion, referring to... Figure 1 This includes steps S1-S6: S1. Solve for the ambiguity vector and covariance matrix.

[0026] In an optional embodiment, carrier phase observation data is acquired, and the floating-point ambiguity vector and the corresponding variance-covariance matrix are calculated using a least squares or filtering algorithm.

[0027] The system acquires internationally recognized standard format files output from the front-end receiver and extracts carrier phase and pseudorange observations for each observation satellite from these files. Simultaneously, it acquires auxiliary observation data output from the inertial measurement unit (IMU), electronic compass, or ranging sensor. This auxiliary observation data serves as state prediction constraints or weak observation supplementary constraints. The auxiliary observation data, along with carrier phase and pseudorange observation data, are input into the least squares or filtering estimation process. The acceleration and angular velocity output from the IMU are used to construct the state prediction equations for position, velocity, and attitude. The heading angle output from the electronic compass or the distance observation from the ranging sensor is used as a measurement update term and, together with the GNSS observation equations, forms a joint observation model. The system calculates the satellite's three-dimensional position and clock bias using broadcast ephemeris data. A dual-frequency ionospheric-free combined model is employed to eliminate the effects of first-order ionospheric delay. Tropospheric delay errors are addressed using tropospheric empirical models, projection functions, or tropospheric residual parameter estimation methods. Tropospheric empirical models include the Saastamoinen model, and projection functions include the Niell mapping function or the global mapping function, thereby establishing a double-difference observation model. The state parameter estimation process is performed by invoking the Extended Kalman Filter algorithm or using the Weighted Least Squares algorithm. After updating the time using the state of the previous epoch, the measurement is updated based on the observations of the current epoch. After convergence over multiple epochs, the ambiguity parameter blocks in the state vector are extracted to form a floating-point ambiguity vector. The row and column submatrices corresponding to the ambiguity parameters in the filter covariance matrix are extracted, and these submatrices are used as the variance-covariance matrix of the floating-point ambiguity vector.

[0028] When using a dual-frequency, ionospheric-free combined model, ionospheric delay is no longer considered an independent parameter in the degree of freedom calculation. If a non-combined model or a non-combined multi-frequency model is used, ionospheric delay can be included as a parameter to be estimated. This embodiment uses an ionospheric-free combined model; therefore, the ionospheric delay parameter is no longer included separately in the subsequent degree of freedom calculation.

[0029] S2. Decompose the matrix to identify ill-conditioned eigenvalues.

[0030] In an optional embodiment, eigenvalue decomposition is performed on the variance-covariance matrix, and a threshold is adaptively set according to the distribution characteristics of the eigenvalues ​​to identify eigenvalues ​​associated with high variance and weakly constrained ill-conditioned structures.

[0031] The extracted variance-covariance matrix is ​​decomposed into eigenvalues ​​to obtain all non-negative eigenvalues ​​and their corresponding orthogonal eigenvector matrices. Since this embodiment analyzes the variance-covariance matrix of floating-point ambiguity, larger eigenvalues ​​indicate greater variance and weaker observation constraints in the corresponding feature direction; this physical meaning corresponds to the meaning of weak constraints corresponding to excessively small eigenvalues ​​in the normal equation matrix and should not be confused. All obtained eigenvalues ​​are sorted in descending order of numerical value. Using the K-means clustering algorithm, all eigenvalues ​​are divided into principal component clusters representing stable models and ill-conditioned component clusters representing high-variance, weakly constrained ill-conditioned structures based on one-dimensional Euclidean distance. The minimum eigenvalue within each ill-conditioned component cluster is extracted as an adaptive decision threshold. All decomposed eigenvalues ​​are iterated over; if an eigenvalue in a certain dimension is greater than or equal to the decision threshold, the dimension containing that eigenvalue is identified as an ill-conditioned structure-related eigenvalue causing high-variance stretching in the search space.

[0032] In one implementation, eigenvalue decomposition is performed on the variance-covariance matrix, and a threshold is adaptively set based on the distribution characteristics of the eigenvalues ​​to identify eigenvalues ​​associated with high-variance, weakly constrained ill-conditioned structures. Specifically, all eigenvalues ​​obtained from the eigenvalue decomposition of the variance-covariance matrix are arranged in descending order of numerical value to obtain a descending eigenvalue sequence; the ratio between adjacent eigenvalues ​​in the descending eigenvalue sequence is calculated to obtain an adjacent ratio sequence; the positions of abrupt change points in the adjacent ratio sequence are detected, and abrupt change points with ratios greater than a preset abrupt change threshold are identified; based on the position of the abrupt change point, the eigenvalue preceding the abrupt change point in the original descending eigenvalue sequence is determined, and this eigenvalue is set as the ill-conditioning threshold for distinguishing high-variance, weakly constrained ill-conditioned structures; all eigenvalues ​​in the descending eigenvalue sequence that are greater than or equal to the ill-conditioning threshold are identified as eigenvalues ​​associated with high-variance, weakly constrained ill-conditioned structures; if no abrupt change point with a ratio greater than the preset abrupt change threshold is detected, it is determined that there are no high-variance, weakly constrained ill-conditioned eigenvalues ​​to be removed from the current variance-covariance matrix, and the eigenvectors corresponding to all eigenvalues ​​are used to construct the basic projection matrix.

[0033] Assuming a multi-system, multi-frequency GNSS observation environment, the floating-point ambiguity variance-covariance matrix calculated using the least squares algorithm has a dimension of 20×20. Eigenvalue decomposition is performed on the covariance matrix to extract 20 eigenvalues, which are then arranged in descending order to form an eigenvalue sequence. Principal component information of the matrix is ​​extracted by arranging the eigenvalues ​​by numerical value, and the eigenvalue sequence is constructed based on this.

[0034] in, The sequence of eigenvalues ​​formed in descending order; to These are the feature values ​​arranged in descending order of their numerical values.

[0035] For example, the first five eigenvalues ​​are significantly biased by the weak observation direction, among which... , , , , ,and The drop to 0.08 indicates same There are obvious spectral discontinuities between them. The preset mutation threshold range for adjacent ratios is set to... to The preset mutation threshold range is obtained by statistical analysis and calibration of the discontinuous distribution of characteristic value spectra of historical multi-system, multi-frequency GNSS observation data under different weak observation environments. In this embodiment, it is set to... .

[0036] Ill-conditioned spectrum discontinuities are identified by the change in the ratio of adjacent eigenvalues, and the ratio of adjacent eigenvalues ​​is constructed based on this.

[0037] in, For the first The ratio of adjacent feature values; The first in the descending eigenvalue sequence One eigenvalue; The first in the descending eigenvalue sequence One eigenvalue; is the index of the eigenvalue.

[0038] When detected At that time, identification The location of the mutation point is determined, and the values ​​preceding the mutation point in the original descending eigenvalue sequence are identified. This is the ill-conditioning threshold used to distinguish ill-conditioning structures with high variance and weak constraints. This threshold is used to determine the ill-conditioning of sequences. to These five larger eigenvalues ​​are eigenvalues ​​related to the ill-conditioned structure with high variance and weak constraints, which solves the problem of excessive stretching of the search space caused by the high variance direction of the covariance matrix under weak observation conditions.

[0039] S3. Construct operators and project them onto a stable subspace.

[0040] In an optional embodiment, a regularization transformation operator is constructed using the eigenvectors corresponding to eigenvalues ​​not identified as being related to ill-conditioned structures, projecting the floating-point ambiguity vector and variance-covariance matrix onto a stable subspace.

[0041] From the eigenvector matrix, identify all eigenvector columns corresponding to eigenvalues ​​associated with high-variance, weakly constrained ill-conditioned structures. For eigenvectors corresponding to healthy eigenvalues ​​not identified as associated with ill-conditioned structures, concatenate and combine them according to their original dimensional order to generate a truncated base projection matrix. To avoid real-number projection disrupting the integer lattice structure of integer ambiguities, this embodiment uses the base projection matrix to construct a stable evaluation subspace and a covariance regularization matrix, instead of directly using the coordinates after real-number projection as independent integer search variables for integer ambiguities. The integer search variables are still represented by the original ambiguity integer vectors or the single-mode integer transformed coordinates formed during the subsequent LDLT decorrelation process.

[0042] Regarding the matrix dimension setting, the spatial mapping relationship is represented using matrices, and the dimensional constraint formula for the basic projection matrix is ​​constructed based on this:

[0043] in, The basic projection matrix; For the real number field; To stabilize the subspace dimension; This represents the original fuzziness dimension.

[0044] The general matrix multiplication algorithm is invoked, and the original floating-point ambiguity vector is multiplied by the basic projection matrix to complete the original parameter space transformation. The mapping of the dimensionality-stable evaluation subspace yields the dimensionality-reduced floating-point fuzzy vector. A congruent transformation of the covariance matrix is ​​performed by left-multiplying the original variance-covariance matrix by the fundamental projection matrix and right-multiplying it by the transpose of the fundamental projection matrix, resulting in a full-rank variance-covariance matrix within the stable subspace. This congruent transformation is referred to as a congruent transformation or projection transformation, not a similarity-equivalence transformation.

[0045] In one implementation, a regularization transformation operator is constructed using eigenvectors corresponding to eigenvalues ​​not identified as being related to the high-variance, weakly constrained ill-conditioned structure. Specifically, eigenvectors corresponding to all eigenvalues ​​not identified as being related to the high-variance, weakly constrained ill-conditioned structure are extracted; the extracted eigenvectors are concatenated in their original order to construct a basic projection matrix; a constant diagonal matrix of a preset dimension is used as the weight matrix, and the weight matrix is ​​multiplied by the basic projection matrix to generate a regularization transformation operator with stable projection and dimensionality reduction capabilities.

[0046] Based on the aforementioned identification results, after removing the feature values ​​associated with pathological structures, the remaining features are extracted. The corresponding feature vectors of the unidentified ill-conditioned feature values. For example, in the original 20-dimensional fuzzy space, if 5 high-variance, weakly constrained ill-conditioned feature values ​​are removed and the remaining features are... For any stable feature values ​​not identified as related to pathological structures, extract the corresponding 15 feature vectors with a column dimension of 20×1. Concatenate these 15 feature vectors column-wise into a 20×15 intermediate matrix, following their original order of eigenvalues ​​from largest to smallest. Construct a projection basis using the transpose of the orthogonal vector group, and then build the basic projection matrix based on this basis.

[0047] in, The basic projection matrix; For the reason An intermediate matrix formed by concatenating feature vectors that were not identified as pathological features; This is the identifier for the matrix transpose operation.

[0048] To maintain the matrix structure of stable feature directions and facilitate subsequent expansion of different weight strategies, a 15×15 constant diagonal matrix is ​​input as the weight matrix. The value of the preset dimension is determined by analyzing the number of stable eigenvalues ​​without spectral discontinuities in the original ambiguity space. In this embodiment, the main diagonal elements of the constant diagonal matrix are uniformly set to 1, and the value of 1 is determined by maintaining the intrinsic condition of unchanged parameter scale before and after the transformation. The off-diagonal elements are all 0 to maintain the intrinsic structure of the stable feature directions of the covariance matrix. If signal-to-noise ratio adaptive weights are used, they should be set separately as another embodiment and should not be mixed with the constant diagonal matrix embodiment. The 15×15 weight matrix is ​​left-multiplied by the 15×20 base projection matrix, and the basis is weighted by the weight matrix. On this basis, a regularization transformation operator is constructed:

[0049] in, For regularization transformation operators; This is the weight matrix; Use the base projection matrix.

[0050] The regularization transformation operator achieves the mapping from the original high-dimensional space to a 15-dimensional stable evaluation subspace by preserving stable feature directions and suppressing the influence of high-variance ill-conditioned directions on matrix inversion and search boundaries.

[0051] In the subsequent recovery phase, the regularization transform operator Moore-Penrose generalized inverse The dimension is 20×15, which can be left-multiplied by a 15-dimensional vector in the stable evaluation subspace and mapped back to the 20-dimensional original fuzzy domain, thus avoiding matrix dimension mismatch. The Moore-Penrose generalized inverse is a generalization of the ordinary matrix inverse, used to obtain stable regressions or approximate solutions that satisfy the least-squares meaning when the matrix is ​​not square or non-invertible. For ill-conditioned directions suppressed by regularization, the true integer values ​​are not recovered solely by the generalized inverse; instead, they are completed and verified in the original fuzzy domain by combining floating-point priors, nearest-neighbor integer constraints, and the original covariance quadratic residuals.

[0052] S4. Decompose the covariance matrix and perform decorrelation.

[0053] In an optional embodiment, LDLT decomposition is performed on the variance-covariance matrix in the stable subspace, and an adaptive integer multiplier selection strategy is used to decorrelate the off-diagonal elements to be reduced.

[0054] The LDLT decomposition module in the Eigen library is invoked to decompose the full-rank variance-covariance matrix in the stable subspace into a product of a unit lower triangular matrix, a diagonal matrix, and the transpose of the unit lower triangular matrix. A decorrelation double loop process is then initiated, scanning the off-diagonal elements of the lower triangular matrix from left to right within each row, starting from the last row. When an off-diagonal element with an absolute value greater than 0.5 is encountered, a candidate multiplier set is constructed centered on the nearest integer of the element to be reduced, and discrete trial calculations are performed around this center with a preset radius. The value of 0.5 is determined using the size reduction condition in LLL lattice reduction theory. If the magnitude of the element to be reduced is large, the interval of the candidate set expands synchronously with the magnitude, ensuring that the candidate set covers integer multipliers that bring the off-diagonal elements close to zero.

[0055] Using a discrete optimization algorithm, the integer that minimizes the sum of squares of the local off-diagonal elements is searched in the candidate integer set and selected as the optimal integer multiplier. Based on this, a Gaussian elementary single-mode integer transformation matrix is ​​constructed. The unit lower triangular matrix, diagonal matrix, reduced-dimensional floating-point ambiguity vector, and cumulative single-mode integer transformation matrix are synchronously multiplied and updated until the absolute values ​​of all off-diagonal elements of the lower triangular matrix are less than or equal to 0.5, or until the preset maximum number of iterations is reached and the sum of squares of the matrix's off-diagonals no longer decreases. The loop then terminates, completing the depth decorrelation processing within the subspace. The Gaussian elementary single-mode integer transformation matrix is ​​a single-mode integer matrix and is not an orthogonal matrix; therefore, this paper does not describe the cumulative transformation matrix as a cumulative orthogonal integer transformation matrix.

[0056] In one implementation, an adaptive integer multiplier selection strategy is used to select the optimal integer multiplier for the off-diagonal element to be reduced. Specifically, a candidate multiplier set containing integers within a preset interval is defined, including zero multipliers and non-zero multipliers. For each off-diagonal element to be reduced, trial calculations are performed by substituting each integer in the candidate multiplier set. The sum of squares of the off-diagonal elements in the row or column containing the off-diagonal element is calculated after each trial calculation, and the integer that minimizes the sum of squares is selected as the optimal integer multiplier. If the optimal integer multiplier is zero or the sum of squares does not decrease after the trial calculation, the integer transformation for this round is not performed.

[0057] In the decorrelation process of generating a lower triangular matrix from the variance-covariance matrix in the stable subspace using LDLT decomposition, a candidate multiplier set is constructed to reduce the correlation between parameters of different dimensions. The candidate multiplier set is determined by preset rules, first calculating the nearest integer of the off-diagonal elements to be reduced. The elements are rounded to obtain the nearest integer, and a formula for calculating the nearest integer is constructed based on this:

[0058] in, The most recent integer; This is the rounding function; The lower triangular matrix is ​​the first... same as the first The non-diagonal elements to be reduced at the intersection of columns.

[0059] Next, we define candidate radii, which are reduced using a single nearest integer, and then construct a set of candidate multipliers based on this:

[0060] in, For the set of candidate multipliers; The most recent integer; The candidate radius is denoted as .

[0061] Set candidate radius For example, setting The value of 5 was obtained by statistically evaluating and calibrating the fluctuation range of off-diagonal elements under weak observation conditions. A candidate multiplier set is formed, and the zero multiplier is retained as a candidate for not performing an update. When the sum of squares of the off-diagonal elements corresponding to the zero multiplier is minimized, or when no non-zero multiplier can decrease the sum of squares, the current round of elementary integer transformations is skipped. For general elements, this can also degenerate into a fixed interval. The set of integers is selected, and zero multipliers are reserved as candidates for not performing an update. For elements with an absolute value greater than 5.5, the candidate interval is automatically expanded; the value of 5.5 is determined by analyzing the upper limit of the statistical tolerance of off-diagonal elements deviating from the nearest integer.

[0062] For the lower triangular matrix, the first... same as the first Any off-diagonal element of the column to be reduced For example, the loop encountered in the current iteration Substitute them into the candidate multiplier set in sequence. Elementary column transformations are performed on each candidate integer in the input array for trial calculation. During the trial calculation phase, for each substituted candidate integer... The matrix column elements will be updated. Elementary column transformations are performed using Gaussian elimination, and based on this, the update formula for the matrix column elements is constructed:

[0063] in, For the updated matrix column elements; These are the column elements of the matrix before the update. The candidate integer to be substituted; The lower triangular matrix is ​​the first... same as the first Elements at the intersection of columns; The row number where the column element is located.

[0064] Synchronous calculation after update The sum of squares of all off-diagonal elements in the column. This is achieved by iterating through the candidate multiplier set multiple times independently and comparing the sums of squares for each candidate integer. Assume the input integer... When the sum of squares of the off-diagonal elements reaches a minimum of 0.15, the integer 3 is considered a pair of elements. The optimal integer multipliers. This discrete integer optimization strategy improves the decorrelation effect of the transformed space.

[0065] Reference Figure 2 The results show the trend of the sum of squares of the corresponding columns after substituting different candidate multipliers. As the candidate multipliers change, the sum of squares exhibits a discrete parabolic trend of first decreasing and then increasing, reaching its minimum extreme point when the candidate multiplier is 3. Based on the convergence trend of the sum of squares with the candidate multiplier, it can be proven that selecting the integer that minimizes the local sum of squares through discrete trial and error as the optimal integer multiplier can reduce the absolute magnitude of the matrix column elements, thereby reducing the correlation between parameters of different dimensions and improving the decorrelation effect.

[0066] S5. Define the search space and generate a candidate set.

[0067] In an optional embodiment, an initial ellipsoidal search space is set based on the chi-square distribution in the transformed ambiguity space; during the dimension-by-dimensional search process, the covariance of all remaining dimensions with other unfixed dimensions is evaluated, the dimension with the strongest independence is fixed first, and a full-dimensional integer ambiguity candidate set is generated in the original ambiguity integer domain or the single-mode integer transformation domain.

[0068] Combining the dimensional degrees of freedom of the stable subspace and a preset 99% confidence level, the chi-square quantile value is calculated using the inverse cumulative distribution function of the chi-square distribution, and this value is used as the boundary square radius of the ellipsoidal search space. The 99% confidence level is determined by balancing the probability of containing correct integer solutions with controlling the volume of the search space. Since the quadratic form norm has been dimensionlessized using the covariance inverse matrix, the boundary square radius is no longer expressed as the ellipsoidal volume, nor is it multiplied by the residual sum of squares. If a unit weight variance factor needs to be introduced, dimensionless normalization should be performed first.

[0069] A depth-first search algorithm based on a tree data structure is employed to traverse the grid nodes. Before entering the search level of each new dimension, the elements of the covariance matrix corresponding to the currently unfixed dimension are accessed, and the sum of the absolute values ​​of the covariances between each unfixed dimension and other unfixed dimensions is calculated. The dimension with the smallest sum of absolute covariances is selected as the most independent node to be searched in the current layer, changing the traditional fixed reverse search structure. This paper no longer uses the combined criterion of maximum entropy and minimum conditional variance to avoid conflict with the criterion of minimum sum of absolute covariances. Along the integer candidate range corresponding to the dimension with the smallest sum of absolute covariances, the nodes to be searched are fixed sequentially based on the Schnorr-Euichner zigzag alternating expansion sequence. The Schnorr-Euichner zigzag alternating expansion sequence refers to a search sequence that takes the nearest integer corresponding to the current floating-point ambiguity estimate as the center and selects integer candidate values ​​alternately on both sides of the center integer in order of increasing distance from the center integer. After obtaining candidate search centers based on stable subspace constraints, a finite number of full-dimensional integer candidate solutions are generated in the original fuzzy integer domain or the full-dimensional integer domain after single-mode integer transformation. For each full-dimensional integer candidate solution, the quadratic distance after projecting the full-dimensional integer candidate solution onto the stable evaluation subspace is calculated. If the quadratic distance is less than the square radius of the boundary, the corresponding full-dimensional integer candidate solution is stored in the integer fuzzy candidate set. The candidate solutions in the integer fuzzy candidate set are full-dimensional integer candidate solutions in the original fuzzy integer domain or the single-mode integer transform domain. The stable subspace is used for projection evaluation, search pruning, and quadratic norm calculation of the full-dimensional integer candidate solutions.

[0070] In one implementation, the initial ellipsoidal search space is set based on the chi-square distribution. Specifically, carrier phase observation data of the station's positioning epoch are acquired, and the degrees of freedom of the chi-square distribution are determined according to the dimension of the stable subspace; a preset chi-square distribution confidence level is selected, and the chi-square distribution quantile values ​​under the corresponding degrees of freedom are found; the chi-square distribution quantile values ​​are used as the boundary square radius of the initial ellipsoidal search space to complete the setting of the initial ellipsoidal search space.

[0071] Since the ellipsoidal search space in this embodiment is used to constrain the quadratic norm within the stable subspace, the chi-squared distribution degrees of freedom are taken from the dimension of the stable evaluation subspace. Assuming the original fuzzy space has a dimension of 20, after removing 5 high-variance, weakly constrained ill-conditioned feature directions, the remaining stable subspace has a dimension of [missing value]. Therefore, 15 is used as the chi-square distribution degree of freedom. At a 99% confidence level, the chi-square distribution quantiles for the corresponding degrees of freedom are obtained by consulting the standard chi-square distribution statistics table or by performing the inverse operation of the cumulative distribution function. These quantiles are then used as the boundary square radius of the initial ellipsoidal search space.

[0072] In one implementation, the covariance of all remaining dimensions with other unfixed dimensions is evaluated, and the dimension with the strongest independence is fixed first. Specifically, after fixing each dimension, the covariance of each remaining dimension with all other unfixed dimensions is extracted; the sum of the absolute values ​​of the covariance of each unfixed dimension is calculated; the sum of the absolute values ​​of the covariance of each unfixed dimension is compared, and the dimension with the smallest sum of absolute values ​​of covariance is selected as the dimension with the strongest independence, and the dimension with the smallest sum of absolute values ​​of covariance is selected as the priority fixing target for the next round of search.

[0073] A conditional covariance evaluation mechanism is utilized in the dimensionality reduction process of fuzzy search. Assuming that in a 15-dimensional stable subspace, the first 3 dimensions have been searched and fixed, there are still 12 unfixed dimensions remaining in the candidate subset. , For each dimension in the candidate subset, access the current reduced variance-covariance matrix and extract all covariance elements at the intersections of the corresponding row and the corresponding columns of the other 11 unfixed dimensions. Take the absolute value of the extracted cross covariance values ​​and sum them to obtain the sum of the absolute values ​​of the covariance for the corresponding dimension.

[0074] Calculate the sum of the absolute values ​​for each of the remaining 12 dimensions. For example, calculate the dimensions... The result is 0.35, dimension The result is 1.24, dimension The result was 0.88. Sorting and comparing the 12 generated values ​​revealed the dimension. The sum of the absolute values ​​of the covariances, 0.35, is the minimum among all unfixed dimensions. In statistics, a minimum value represents the minimum value of a dimension. The linear correlation with all other parameters is weak, i.e., the dimension is relatively small. The dimension with the weakest coupling to other dimensions is suitable for priority fixing to reduce the number of subsequent search branches; therefore, the dimension is... The dimension is determined to be the most independent at present. In the next round of search, the original preset order is skipped, and the dimension is... The extracted targets are prioritized and fixed. The evaluation mechanism reduces the conditional variance coupling of deep nodes and decreases the traversal volume of the internal mesh search space.

[0075] Reference Figure 3 The figure presents the distribution trend of the sum of absolute covariances as a function of dimension indices. As can be seen from the figure, the sum of absolute covariances for different dimension indices exhibits a fluctuating pattern, reaching a global minimum at a certain dimension index. Based on the distributional differences in the sum of absolute covariances, it can be demonstrated that by calculating and comparing the sum of absolute covariances for each unfixed dimension in real time, the dimension with the weakest linear correlation to other parameters can be effectively identified. Prioritizing the fixing of the dimension with the smallest sum of absolute covariance can significantly reduce the conditional variance coupling of deep nodes, thereby reducing the traversal volume of the grid search space and improving search efficiency.

[0076] S6. Calculate the quadratic norm and solve for the coordinates.

[0077] In an optional embodiment, the integer ambiguity candidate set is transformed back to the stable subspace, and the quadratic norm of each candidate solution is calculated to identify the optimal and suboptimal solutions. The optimal solution is confirmed when the ratio of the quadratic norm of the suboptimal solution to the quadratic norm of the optimal solution is greater than a preset threshold. The preset threshold is determined by statistical analysis and calibration of the norm ratio distribution characteristics of massive historical observation data under different weak observation environments. The generalized inverse of the regularization transform operator is used to obtain the reflected reference value of the original ambiguity domain, and the integer ambiguity solution is determined by combining the consistency verification of the original ambiguity domain, which is then used to solve for high-precision positioning coordinates.

[0078] The algorithm iterates through the integer ambiguity candidate set in the memory stack, calling a general matrix operation library. When a candidate solution is located in the single-mode integer transform domain, the inverse of the cumulative single-mode integer transform matrix generated in step 4 is used to map the candidate solution back to the original integer ambiguity domain. If the candidate solution is already located in the original integer ambiguity domain, the algorithm directly proceeds to subsequent projection evaluation. Then, a regularization transform operator is used to project the full-dimensional integer candidate solution onto the stable subspace, obtaining candidate vectors that are in the same evaluation coordinate system as the floating-point ambiguity vector in the stable subspace. The residual difference vector between each candidate solution vector and the reduced-dimensional floating-point ambiguity vector is calculated. Using the inverse of the variance covariance matrix in the stable subspace as weighting coefficients, the squared Mahalanobis distance of each residual difference vector is calculated as its quadratic norm. The Mahalanobis distance is a statistical distance obtained by weighting the difference vector using the variance covariance matrix, used to characterize the deviation between the candidate solution and the floating-point solution under the relevant noise distribution. The quicksort algorithm is used to sort all quadratic form norm values ​​in ascending order. The vector corresponding to the smallest first value is selected as the optimal solution, and the vector corresponding to the second smallest value is selected as the suboptimal solution. The ratio of the quadratic form norm of the suboptimal solution to that of the optimal solution is calculated, and this ratio is compared with a predetermined empirical constant threshold 3. The value of threshold 3 is obtained by experimental verification and calibration of the statistical distance difference boundary between the suboptimal and optimal solutions in the observation grid. When the ratio is greater than the threshold, the optimal solution is deemed to have passed the reliability test.

[0079] The singular value decomposition module is invoked to calculate the Moore-Penrose generalized inverse matrix of the regularization transform operator generated in step 3. Since the regularization transform operator in this embodiment... For a 15×20 dimensional matrix, the generalized inverse of the regularization transform operator. It is a 20×15 dimensional matrix, therefore it can be derived from... The optimal solution vector in the 15-dimensional stable subspace, which has passed verification, is then mapped back from the stable evaluation subspace to the 20-dimensional original fuzziness domain to avoid dimensionality mismatch. For complement space components not directly searched in the stable subspace, instead of using a single rounding result, a finite neighborhood integer set is determined based on the standard deviation of the original floating-point fuzziness components. Multiple complement space candidate combinations are then enumerated within this finite neighborhood integer set. Centered on the original fuzziness domain reference value obtained from the generalized inverse mapping, multiple full-dimensional integer candidate solutions are generated by combining the finite neighborhood integer combinations of the complement space. These solutions are then substituted into the quadratic residual and Ratio test of the original fuzziness domain for consistency verification. For the slight floating-point truncation error generated during the pull-back mapping process, only the full-dimensional integer candidate solutions that have passed the quadratic residual and Ratio test of the original fuzziness domain are corrected for integer consistency to obtain integer fuzziness solutions. Using integer ambiguity solutions as fixed ambiguity constraints, we substitute them back into the state update equations of the extended Kalman filter to reconstruct the measurement matrix, solve for the feedback gain of the baseline position parameter covariance matrix, perform forward filtering and smoothing, and calculate high-precision 3D positioning coordinates that converge to the centimeter level.

[0080] In one implementation, the integer ambiguity candidate set is transformed back to the stable subspace, and the quadratic norm of each candidate solution is calculated to identify the optimal and suboptimal solutions. Specifically, each full-dimensional integer ambiguity candidate solution in the integer ambiguity candidate set is projected onto the stable subspace using a regularization transformation operator to obtain a stable subspace candidate vector; the difference vector between the stable subspace candidate vector and the floating-point ambiguity vector in the stable subspace is calculated; the quadratic norm of the difference vector is calculated using the inverse matrix of the variance-covariance matrix in the stable subspace; the quadratic norms of all candidate solutions are sorted in ascending order of numerical value; the candidate solution corresponding to the smallest quadratic norm in the first position is set as the optimal solution, and the candidate solution corresponding to the second smallest quadratic norm is set as the suboptimal solution; the quadratic norm corresponding to the suboptimal solution is divided by the quadratic norm corresponding to the optimal solution to obtain the norm ratio used to confirm the optimal solution.

[0081] Each full-dimensional integer ambiguity candidate solution in the integer ambiguity candidate set is extracted sequentially. A regularization transformation operator is used to project the full-dimensional integer ambiguity candidate solution onto the stable subspace, obtaining the corresponding stable subspace candidate vector. The stable subspace candidate vector is compared with the stable subspace floating-point ambiguity vector to obtain the difference vector between them. The variance-covariance matrix corresponding to the stable subspace is called, and the inverse of the variance-covariance matrix is ​​used to perform a weighted operation on the difference vector to obtain the quadratic norm used to characterize the degree of deviation of the candidate solutions. The smaller the quadratic norm, the closer the corresponding candidate solution is to the stable subspace floating-point ambiguity vector in a statistical sense.

[0082] For example, the quadratic form norm results of the first three solutions in the candidate set obtained by calculating the quadratic form norm are as follows: , , Sort all calculated quadratic form norms in ascending order to generate an increasing sequence from smallest to largest. According to the least squares objective function principle, the smaller the quadratic norm value, the closer the probability distribution of the integer solution set is to the floating-point solution of the observed epoch. Therefore, the candidate solution for integer ambiguity corresponding to the first and smallest quadratic norm in the sequence, i.e., 1.12, is set as the optimal solution. Similarly, the candidate solution for integer ambiguity with the second largest norm value, 3.7, is set as the suboptimal solution. Subsequently, a preset Ratio threshold can be used, for example, to determine... The optimal solution is confirmed through double verification and restored to the original fuzzy domain for coordinate calculation.

[0083] The experiment was conducted in a weak urban observation environment, continuously acquiring 12 hours of multi-frequency, multi-system carrier phase data at a sampling rate of 1 Hz. Multipath noise was introduced to simulate the high-variance ill-conditioned characteristics of the floating-point ambiguity covariance matrix. Three sets of comparative models were set up: the baseline group used a traditional fixed ambiguity model with conventional Gaussian rounding, decorrelation, and fixed-order search; the ablation group retained the covariance evaluation search mechanism but removed the adaptive integer multiplier selection strategy; and the complete group used the complete fixed model proposed in this paper, which includes eigenvalue regularization screening, covariance evaluation search, and integer multiplier optimization.

[0084] The experiment processed a total of 43,200 epochal data points. Results showed that the baseline group achieved a 72.5% ambiguity fixation success rate, with an average single-epoch search time of 56 ms and a first convergence time of 29 s for 3D coordinates. The ablation group saw an improved fixation success rate of 83.2%, a reduction in the average single-epoch search time to 31 ms, and a first convergence time shortened to 20 s. The complete group performed best, achieving a fixation success rate of 94.8%, an average single-epoch search time of 14 ms, and a first convergence time shortened to 8 s.

[0085] Compared to the baseline group, the ablation group shows significant improvements in search time and convergence speed, indicating that the covariance evaluation search mechanism can reduce the conditional variance coupling of deep nodes and shrink the grid search space. The complete group further improves the fixation success rate and reduces search time based on the ablation group, demonstrating that the adaptive integer multiplier trial-and-error strategy can reduce rounding error accumulation and improve the decorrelation effect of the transformed space. The modules work together to improve the problems of floating-point dedivergence and low ambiguity fixation efficiency under weak observation conditions.

[0086] It should be noted that those skilled in the art can make various modifications and improvements without departing from the inventive concept, and these all fall within the scope of protection of this invention. Therefore, the scope of protection of this patent should be determined by the appended claims.

Claims

1. A mapping and positioning method based on multi-sensor fusion, characterized in that, include: S1. Acquire multi-sensor observation data and calculate the floating-point ambiguity vector and variance-covariance matrix; S2. Perform eigenvalue decomposition on the variance-covariance matrix and extract ill-conditioned eigenvalues ​​that are greater than the judgment threshold. S3. Construct a regularization transformation operator from the eigenvectors corresponding to the remaining stable eigenvalues, and use the regularization transformation operator to project the floating-point ambiguity vector and the variance-covariance matrix to the stable subspace. S4. Perform LDLT decomposition on the variance-covariance matrix in the stable subspace to obtain a lower triangular matrix. Traverse the off-diagonal elements of the lower triangular matrix and substitute them into a preset set of candidate multipliers for discrete trial calculation. Select the integer that minimizes the sum of squares of the local off-diagonal elements as the optimal integer multiplier to complete decorrelation. S5. Determine the chi-square distribution quantile value based on the dimension of the stable subspace as the boundary square radius of the initial ellipsoid search space. Calculate the sum of the absolute values ​​of the covariances between each unfixed dimension and other unfixed dimensions. Prioritize node search along the dimension with the smallest sum of absolute values ​​of the covariances to generate a full-dimensional integer candidate solution to form an integer ambiguity candidate set. S6. Project the integer ambiguity candidate set back to the stable subspace using the regularization transformation operator to calculate the quadratic norm. Extract the optimal solution corresponding to the minimum quadratic norm and use the optimal solution to calculate the three-dimensional positioning coordinates.

2. The mapping and positioning method based on multi-sensor fusion according to claim 1, characterized in that, The step of acquiring multi-sensor observation data and calculating the floating-point ambiguity vector and variance covariance matrix includes: extracting carrier phase data, pseudorange data and auxiliary observation data from the multi-sensor observation data to construct a joint observation model; calling the Kalman filter algorithm to perform state parameter estimation based on the joint observation model, extracting ambiguity parameter blocks from the state vector to form the floating-point ambiguity vector, and extracting the row and column submatrices corresponding to the ambiguity parameters in the filter covariance matrix as the variance covariance matrix.

3. The mapping and positioning method based on multi-sensor fusion according to claim 1, characterized in that, The step of extracting ill-conditioned feature values ​​that are greater than a judgment threshold includes: arranging all feature values ​​obtained by eigenvalue decomposition of the variance-covariance matrix in descending order of numerical value to obtain a descending feature value sequence; calculating the ratio between adjacent feature values ​​in the descending feature value sequence to obtain an adjacent ratio sequence; detecting the mutation point position in the adjacent ratio sequence where the ratio is greater than a preset mutation threshold, and setting the feature value in the descending feature value sequence that is before the mutation point position as the judgment threshold.

4. The mapping and positioning method based on multi-sensor fusion according to claim 3, characterized in that, The step of constructing a regularization transformation operator from the eigenvectors corresponding to the remaining stable eigenvalues ​​includes: extracting the eigenvectors corresponding to the remaining stable eigenvalues ​​and concatenating them in order to construct a basic projection matrix; using a constant diagonal matrix of a preset dimension as a weight matrix; and multiplying the weight matrix with the basic projection matrix to generate the regularization transformation operator.

5. The mapping and positioning method based on multi-sensor fusion according to claim 1, characterized in that, The process of traversing the off-diagonal elements of the lower triangular matrix and substituting them into a preset set of candidate multipliers for discrete trial calculations, selecting the integer that minimizes the sum of squares of the local off-diagonal elements as the optimal integer multiplier, includes: calculating the nearest integer of the off-diagonal elements to be reduced; constructing a set of candidate multipliers containing zero multipliers and non-zero multipliers with the nearest integer as the center and a set of candidate radii; substituting each candidate multiplier into the column element update formula of the lower triangular matrix for trial calculations, calculating the sum of squares of all off-diagonal elements in the column after each trial calculation, and extracting the candidate multiplier that minimizes the sum of squares as the optimal integer multiplier.

6. The mapping and positioning method based on multi-sensor fusion according to claim 1 or 5, characterized in that, The step of prioritizing node search along the dimension with the smallest sum of absolute covariance values ​​includes: determining the dimension with the smallest sum of absolute covariance values ​​as the priority fixed target, calculating the nearest integer of the floating-point fuzziness estimate corresponding to the priority fixed target; using the nearest integer as the center, alternately selecting integer candidate values ​​on both sides of the center in order from near to far from the center to generate a search sequence, and performing node search based on the search sequence.

7. The mapping and positioning method based on multi-sensor fusion according to claim 1, characterized in that, The step of projecting the integer ambiguity candidate set back to the stable subspace using the regularization transform operator to calculate the quadratic norm and extracting the optimal solution corresponding to the minimum quadratic norm includes: projecting the full-dimensional integer candidate solutions onto the stable subspace using the regularization transform operator to obtain candidate vectors; calculating the difference vector between the candidate vectors and the floating-point ambiguity vectors in the stable subspace; calculating the squared Mahalanobis distance of the difference vectors using the inverse matrix of the variance-covariance matrix in the stable subspace as the quadratic norm; extracting the candidate solution corresponding to the minimum quadratic norm as the optimal solution, and extracting the candidate solution corresponding to the second smallest quadratic norm as the suboptimal solution; and confirming that the extracted optimal solution passes the reliability test when the ratio of the quadratic norm of the suboptimal solution to the quadratic norm of the optimal solution is greater than a preset threshold.

8. The mapping and positioning method based on multi-sensor fusion according to claim 1, characterized in that, The extraction of ill-conditioned feature values ​​that are greater than the judgment threshold includes: using the K-means clustering algorithm to divide all feature values ​​obtained by eigenvalue decomposition of the variance-covariance matrix into principal component clusters and ill-conditioned component clusters; and extracting the minimum feature value within the ill-conditioned component cluster as the judgment threshold.

9. The mapping and positioning method based on multi-sensor fusion according to claim 1, characterized in that, The step of traversing the off-diagonal elements of the lower triangular matrix includes: scanning the off-diagonal elements of the lower triangular matrix in an order from left to right within each row, starting from the last row of the lower triangular matrix upwards; performing discrete trial calculations and selecting the optimal integer multiplier only on off-diagonal elements whose absolute values ​​are greater than a preset value, until the absolute values ​​of all off-diagonal elements of the lower triangular matrix are not greater than the preset value.

10. The mapping and positioning method based on multi-sensor fusion according to claim 1, characterized in that, The step of calculating the three-dimensional positioning coordinates using the optimal solution includes: determining the integer ambiguity solution by combining the confirmed optimal solution with the generalized inverse mapping of the regularization transform operator back to the original ambiguity domain; substituting the integer ambiguity solution as a fixed ambiguity constraint value back into the state update equation set to reconstruct the measurement matrix; and performing forward filtering and smoothing processing to calculate the three-dimensional positioning coordinates.

Images