Centralized Beidou pseudolite indoor dynamic positioning method based on known initial point

Abstract

The invention discloses a centralized Beidou pseudo-satellite indoor dynamic positioning method based on a known initial point, which is used under the conditions of large indoor pseudo-range error and limited observation data, and comprises the following contents: 1) a receiver collects pseudo-satellite observation data from the initial point, and calculates the position of each pseudo-satellite transmitting antenna; 2) calibrating differential carrier phase integer ambiguity and relative satellite channel delay combination quantity of other satellites by taking one satellite as a reference through a known initial point; 3) storing the difference integer ambiguity and the relative satellite channel time delay combination amount, and substituting the difference integer ambiguity and the relative satellite channel time delay combination amount into a satellite observation equation corresponding to a subsequent epoch for data correction; and 4) performing position calculation by using a combined weight weighted least square method to obtain a high-precision receiver coordinate, and comparing different weighted combined positioning precisions.

Description

Technical Field

[0001] The invention belongs to the field of indoor positioning, and in particular relates to a centralized Beidou pseudo-satellite indoor dynamic positioning method based on a known initial point. Background Art

[0002] With the development of the global navigation satellite system (GNSS), satellite navigation technology has been widely used in various positioning, navigation and timing (PNT) applications. Beidou Satellite Navigation System (BDS), as a global satellite navigation system independently developed by my country, has become one of the four major global navigation systems along with GPS, GLONASS and Galileo. However, the application of traditional GNSS signals in indoor environments is severely limited. The main reason is that satellite signals will be significantly attenuated when penetrating buildings, or even completely blocked, resulting in low indoor positioning accuracy or even inability to achieve positioning. In response to the above problems, the concept of pseudo-satellite was proposed and gradually applied to indoor positioning systems. Pseudolite is a ground-based transmitting device that enhances or replaces satellite signals by simulating real satellite signals, thereby achieving positioning indoors or in other areas where GNSS signals are difficult to cover. However, in indoor environments, multipath error is a major problem that hinders the realization of pseudo-satellite positioning. Existing pseudo-satellite positioning research often ignores multipath error for simulation research and lacks practical application value. In practical applications, especially in the fields of industrial automation, warehousing logistics and intelligent robots, the demand for high-precision positioning in dynamic environments is becoming more and more urgent. Therefore, it is of great practical significance and application value to study a Beidou pseudo-satellite indoor dynamic positioning method based on a known initial point that can overcome multipath errors. Summary of the invention

[0003] In view of the above technical problems, the present invention proposes a centralized Beidou pseudo-satellite indoor dynamic positioning method based on a known initial point to overcome the shortcomings of the existing technology, which is used for indoor positioning under the condition of large pseudo-range error and limited observation data. Based on the known initial point, the reference satellite is determined, and the differential carrier phase integer ambiguity and relative satellite channel delay combination of other satellites are calibrated by inter-satellite differential. After the differential integer ambiguity and relative satellite channel delay combination are stored, they are brought into the satellite observation equation corresponding to the subsequent epoch for data correction, and finally a combined weighted least squares method is used to solve the position and obtain high-precision receiver coordinates.

[0004] To achieve the above object, the present invention comprises the following steps:

[0005] Step 1: The receiver starts to collect pseudolite observation data from the initial point and calculates the position of each pseudolite transmitting antenna;

[0006] Step 2: Using a known initial point and taking one satellite as a reference, calibrate the differential carrier phase integer ambiguity and relative satellite channel delay combination of other satellites;

[0007] Step 3: After storing the differential integer ambiguity and the relative satellite channel delay combination, bring it into the satellite carrier phase observation equation corresponding to the subsequent epoch to perform data adjustment;

[0008] Step 4: Calculate the altitude angle of each pseudo-satellite based on the latest receiver position;

[0009] Step 5: Use a weighted least square method based on the median value of the mutual difference, satellite altitude angle, and C / N0 to solve the position and obtain high-precision receiver coordinates, and compare the positioning errors of different weighted combinations.

[0010] Furthermore, in the step 1, the receiver collects pseudo-satellite observation data from the initial point and solves the position of each pseudo-satellite transmitting antenna, specifically: the UM980 receiver module of Unicore Starlink is used to collect Beidou pseudo-satellite original observation data and ephemeris data from the initial point, and the carrier phase observation value and C / N0 of each satellite are stored according to the epoch, and the transmitting antenna position of each satellite is calculated through the ephemeris file.

[0011] Furthermore, in step 2, the differential carrier phase integer ambiguity of other satellites and the relative satellite channel delay combination are calibrated with one satellite as the reference through the known initial point, specifically:

[0012] (1) Observe n (n>=5) satellites at a known initial point, list the carrier phase observation equations of n satellites, and select one of the satellites as the reference satellite. The carrier phase equations of n satellites are as follows:

[0013]

[0014] Where L represents the carrier phase observation value, with the satellite number as the superscript and the epoch number as the subscript, N represents the integer ambiguity, with the satellite number as the superscript, r is the geometric distance between the satellite and the receiver, with the satellite number as the superscript and the epoch number as the subscript, δt r and δt s are the receiver clock error and satellite clock error respectively, ε is other random error terms, λ and C are the carrier wavelength and speed of light respectively.

[0015] (2) Considering that this centralized pseudo-satellite uses the same source clock, the error generated by the clock source in the satellite clock error is the same, so the observation equation is re-expressed as follows:

[0016]

[0017] Among them, ti is the clock error generated by the clock source of the i-th epoch, T is the satellite channel delay, and the first satellite is taken as the reference satellite. Since the other random error terms of the carrier phase observation are very small, their influence is ignored here. The carrier phase observation equations of the remaining n-1 satellites are subtracted from the carrier phase observation equations of the reference satellite, and the differential carrier phase integer ambiguity of the n-1 satellites and the relative satellite channel delay combination can be obtained:

[0018]

[0019] Furthermore, in step 3, after the differential integer ambiguity and the relative satellite channel delay combination are stored, they are brought into the satellite observation equation corresponding to the subsequent epoch to perform data correction, which is specifically:

[0020] Assume that we calibrate the combination of differential integer ambiguity and relative satellite channel delay through the first epoch:

[0021]

[0022] (T 1 -T n )+λ(N 1 -N n ) is brought into the subsequent epoch carrier phase observation equation, taking the second epoch as an example, we have:

[0023]

[0024] Furthermore, in step 4, the altitude angle of each pseudo-satellite is calculated with the latest receiver position for subsequent combined weighted least squares algorithm to solve the satellite position, which is specifically as follows:

[0025] (1) If the pseudo-satellite system uses an indoor relative coordinate system, assuming that the relative coordinates of the satellite are (x s ,y s ,z s )The relative coordinates of the receiver are (x r ,y r ,z r ), the relative coordinate system altitude angle calculation method is used.

[0026] Calculate the satellite's position vector relative to the receiver:

[0027] Δx=x s -x r

[0028] Δy=y s -y r

[0029] Δz=z s -zr

[0030] Calculate the horizontal distance d xy , which is the projection length of the relative position on the horizontal plane:

[0031]

[0032] Calculate the altitude angle:

[0033]

[0034] (2) If the pseudo-satellite system uses an absolute coordinate system, assuming that the ECEF coordinates of the satellite are (x s ,y s ,z s )The ECEF coordinates of the receiver are (x r ,y r ,z r ), the geodetic coordinates are (φ, λ, h), and the absolute coordinate system altitude angle calculation method is used.

[0035] The observation vector from the user to the satellite is:

[0036]

[0037] The unit observation vector I of the satellite at the receiver coordinate point is:

[0038]

[0039] Observation vector [Δx Δy Δz] T It can be equivalently expressed as the vector [Δe Δn Δu] in the station center coordinate system with the receiver coordinate point as the origin T , and its transformation relationship is:

[0040]

[0041] Among them, the coordinate transformation matrix S is:

[0042]

[0043] With the satellite observation vector at the receiver [Δe Δn Δu] T , the satellite's elevation angle θ relative to the receiver is:

[0044]

[0045] Furthermore, in step 5, a combined weighted least square method based on the median value of the mutual difference, the satellite altitude angle, and the carrier-to-noise ratio is used to perform position calculation to obtain high-precision receiver coordinates, which is specifically:

[0046] (1) The observation equation after data adjustment is:

[0047]

[0048] Use the adjusted observation term on the left side of the equal sign in the above formula Indicated by TN, n is the satellite number, i is the epoch number, and the items on the right except the geometric distance between the satellite and the receiver are used as residual items. i Expressed as, and the geometric distance term is expanded, that is:

[0049]

[0050] in represents the receiver coordinates at the i-th epoch, Indicates the coordinates of the nth satellite, TN i is the residual term of the i-th epoch.

[0051] (2) In the first epoch, the differential carrier phase integer ambiguity and the relative satellite channel delay are combined. The receiver coordinates of the first epoch are the known initial coordinates. Starting from the second epoch, the combined weighted least squares algorithm is used to solve the receiver coordinates and residual terms. The receiver initial coordinates (x0, y0, z0) are set to the receiver coordinates of the previous epoch, and the initial residual term TN0 is set to 0. k represents the number of iterations currently being performed, and k-1 represents the number of iterations currently completed. When k = 1, the initial coordinates and TN0 are used as the initial values ​​of the iteration, otherwise the results of k-1 iterations are used as the initial values ​​of the kth iteration. For the kth iteration, the observation equation is (x k-1 ,y k-1 ,z k-1 ), TN k-1 Linearize at . Take the nth observation equation as an example to linearize it, ignoring the epoch number, and find the nonlinear term r n x r The partial derivative of , we get:

[0052]

[0053] Similarly, we can find r n For r , z r The partial derivative of , we get:

[0054]

[0055] Note that the right side of the observation equation is already a linear function of TN, or the partial derivative of the left side of the equation with respect to TN is 1. The observation equation group is in (x k-1 ,y k-1 ,z k-1 ), TNk-1 The result of linearization is:

[0056]

[0057] To simplify the equation, introduce matrices G and b:

[0058]

[0059]

[0060] The linearized equation can be simplified to:

[0061]

[0062] Considering the minimum sum of squares of the differences on both sides of the equation, we can get the least squares solution of the matrix equation:

[0063]

[0064] Update the receiver's position coordinates and residual items according to the above formula:

[0065]

[0066] The iteration result is compared with a determined threshold. If the result is greater than the threshold, the iteration continues. Otherwise, the iteration can be terminated and the solution of the current nonlinear equation group is used as the position coordinates of the receiver obtained by positioning and the residual term.

[0067] The idea of ​​weighted least squares is to give different weights to the observations of each satellite. The observations with larger weights play a greater role in solving the equation. The weighted model can reduce the weight of satellites with poor observation quality in positioning. The weighted model can also suppress the multipath effect to a certain extent and provide better positioning solution results. The steps of weighted least squares positioning solution are basically the same as those of least squares positioning solution. The difference is mainly reflected in the step of solving the matrix equation. The weighted least squares method introduces the weight matrix W to solve the matrix equation. Since there is no correlation between the measurements, the weight matrix is ​​usually a diagonal matrix, as follows:

[0068]

[0069] After setting the weights corresponding to the satellite measurements, multiply each item in the simplified linearized equation by the corresponding weight coefficient, and the corresponding matrix equation is:

[0070]

[0071] The least squares method is used to solve the matrix equation and the solution of the equation is obtained:

[0072]

[0073] (3) Considering the design of the weight matrix, a combined weighted least squares method based on the median value of the mutual difference, satellite altitude angle, and C / N0 is proposed, considering the following formula:

[0074]

[0075] It is in (x k-1 ,y k-1 ,z k-1 ), TN k-1 The residual term at is the deviation between the pseudorange value corresponding to the actual position of the receiver and the corresponding pseudorange value at the linearization point. Other errors that are ignored are also included in TN and b i TN k-1 The proportion of b in the whole item is much greater than that of other items, so i It cannot accurately reflect the size of the error and deviation. i Arrange in ascending order to get b asc ,Right now Let mid = fix((n+1) / 2), where fix is ​​the rounding function, and take b asc The middle value of Subtract b from each element in the original array b mid And take the absolute value to get the new array a=[a1, a2, ..., a n ]

[0076] in:

[0077] For the same linearization point, after calculating the mutual difference, the common errors such as the residual term and receiver noise error shared by different satellites at the same epoch are eliminated. In practice, due to the influence of various errors, b is positive or negative, so we take the b in the middle of the b array mid The corresponding satellite has the largest weight, which is set to c mid =(0.55-n / 20), which changes according to the number n of satellites involved in positioning.

[0078] The remaining weight coefficients are:

[0079] The weight coefficient matrix ω1 can be constructed:

[0080]

[0081] Next, consider the weight coefficient matrix ω2 based on the altitude angle. Satellites with high altitude angles often have smaller multipath errors. The weight coefficient matrix is ​​constructed using the sine function model:

[0082]

[0083] in:

[0084]

[0085] where α i = sin(θ i ), θ i is the satellite altitude angle calculated previously.

[0086] Then construct the weight coefficient matrix ω3 based on C / N0. Often, satellites with high C / N0 have better signal quality and smaller errors. According to experiments, when the transmission power and gain are fixed to the set value, the size of C / N0 is between 40dbHZ and 60dbHZ. When C / N0 is less than 40dbHZ, the satellite is considered unavailable. Therefore, we can make:

[0087]

[0088] The weight matrix is ​​designed as:

[0089] in:

[0090]

[0091] Combining the first three weighting methods to design the combined weight, the combined weight coefficient is designed as follows:

[0092]

[0093] The combined weight matrix is:

[0094]

[0095] High-precision positioning can be achieved by taking ω as the weight matrix and bringing it into weighted least squares.

[0096] Finally, a dynamic experiment was designed, and different weighted combinations were analyzed and compared. The receiver moved on the trajectory and stayed at known point 1 and known point 2 for 60 epochs. The positioning errors in all directions of the coordinates solved by the least square method using different weighted combinations were compared and analyzed, and the average error and standard deviation were calculated.

[0097] The present invention provides a centralized Beidou pseudo-satellite indoor dynamic positioning method based on a known initial point. Based on the known initial point, by determining the reference satellite, the differential carrier phase integer ambiguity of other satellites and the relative satellite channel delay combination are calibrated by inter-satellite difference, the differential integer ambiguity and the relative satellite channel delay combination are stored, and then brought into the satellite observation equation corresponding to the subsequent epoch to perform data correction, and finally a combined weighted least square method is used to perform position solution to obtain high-precision receiver coordinates. This method can effectively solve the problem of large indoor multipath error, poor pseudo-range observation value accuracy, and the inability to fix the integer ambiguity. BRIEF DESCRIPTION OF THE DRAWINGS

[0098] Figure 1 It is a flow chart of the centralized Beidou pseudo-satellite indoor dynamic positioning method based on a known initial point of the present invention;

[0099] Figure 2 This is the x-direction positioning accuracy diagram for 60 epochs at points 1 and 2 in the dynamic experimental trajectory;

[0100] Figure 3 This is the y-direction positioning accuracy map of 60 epochs for points 1 and 2 in the dynamic experimental trajectory;

[0101] Figure 4 This is the z-direction positioning accuracy diagram for 60 epochs at points 1 and 2 in the dynamic experimental trajectory;

[0102] Figure 5 This is the mean value of the x-direction positioning error of 60 epochs at points 1 and 2 in the dynamic experimental trajectory;

[0103] Figure 6 This is the mean value of the y-direction positioning error of point 1 and point 2 in the dynamic experimental trajectory for 60 epochs;

[0104] Figure 7 This is the mean value of the z-direction positioning error of point 1 and point 2 in the dynamic experimental trajectory for 60 epochs; DETAILED DESCRIPTION

[0105] according to Figure 1 The steps shown are used to describe in detail a centralized Beidou pseudo-satellite indoor dynamic positioning method based on a known initial point of the present invention.

[0106] Step 1: The receiver collects pseudo-satellite observation data from the initial point and calculates the position of each pseudo-satellite transmitting antenna. Specifically, the UM980 receiver module of Unicore Starlink is used to collect Beidou pseudo-satellite raw observation data and ephemeris data from the initial point, store the carrier phase observation value and C / N0 of each satellite by epoch, and calculate the transmitting antenna position of each satellite through the ephemeris file.

[0107] Step 2: Using a known initial point and taking one satellite as a reference, calibrate the differential carrier phase integer ambiguity of other satellites and the relative satellite channel delay combination. The detailed steps are as follows:

[0108] (1) Observe n (n>=5) satellites at a known initial point, list the carrier phase observation equations of n satellites and select one of them as the reference satellite. The carrier phase equations of n satellites are as follows:

[0109]

[0110] Where L represents the carrier phase observation value, with the satellite number as the superscript and the epoch number as the subscript, N represents the integer ambiguity, with the satellite number as the superscript, r is the geometric distance between the satellite and the receiver, with the satellite number as the superscript and the epoch number as the subscript, δt r and δt s are the receiver clock error and satellite clock error respectively, ε is other random error terms, λ and C are the carrier wavelength and speed of light respectively.

[0111] (2) Considering that this centralized pseudo-satellite uses the same source clock, the error generated by the clock source in the satellite clock error is the same, so the observation equation is re-expressed as follows:

[0112]

[0113] Among them, t i is the clock error generated by the clock source of the i-th epoch, T is the satellite channel delay, and the first satellite is taken as the reference satellite. Since the other random error terms of the carrier phase observation are very small, their influence is ignored here. The carrier phase observation equations of the remaining n-1 satellites are subtracted from the carrier phase observation equations of the reference satellite, and the differential carrier phase integer ambiguity of the n-1 satellites and the relative satellite channel delay combination can be obtained:

[0114]

[0115] Step 3: After storing the differential integer ambiguity and the relative satellite channel delay combination, bring it into the satellite carrier phase observation equation corresponding to the subsequent epoch to perform data calibration. The detailed steps are as follows:

[0116] Assume that we calibrate the combination of differential integer ambiguity and relative satellite channel delay through the first epoch:

[0117]

[0118] (T 1 -T n )+λ(N 1 -N n) is brought into the subsequent epoch carrier phase observation equation, taking the second epoch as an example, we have:

[0119]

[0120] Step 4: Calculate the altitude angle of each pseudo-satellite based on the latest receiver position, which is then used to calculate the satellite position using the combined weighted least squares algorithm. The detailed steps are as follows:

[0121] (1) If the pseudo-satellite system uses an indoor relative coordinate system, assuming that the relative coordinates of the satellite are (x s ,y s ,z s )The relative coordinates of the receiver are (x r ,y r ,z r ), the relative coordinate system altitude angle calculation method is used.

[0122] Calculate the satellite's position vector relative to the receiver:

[0123] Δx=x s -x r

[0124] Δy=y s -y r

[0125] Δz=z s -z r

[0126] Calculate the horizontal distance d xy , which is the projection length of the relative position on the horizontal plane:

[0127]

[0128] Calculate the altitude angle:

[0129]

[0130] (2) If the pseudo-satellite system uses an absolute coordinate system, assuming that the ECEF coordinates of the satellite are (x s ,y s ,z s )The ECEF coordinates of the receiver are (x r ,y r ,z r ), the geodetic coordinates are (φ, λ, h), and the absolute coordinate system altitude angle calculation method is used.

[0131] The observation vector from the user to the satellite is:

[0132]

[0133] The unit observation vector I of the satellite at the receiver coordinate point is:

[0134]

[0135] Observation vector [Δx Δy Δz] T It can be equivalently expressed as the vector [Δe Δn Δu] in the station center coordinate system with the receiver coordinate point as the origin T , and its transformation relationship is:

[0136]

[0137] Among them, the coordinate transformation matrix S is:

[0138]

[0139] With the satellite observation vector at the receiver [Δe Δn Δu] T , the satellite's elevation angle θ relative to the receiver is:

[0140]

[0141] Step 5: Use a weighted least square method based on the median value of the mutual difference, satellite altitude angle, and C / N0 to solve the position and obtain high-precision receiver coordinates. Compare the positioning errors of different weighted combinations. The detailed steps are as follows:

[0142] (1) The observation equation after data adjustment is:

[0143]

[0144] Use the adjusted observation term on the left side of the equal sign in the above formula Indicated by TN, n is the satellite number, i is the epoch number, and the items on the right except the geometric distance between the satellite and the receiver are used as residual items. i Expressed as, and the geometric distance term is expanded, that is:

[0145]

[0146] in represents the receiver coordinates at the i-th epoch, Indicates the coordinates of the nth satellite, TN i is the residual term of the i-th epoch.

[0147] (2) In the first epoch, the differential carrier phase integer ambiguity and the relative satellite channel delay are combined. The receiver coordinates of the first epoch are the known initial coordinates. Starting from the second epoch, the combined weighted least squares algorithm is used to solve the receiver coordinates and residual terms. The receiver initial coordinates (x0, y0, z0) are set to the receiver coordinates of the previous epoch, and the initial residual term TN0 is set to 0. k represents the number of iterations currently being performed, and k-1 represents the number of iterations currently completed. When k = 1, the initial coordinates and TN0 are used as the initial values ​​of the iteration, otherwise the results of k-1 iterations are used as the initial values ​​of the kth iteration. For the kth iteration, the observation equation is (x k-1 ,y k-1 ,z k-1 ), TN k-1 Linearize at . Take the nth observation equation as an example to linearize it, ignoring the epoch number, and find the nonlinear term r n x r The partial derivative of , we get:

[0148]

[0149] Similarly, we can find r n For r , z r The partial derivative of , we get:

[0150]

[0151] Note that the right side of the observation equation is already a linear function of TN, or the partial derivative of the left side of the equation with respect to TN is 1. The observation equation group is in (x k-1 ,y k-1 ,z k-1 ), TN k-1 The result of linearization is:

[0152]

[0153] To simplify the equation, introduce matrices G and b:

[0154]

[0155]

[0156] The linearized equation can be simplified to:

[0157]

[0158] Considering the minimum sum of squares of the differences on both sides of the equation, we can get the least squares solution of the matrix equation:

[0159]

[0160] Update the receiver's position coordinates and residual items according to the above formula:

[0161]

[0162] The iteration result is compared with a determined threshold. If the result is greater than the threshold, the iteration continues. Otherwise, the iteration can be terminated and the solution of the current nonlinear equation group is used as the position coordinates of the receiver obtained by positioning and the residual term.

[0163] The idea of ​​weighted least squares is to give different weights to the observations of each satellite. The observations with larger weights play a greater role in solving the equation. The weighted model can reduce the weight of satellites with poor observation quality in positioning. The weighted model can also suppress the multipath effect to a certain extent and provide better positioning solution results. The steps of weighted least squares positioning solution are basically the same as those of least squares positioning solution. The difference is mainly reflected in the step of solving the matrix equation. The weighted least squares method introduces the weight matrix W to solve the matrix equation. Since there is no correlation between the measurements, the weight matrix is ​​usually a diagonal matrix, as follows:

[0164]

[0165] After setting the weights corresponding to the satellite measurements, multiply each item in the simplified linearized equation by the corresponding weight coefficient, and the corresponding matrix equation is:

[0166]

[0167] The least squares method is used to solve the matrix equation and the solution of the equation is obtained:

[0168]

[0169] (3) Considering the design of the weight matrix, a combined weighted least squares method based on the median value of the mutual difference, satellite altitude angle, and C / N0 is proposed, considering the following formula:

[0170]

[0171] It is in (x k-1 ,y k-1 ,z k-1 ), TN k-1 The residual term at is the deviation between the pseudorange value corresponding to the actual position of the receiver and the corresponding pseudorange value at the linearization point. Other errors that are ignored are also included in TN and b i TN k-1The proportion of b in the whole item is much greater than that of other items, so i It cannot accurately reflect the size of the error and deviation. i Arrange in ascending order to get b asc ,Right now Let mid = fix((n+1) / 2), where fix is ​​the rounding function, and take b asc The middle value of Subtract b from each element in the original array b mid And take the absolute value to get the new array a=[a1, a2, ..., a n ]

[0172] in:

[0173] For the same linearization point, after calculating the mutual difference, the common errors such as the residual term and receiver noise error shared by different satellites at the same epoch are eliminated. In practice, due to the influence of various errors, b is positive or negative, so we take the b in the middle of the b array mid The weight of the corresponding satellite is the largest, set it to c mid =(0.55-n / 20), which changes according to the number n of satellites involved in positioning.

[0174] The remaining weight coefficients are:

[0175] The weight coefficient matrix ω1 can be constructed:

[0176]

[0177] Next, consider the weight coefficient matrix ω2 based on the altitude angle. Satellites with high altitude angles often have smaller multipath errors. The weight coefficient matrix is ​​constructed using the sine function model:

[0178]

[0179] in:

[0180]

[0181] where α i = sin(θ i ), θ i is the satellite altitude angle calculated previously.

[0182] Then construct the weight coefficient matrix ω3 based on C / N0. Often, satellites with high C / N0 have better signal quality and smaller errors. According to experiments, when the transmission power and gain are fixed to the set value, the size of C / N0 is between 40dbHZ and 60dbHZ. When C / N0 is less than 40dbHZ, the satellite is considered unavailable. Therefore, we can make:

[0183]

[0184] The weight matrix is ​​designed as:

[0185] in:

[0186]

[0187] Finally, the combination weights are designed by combining the first three weighting methods, and the combination weight coefficient is designed as follows:

[0188]

[0189] The combined weight matrix is:

[0190]

[0191] High-precision positioning can be achieved by taking ω as the weight matrix and bringing it into weighted least squares.

[0192] (4) Design a dynamic experiment and compare and analyze different weighted combinations. The receiver moves on the trajectory and stays at known point 1 and known point 2 for 60 epochs. The positioning errors in all directions of the coordinates of the two points calculated using the least squares method with different weighted combinations are compared and analyzed, and the average error and standard deviation are calculated.

[0193] The present invention provides a centralized Beidou pseudo-satellite indoor dynamic positioning method based on a known initial point. Based on the known initial point, by determining the reference satellite, the differential carrier phase integer ambiguity of other satellites and the relative satellite channel delay combination are calibrated by inter-satellite difference, the differential integer ambiguity and the relative satellite channel delay combination are stored, and then brought into the satellite observation equation corresponding to the subsequent epoch to perform data correction, and finally a combined weighted least square method is used to perform position solution to obtain high-precision receiver coordinates. This method can effectively solve the problem of large indoor multipath error, poor pseudo-range observation value accuracy, and the inability to fix the integer ambiguity.

[0194] The above description is only a specific implementation of the present invention, and the protection scope of the present invention is not limited thereto. Any person familiar with the technology can understand and think of any changes or substitutions within the technical scope disclosed by the present invention, which should be included in the scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims

1. A centralized Beidou pseudo-satellite indoor dynamic positioning method based on a known initial point, characterized in that: The steps include: Step 1: The receiver starts to collect pseudolite observation data from the initial point and calculates the position of each pseudolite transmitting antenna; Step 2: Using a known initial point and taking one satellite as a reference, calibrate the differential carrier phase integer ambiguity and relative satellite channel delay combination of other satellites; Step 3: After storing the differential integer ambiguity and the relative satellite channel delay combination, bring it into the satellite carrier phase observation equation corresponding to the subsequent epoch to perform data adjustment; Step 4: Calculate the altitude angle of each pseudo-satellite based on the latest receiver position; Step 5: Use a weighted least square method based on the median value of the mutual difference, satellite altitude angle, and C / N0 to solve the position and obtain high-precision receiver coordinates. Finally, design a dynamic experiment to analyze and compare the different weighted combinations.

2. The method for indoor dynamic positioning of a centralized Beidou pseudo-satellite based on a known initial point according to claim 1, characterized in that: In the step 1, the receiver collects pseudo-satellite observation data from the initial point and solves the position of each pseudo-satellite transmitting antenna, specifically: the UM980 receiver module of Unicore Starlink is used to collect Beidou pseudo-satellite original observation data and ephemeris data from the initial point, the carrier phase observation value and C / N0 of each satellite are stored according to the epoch, and the transmitting antenna position of each satellite is calculated through the ephemeris file.

3. The method for indoor dynamic positioning of a centralized Beidou pseudo-satellite based on a known initial point according to claim 1, characterized in that: In step 2, the differential carrier phase integer ambiguity and relative satellite channel delay combination of other satellites are calibrated with one satellite as the reference by using the known initial point, specifically: (1) Observe n (n>=5) satellites at a known initial point, list the carrier phase observation equations of n satellites, and select one of the satellites as the reference satellite. The carrier phase equations of n satellites are as follows: Where L represents the carrier phase observation value, with the satellite number as the superscript and the epoch number as the subscript, N represents the integer ambiguity, with the satellite number as the superscript, r is the geometric distance between the satellite and the receiver, with the satellite number as the superscript and the epoch number as the subscript, δt r and δt s are the receiver clock error and satellite clock error respectively, ε is other random error terms, λ and C are the carrier wavelength and speed of light respectively. (2) Considering that this centralized pseudo-satellite uses the same source clock, the error generated by the clock source in the satellite clock error is the same, so the observation equation is re-expressed as follows: Among them, t i is the clock error generated by the clock source of the i-th epoch, T is the satellite channel delay, and the first satellite is taken as the reference satellite. Since the other random error terms of the carrier phase observation are very small, their influence is ignored here. The carrier phase observation equations of the remaining n-1 satellites are subtracted from the carrier phase observation equations of the reference satellite, and the differential carrier phase integer ambiguity of the n-1 satellites and the relative satellite channel delay combination can be obtained:

4. The method for indoor dynamic positioning of a centralized Beidou pseudo-satellite based on a known initial point according to claim 1, characterized in that: In step 3, after storing the differential integer ambiguity and the relative satellite channel delay combination, it is brought into the satellite observation equation corresponding to the subsequent epoch to perform data correction, which is specifically: Assume that we calibrate the combination of differential integer ambiguity and relative satellite channel delay through the first epoch: (T 1 -T n )+λ(N 1 -N n ) is brought into the subsequent epoch carrier phase observation equation, taking the second epoch as an example, we have:

5. The method for indoor dynamic positioning of centralized Beidou pseudo-satellite based on a known initial point according to claim 1, characterized in that: In step 4, the altitude angle of each pseudo-satellite is calculated with the latest receiver position, which is used for the subsequent combined weighted least squares algorithm to solve the satellite position, which is specifically as follows: (1) If the pseudo-satellite system uses an indoor relative coordinate system, assuming that the relative coordinates of the satellite are (x s ,y s ,z s )The relative coordinates of the receiver are (x r ,y r ,z r ), the relative coordinate system altitude angle calculation method is used. Calculate the satellite's position vector relative to the receiver: Δx=x s -x r Δy=y s -y r Δz=z s -with r Calculate the horizontal distance d xy , which is the projection length of the relative position on the horizontal plane: Calculate the altitude angle: (2) If the pseudo-satellite system uses an absolute coordinate system, assuming that the ECEF coordinates of the satellite are (x s ,y s ,z s ) The ECEF coordinates of the receiver are (x r ,y r ,z r ), the geodetic coordinates are (φ, λ, h), and the absolute coordinate system altitude angle calculation method is used. The observation vector from the user to the satellite is: The unit observation vector I of the satellite at the receiver coordinate point is: Observation vector [Δx Δy Δz] T It can be equivalently expressed as the vector [Δe Δn Δu] in the station center coordinate system with the receiver coordinate point as the origin T , and its transformation relationship is: Among them, the coordinate transformation matrix S is: With the satellite observation vector at the receiver [Δe Δn Δu] T , the satellite's elevation angle θ relative to the receiver is:

6. The method for indoor dynamic positioning of centralized Beidou pseudo-satellite based on a known initial point according to claim 1, characterized in that: In step 5, a weighted least square method based on the median value of the mutual difference, the satellite altitude angle and the carrier-to-noise ratio is used to solve the position and obtain the high-precision receiver coordinates. Finally, a dynamic experiment is designed to analyze and compare the different weighted combinations. Specifically: (1) The observation equation after data adjustment is: Use the adjusted observation term on the left side of the equal sign in the above formula Indicated by TN, n is the satellite number, i is the epoch number, and the items on the right except the geometric distance between the satellite and the receiver are used as residual items. i Expressed as, and the geometric distance term is expanded, that is: in represents the receiver coordinates at the i-th epoch, Indicates the coordinates of the nth satellite, TN i is the residual term of the i-th epoch. (2) In the first epoch, the differential carrier phase integer ambiguity and the relative satellite channel delay are combined. The receiver coordinates of the first epoch are the known initial coordinates. Starting from the second epoch, the combined weighted least squares algorithm is used to solve the receiver coordinates and residual terms. The receiver initial coordinates (x0, y0, z0) are set to the receiver coordinates of the previous epoch, and the initial residual term TN0 is set to 0. k represents the number of iterations currently being performed, and k-1 represents the number of iterations currently completed. When k = 1, the initial coordinates and TN0 are used as the initial values ​​of the iteration, otherwise the results of k-1 iterations are used as the initial values ​​of the kth iteration. For the kth iteration, the observation equation is (x k-1 ,y k-1 ,z k-1 ), TN k-1 Linearize at . Take the nth observation equation as an example to linearize it, ignoring the epoch number, and find the nonlinear term r n x r The partial derivative of , we get: Similarly, we can find r n For r , z r The partial derivative of , we get: Note that the right side of the observation equation is already a linear function of TN, or the partial derivative of the left side of the equation with respect to TN is 1. The observation equation group is in (x k-1 ,y k-1 ,z k-1 ), TN k-1 The result of linearization is: To simplify the equation, introduce matrices G and b: The linearized equation can be simplified to: Considering the minimum sum of squares of the differences on both sides of the equation, we can get the least squares solution of the matrix equation: Update the receiver's position coordinates and residual items according to the above formula: The iteration result is compared with a determined threshold. If the result is greater than the threshold, the iteration continues. Otherwise, the iteration can be terminated and the solution of the current nonlinear equation group is used as the position coordinates of the receiver obtained by positioning and the residual term. The idea of ​​weighted least squares is to give different weights to the observations of each satellite. The observations with larger weights play a greater role in solving the equation. The weighted model can reduce the weight of satellites with poor observation quality in positioning. The weighted model can also suppress the multipath effect to a certain extent and provide better positioning solution results. The steps of weighted least squares positioning solution are basically the same as those of least squares positioning solution. The difference is mainly reflected in the step of solving the matrix equation. The weighted least squares method introduces the weight matrix W to solve the matrix equation. Since there is no correlation between the measurements, the weight matrix is ​​usually a diagonal matrix, as follows: After setting the weights corresponding to the satellite measurements, multiply each item in the simplified linearized equation by the corresponding weight coefficient, and the corresponding matrix equation is: The least squares method is used to solve the matrix equation and the solution of the equation is obtained: (3) Considering the design of the weight matrix, a combined weighted least squares method based on the median value of the mutual difference, satellite altitude angle, and C / N0 is proposed, considering the following formula: It is in (x k-1 ,y k-1 ,z k-1 ), TN k-1 The residual term at is the deviation between the pseudorange value corresponding to the actual position of the receiver and the corresponding pseudorange value at the linearization point. Other errors that are ignored are also included in TN and b i TN k-1 The proportion of b in the whole item is much greater than that of other items, so i It cannot accurately reflect the size of the error and deviation. i Arrange in ascending order to get b asc ,Right now Let mid = fix((n+1) / 2), where fix is ​​the rounding function, and take b asc The middle value of Subtract b from each element in the original array b mid And take the absolute value to get the new array a=[a1, a2, ..., a n ] in: For the same linearization point, after calculating the mutual difference, the common errors such as the residual term and receiver noise error shared by different satellites at the same epoch are eliminated. In practice, due to the influence of various errors, b is positive or negative, so we take the b in the middle of the b array mid The corresponding satellite has the largest weight, which is set to c mid =(0.55-n / 20), which changes according to the number n of satellites involved in positioning. The remaining weight coefficients are: The weight coefficient matrix ω1 can be constructed: Next, consider the weight coefficient matrix ω2 based on the altitude angle. Satellites with high altitude angles often have smaller multipath errors. The weight coefficient matrix is ​​constructed using the sine function model: in: where α i = sin(θ i ), θ i is the satellite altitude angle calculated previously. Then construct the weight coefficient matrix ω3 based on C / N0. Often, satellites with high C / N0 have better signal quality and smaller errors. According to experiments, when the transmission power and gain are fixed to the set value, the size of C / N0 is between 40dbHZ and 60dbHZ. When C / N0 is less than 40dbHZ, the satellite is considered unavailable. Therefore, we can make: The weight matrix is ​​designed as: in: Finally, the combination weights are designed by combining the first three weighting methods, and the combination weight coefficient is designed as follows: The combined weight matrix is: High-precision positioning can be achieved by taking ω as the weight matrix and bringing it into weighted least squares. (4) Design a dynamic experiment and compare and analyze different weighted combinations. The receiver moves on the trajectory and stays at known point 1 and known point 2 for 60 epochs. The positioning errors in all directions of the coordinates of the two points calculated using the least squares method with different weighted combinations are compared and analyzed, and the average error and standard deviation are calculated.

Images