-
摘要: 雪水当量的监测对于气候变化的预测、水资源管理、农业生产规划具有重要意义。GPS干涉反射(GPS interferometric reflectometry, GPS-IR)技术是一种十分有效的地表积雪监测技术，基于GPS-IR技术提出了一种雪水当量的快速估计方法。首先基于GPS-IR技术获取美国板块边界观测（plate boundary observatory，PBO）GPS站的雪深时间序列；然后利用美国积雪遥测（SNowTELemetry, SNOTEL）站观测数据构建雪水当量转换模型；最后以北美历史与预测气候数据项目（historical and projected climate data for North America，ClimateNA）的气候预测数据作为参数约束，将GPS日雪深快速转化为雪水当量，并对雪水当量估计与验证过程的影响因素进行评价。实验结果表明，基于GPS-IR技术得到的雪深序列具有良好可靠性，与观测值的相关系数（R2）达到0.98，均方根误差（root mean square error, RMSE）为11.1 cm，偏差（Bias）为-3.7 cm；快速转化模型对雪水当量估计具有较高精度（R2=0.98，RMSE=4.2 cm，Bias=-2.5 cm）与稳定性；转化模型时空稳定性较高，残差集中在5 cm内；气候预测数据的引入、积雪分布差异对雪水当量估计与验证影响较小。所提方法在积雪监测设备缺乏区域可实现雪水当量快速估计，同时也为现有积雪观测网络增强、积雪产品改善等研究提供参考。Abstract:Objectives Snow water equivalent (SWE) plays a key role in climate change prediction, water resource management and agricultural production planning. GPS interferometric reflectometry (GPS-IR) has been proven to be a powerful tool to monitor snow depth. An efficient framework is presented to rapidly estimate SWE from GPS-IR derived snow depth.Methods Firstly, the daily snow depth product is obtained using GPS observations with GPS-IR technique. Secondly, an SWE conversion model is constructed using snow depth, SWE and climate observations from snow telemetry (SNOTEL) stations in the study region. Finally, using the climate forecast data provided by the historical and projected climate data for North America (ClimateNA) project as parameter constraints, the daily GPS snow depth product is converted into SWEs.Results The application of the proposed framework to GPS data from the plate boundary observatory (PBO), USA shows that GPS-IR derived snow depth product is reliable (R2 = 0.98, RMSE(root mean square error)= 11.1 cm, Bias =-3.7 cm).The daily GPS snow depth can be converted into SWE with high reliability (R2 = 0.98, RMSE = 4.2 cm, Bias =-2.5 cm) and stability. The rapid conversion model has high spatiotemporal stability with most of the residuals lying with the range of [-5, 5] cm. The additional uncertainties introduced by the climate forecast data and the spatial variations of snow depths have limited impacts on the estimation of SWE.Conclusions It is believed that the proposed framework can not only provide guidance to rapidly estimate SWE inregions lacking snow monitoring equipment, but also provide a reference to enhance existing snow observation network and improve accumulated snow products.
-
近些年无人机和自动驾驶等技术的快速发展引起了社会的极大关注，这些技术对实时位置的精度要求也越来越高。尤其是对于无人机机群和自动驾驶车队而言，载体精确、稳定的相对位置信息是保持队形的关键[1]。全球导航卫星系统（global navigation satellite system，GNSS）作为一种全球全天候的无线电定位技术，已广泛应用于各类位置信息获取的场景中[2-3]。如何利用GNSS数据获取移动载体精确、稳定的实时动态位置信息，是目前卫星导航领域的研究热点之一。
相较于经典最小二乘解算方法，卡尔曼滤波具有计算速度快、抗干扰能力强、相邻历元间位置信息相关联等优势[4-5]。将卡尔曼滤波应用于动态单点定位算法，可获得更加平滑和精确的定位结果，尤其是在观测卫星数较少的情况下定位精度改善明显[6]。在实时动态定位过程中，滤波模型的选择对定位结果具有较大影响。联合伪距与多普勒观测值，文献[7]利用常速模型来近似描述载体的运动状态，实现了米级的实时动态定位，但该模型的精度限制了滤波的最终效果。文献[8]利用载波相位时间差分得到的载体位置变化量来构建状态方程，对伪距观测值进行滤波更新，实现了对载体状态变化的准确刻画，但因只使用了伪距观测值进行绝对定位，最终定位精度仍有提升空间。等价消去原理[9-10]可通过对观测方程的等价变换实现对特定误差参数的消去。通过对求解问题无关的若干独立参数进行移除，可达到简化观测方程和提高计算效率的目的[11]。利用该原理并结合精密单点定位模型，文献[12]首先通过消去位置及钟差参数获得了仅含模糊度参数的等价观测方程，然后利用准静态滤波实时估计模糊度的浮点解，最终获得了移动载体的实时位置。由于采用了精密单点定位模型，该方法依赖于精密星历和钟差产品，且收敛速度较慢，实时性较差，不太适用于实时定位。
针对动态定位过程中载体运动状态不易准确描述的问题，同时考虑到在诸如车辆或无人机编队等应用场景中用户更关注载体间的运动轨迹是否平滑而对其绝对位置精度要求不高，本文提出了一种基于等价消去原理的实时动态定位方法，利用时间差分观测值和准静态卡尔曼滤波来实现移动载体的高精度轨迹确定。将精密单点定位模型中的观测量进行时间差分，则滤波方程可转化为对历元间模型残差的估计。通过随机游走对模型残差的随机特性进行描述，并利用滤波后的模型残差对载体的位置变化进行修正，仅采用广播星历就实现了滤波的快速收敛和载体位置变化的精确估计，极大削弱了模型残差随时间累积造成的滤波发散风险。
本文首先简要介绍了基于伪距及多普勒观测值的常速滤波模型、基于相位时间差分的位置变化量滤波模型，并提出了基于等价消去原理的单点准静态滤波模型；然后分别利用上述3种滤波模型，对无人机平台和地面移动平台采集的GPS和BDS观测数据进行动态定位数据处理，并对3种动态滤波模型的实验结果进行对比分析。
1 滤波模型
离散线性系统的卡尔曼滤波主要包括状态方程和观测方程两个部分，计算公式分别为[13]：
$$ {\mathit{\boldsymbol{X}}}_{k}={\mathit{\boldsymbol{\phi }}}_{k, k-1}{\mathit{\boldsymbol{X}}}_{k-1}+{\mathit{\boldsymbol{U}}}_{k-1}+{\mathit{\boldsymbol{W}}}_{k} $$ (1) $$ {\mathit{\boldsymbol{Z}}}_{k}={\mathit{\boldsymbol{H}}}_{k}{\mathit{\boldsymbol{X}}}_{k}+{\mathit{\boldsymbol{\omega }}}_{k} $$ (2) 式中，下角标$ k $和$ k-1 $表示相邻两个历元；$ \mathit{\boldsymbol{X}} $和$ \mathit{\boldsymbol{\phi }} $分别为移动载体的状态向量和状态转移矩阵；$ \mathit{\boldsymbol{U}} $为控制输入向量；$ \mathit{\boldsymbol{W}} $为过程噪声；$ \mathit{\boldsymbol{Z}} $为观测值向量；$ \mathit{\boldsymbol{H}} $为观测方程的系数矩阵；$ \mathit{\boldsymbol{\omega }} $为测量噪声。为方便表述，后文不再区分下角标。
为简化滤波模型，本文的观测方程中卫星轨道及卫星钟差直接采用广播星历进行确定，电离层延迟误差采用Klobuchar模型改正[14-16]，对流层延迟误差采用Hopfield模型改正[17]，其余误差包括硬件延迟偏差、潮汐改正等均被纳入非模型化残余误差。此外，对于观测值的随机模型，假设各卫星的伪距及载波相位观测值均相互独立且不考虑其时间相关性，其测量噪声的方差σ2利用高度角模型定义为[18]：
$$ {\sigma }^{2}=\left\{\begin{array}{ll}{\sigma }_{0}^{2}/\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}E, 0° < E < 30°& \\ {\sigma }_{0}^{2}, E\ge 30°& \end{array}\right. $$ (3) 式中，$ E $为卫星高度角；$ {\sigma }_{0} $为观测值的标称精度，对于伪距和载波相位观测值，按经验分别取3 m和2 mm[19]。差分观测值的测量噪声可在此基础上通过误差传播定律确定。多普勒观测值的噪声方差按文献[7]的推荐取固定值0.18（周/s）2。
1.1 基于伪距及多普勒观测值的常速滤波模型
假设移动载体在运动时速度较慢且观测数据采样率较高，其运动状态则可通过常速模型进行描述，此时卡尔曼滤波模型的状态方程式（1）中状态转移矩阵的表达式固定不变且输入项$ \mathit{\boldsymbol{U}} $的元素均为0。具体到伪距及多普勒观测方程，根据移动载体的运动学方程，此时状态向量$ \mathit{\boldsymbol{X}} $及其转移矩阵$ \mathit{\boldsymbol{\phi }} $分别表示为：
$$ \mathit{\boldsymbol{X}}={\left[\begin{array}{cccc}\mathit{\boldsymbol{x}}& {t}_{r}& \mathit{\boldsymbol{v}}& {\dot{t}}_{r}\end{array}\right]}^{\mathrm{T}} $$ (4) $$ \mathit{\boldsymbol{\phi }}=\left[\begin{array}{cccc}\mathit{\boldsymbol{I}}& 0& \mathrm{\Delta }t\bullet \mathit{\boldsymbol{I}}& 0\\ 0& 1& 0& \mathrm{\Delta }t\\ 0& 0& \mathit{\boldsymbol{I}}& 0\\ 0& 0& 0& \mathit{\boldsymbol{I}}\end{array}\right] $$ (5) 式中，$ \mathit{\boldsymbol{x}} $和$ \mathit{\boldsymbol{v}} $分别为移动载体的位置向量和速度向量；$ {t}_{r} $和$ {\dot{t}}_{r} $分别表示以距离为单位的接收机钟差和以速度为单位的接收机钟速；$ \mathit{\boldsymbol{I}} $表示单位阵，维数为3×3；$ \mathrm{\Delta }t $为采样间隔，在动态定位中一般设为1 s。上述状态方程的过程噪声方差阵$ {\mathit{\boldsymbol{D}}}_{W} $按文献[7]进行取值。
在忽略硬件延迟影响的情况下，联合线性化后的伪距及多普勒观测方程，即可得到如式（2）形式的滤波观测方程。方程中观测值向量$ \mathit{\boldsymbol{Z}} $和系数矩阵$ \mathit{\boldsymbol{H}} $可具体表示为：
$$ \mathit{\boldsymbol{Z}}={\left[\begin{array}{cc}P& -\lambda D\end{array}\right]}^{\mathrm{T}} $$ (6) $$ \mathit{\boldsymbol{H}}=\left[\begin{array}{cc}\mathit{\boldsymbol{B}}& 0\\ 0& \mathit{\boldsymbol{B}}\end{array}\right] $$ (7) 式中，$ P $为伪距观测值；$ D $为多普勒观测值；$ \lambda $为多普勒观测值对应的波长；$ \mathit{\boldsymbol{B}}=\left[l\ m\ n\ 1\right] $为伪距及多普勒观测方程的系数矩阵，其中$ l、m $和$ n $为测站到卫星的方向余弦。
1.2 基于相位时间差分的位置变化量滤波模型
载波相位观测值比伪距具有更高的观测精度，因此常用在高精度定位场景中。在无周跳发生的情况下，载波相位观测方程中的整周模糊度参数可通过历元间差分予以消除。同时，大气延迟误差在较短时间间隔的相邻历元中可视为不变，该类误差亦可通过历元间差分予以消除。因此，当相邻历元的方向余弦向量视为不变时，线性化后的载波相位时间差分观测方程可表示为[20-21]：
$$ \lambda \Delta \varphi =\mathit{\boldsymbol{B}}{\left[\begin{array}{cc}\Delta \mathit{\boldsymbol{x}}& \Delta {t}_{r}\end{array}\right]}^{\mathrm{T}}+\Delta {\rho }^{0}-\Delta {t}^{s}+{\epsilon }_{\Delta \varphi } $$ (8) 式中，$ \Delta $为时间差分算子；$ \varphi $为载波相位观测值；$ {\rho }^{0} $为地卫距的近似值；$ {\epsilon }_{\Delta \varphi } $为非模型化误差；此时方程中的参数$ \Delta \mathit{\boldsymbol{x}} $即为移动载体在相邻历元的位置变化量；$ \Delta {t}_{r} $为接收机的钟差变化量；$ \Delta {t}^{s} $为卫星的钟差变化量。
对式（8）进行最小二乘平差解算，即可获得$ \mathrm{\Delta }\mathit{\boldsymbol{x}} $、$ \mathrm{\Delta }{t}_{r} $的估值及其方差信息。此时，再联合伪距和载波相位观测方程，即可构建基于相位时间差分的位置变化量滤波模型。模型中滤波状态方程的状态转移矩阵$ \mathit{\boldsymbol{\phi }} $变为单位阵，状态向量$ \mathit{\boldsymbol{X}} $及输入向量$ \mathit{\boldsymbol{U}} $分别表示为：
$$ \mathit{\boldsymbol{X}}={\left[\begin{array}{ccc}\mathit{\boldsymbol{x}}& {t}_{r}& \mathit{\boldsymbol{N}}\end{array}\right]}^{\mathrm{T}} $$ (9) $$ \mathit{\boldsymbol{U}}={\left[\begin{array}{ccc}\mathrm{\Delta }\mathit{\boldsymbol{x}}& \mathrm{\Delta }{t}_{r}& \mathrm{\Delta }\mathit{\boldsymbol{N}}\end{array}\right]}^{\mathrm{T}} $$ (10) 式中，$ \mathit{\boldsymbol{N}} $为模糊度向量，在无周跳情况下其状态变化$ \mathrm{\Delta }\mathit{\boldsymbol{N}} $视为0。过程噪声矩阵$ \mathit{\boldsymbol{W}} $中位置与钟差部分的噪声方差直接采用上述平差求得的验后方差，模糊度部分各参数均视为独立且其噪声设为0.005周。
在式（2）形式的伪距和载波相位滤波观测方程中，观测值向量$ \mathit{\boldsymbol{Z}} $和系数矩阵$ \mathit{\boldsymbol{H}} $可具体表示为：
$$ \mathit{\boldsymbol{Z}}={\left[\begin{array}{cc}P& \lambda \varphi \end{array}\right]}^{\mathrm{T}} $$ (11) $$ \mathit{\boldsymbol{H}}=\left[\begin{array}{cc}\mathit{\boldsymbol{B}}& 0\\ \mathit{\boldsymbol{B}}& \lambda \mathit{\boldsymbol{I}}\end{array}\right] $$ (12) 此时单位阵$ \mathit{\boldsymbol{I}} $的维数为$ n\times n $，$ n $为观测卫星的个数。
值得说明的是，根据文献[22]，在利用式（8）计算移动载体的位置与钟差变化量的同时，还可以根据获得的验后单位权中误差并基于抗差最小二乘获得的观测值残差，对相位观测值的周跳进行探测与修复。对于发生周跳的卫星，状态方程中$ \mathrm{\Delta }\mathit{\boldsymbol{N}} $需对应修改为所确定的周跳值。
1.3 基于等价消去原理的单点准静态滤波模型
根据等价消去原理[10, 12]，在伪距及载波相位联合观测方程中，未知参数包括测站位置和接收机钟差等时变参数$ {\mathit{\boldsymbol{X}}}_{1}={\left[\begin{array}{cc}\mathit{\boldsymbol{x}}& {t}_{r}\end{array}\right]}^{\mathrm{T}} $及不随时间变化的模糊度参数$ {\mathit{\boldsymbol{X}}}_{2}={\mathit{\boldsymbol{N}}}^{\mathrm{T}} $。将式（11）、式（12）代入式（2），可重新将观测方程表述为：
$$ \mathit{\boldsymbol{Z}}={\mathit{\boldsymbol{H}}}_{1}{\mathit{\boldsymbol{X}}}_{1}+{\mathit{\boldsymbol{H}}}_{2}{\mathit{\boldsymbol{X}}}_{2}+\mathit{\boldsymbol{\omega }} $$ (13) 式中，$ {\mathit{\boldsymbol{H}}}_{1}=\left[\begin{array}{c}\mathit{\boldsymbol{B}}\\ \mathit{\boldsymbol{B}}\end{array}\right];{\mathit{\boldsymbol{H}}}_{2}=\left[\begin{array}{c}0\\ \lambda \mathit{\boldsymbol{I}}\end{array}\right] $。设观测值的权阵为$ \mathit{\boldsymbol{P}} $，则式（13）的法方程可表示为：
$$ \ \left[\begin{array}{cc}{\mathit{\boldsymbol{H}}}_{1}^{\mathrm{T}}\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{H}}}_{1}& {\mathit{\boldsymbol{H}}}_{1}^{\mathrm{T}}\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{H}}}_{2}\\ {\mathit{\boldsymbol{H}}}_{2}^{\mathrm{T}}\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{H}}}_{1}& {\mathit{\boldsymbol{H}}}_{2}^{\mathrm{T}}\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{H}}}_{2}\end{array}\right]\left[\begin{array}{c}{\mathit{\boldsymbol{X}}}_{1}\\ {\mathit{\boldsymbol{X}}}_{2}\end{array}\right] = \\ \left[\begin{array}{cc}{\mathit{\boldsymbol{M}}}_{11}& {\mathit{\boldsymbol{M}}}_{12}\\ {\mathit{\boldsymbol{M}}}_{21}& {\mathit{\boldsymbol{M}}}_{22}\end{array}\right]\left[\begin{array}{c}{\mathit{\boldsymbol{X}}}_{1}\\ {\mathit{\boldsymbol{X}}}_{2}\end{array}\right]=\left[\begin{array}{c}{\mathit{\boldsymbol{H}}}_{1}^{T}\mathit{\boldsymbol{P}}\mathit{\boldsymbol{Z}}\\ {\mathit{\boldsymbol{H}}}_{2}^{T}\mathit{\boldsymbol{P}}\mathit{\boldsymbol{Z}}\end{array}\right] $$ (14) 式中，$ {\mathit{\boldsymbol{M}}}_{ij} $（i、j可分别取值为1和2）为过渡的系数矩阵。经过消参变换[12, 23]后，式（14）可拆分成两个新的法方程为：
$$ \left\{\begin{array}{l}{\mathit{\boldsymbol{M}}}_{11}{\mathit{\boldsymbol{X}}}_{1}+{\mathit{\boldsymbol{M}}}_{12}{\mathit{\boldsymbol{X}}}_{2}={\mathit{\boldsymbol{H}}}_{1}^{\mathrm{T}}\mathit{\boldsymbol{P}}\mathit{\boldsymbol{Z}}\\ {\stackrel{-}{\mathit{\boldsymbol{H}}}}_{2}^{\mathrm{T}}\mathit{\boldsymbol{P}}{\stackrel{-}{\mathit{\boldsymbol{H}}}}_{2}{\mathit{\boldsymbol{X}}}_{2}={\stackrel{-}{\mathit{\boldsymbol{H}}}}_{2}^{\mathrm{T}}\mathit{\boldsymbol{P}}\mathit{\boldsymbol{Z}}\end{array}\right. $$ (15) 式中，$ {\stackrel{-}{\mathit{\boldsymbol{H}}}}_{2}={\mathit{\boldsymbol{H}}}_{2}-{\mathit{\boldsymbol{H}}}_{1}{\mathit{\boldsymbol{M}}}_{11}^{-1}{\mathit{\boldsymbol{H}}}_{1}^{\mathrm{T}}\mathit{\boldsymbol{P}}{\mathit{\boldsymbol{H}}}_{2} $。
式（15）中第2个法方程对应一个新的观测方程$ \mathit{\boldsymbol{Z}}={\stackrel{-}{\mathit{\boldsymbol{H}}}}_{2}{\mathit{\boldsymbol{X}}}_{2} $，此时未知参数仅包含不随时间变化的模糊度参数。同时，法方程式（15）中观测值$ \mathit{\boldsymbol{Z}} $及其权阵$ \mathit{\boldsymbol{P}} $均与式（14）保持一致，说明消参后的观测方程与原观测方程等价。此时，可利用卡尔曼滤波对方程中不随时间变化的模糊度参数进行求解。
在本文数据处理中，由于仅采用广播星历和大气误差的经验模型改正，伪距和载波相位观测方程中均存在较大残余误差。为削弱残余误差的影响，本文采用历元间伪距及载波相位差分观测值$ \Delta \mathit{\boldsymbol{Z}} $，此时时变参数和时不变参数分别调整为测站位置和接收机钟差变化$ \Delta {\mathit{\boldsymbol{X}}}_{1} $及模糊度变化$ \Delta {\mathit{\boldsymbol{X}}}_{2} $。无周跳发生时$ \Delta {\mathit{\boldsymbol{X}}}_{2} $应为0，但考虑到观测方程中非模型化残余误差的影响，状态方程将采用准静态滤波的随机游走过程加以描述，即$ \mathit{\boldsymbol{\phi }} $为单位阵且无输入项$ \mathit{\boldsymbol{U}} $，此时系统状态向量为：
$$ \mathit{\boldsymbol{X}}=\Delta {\mathit{\boldsymbol{X}}}_{2} $$ (16) 其过程噪声$ \mathit{\boldsymbol{W}} $包含有各类非模型化误差及大气延迟的历元间残差，综合噪声设为0.02 m。对应的观测方程为：
$$ \Delta \mathit{\boldsymbol{Z}}={\stackrel{-}{\mathit{\boldsymbol{H}}}}_{2}\Delta {\mathit{\boldsymbol{X}}}_{2} $$ (17) 利用上述状态方程与观测方程组成卡尔曼滤波系统，即可对参数$ \Delta {\mathit{\boldsymbol{X}}}_{2} $进行滤波，此时该参数本质上已转化为模型残差。当获得模型残差$ \Delta {\mathit{\boldsymbol{X}}}_{2} $的滤波值后，即可逐历元计算接收机的位置增量及钟差增量，计算公式为：
$$ \Delta {\mathit{\boldsymbol{X}}}_{1}={\mathit{\boldsymbol{M}}}_{11}^{-1}\left({\mathit{\boldsymbol{H}}}_{1}^{\mathrm{T}}\mathit{\boldsymbol{P}}\Delta \mathit{\boldsymbol{Z}}-{\mathit{\boldsymbol{M}}}_{12}\Delta {\mathit{\boldsymbol{X}}}_{2}\right) $$ (18) 在获得接收机的位置增量后，即可从首历元开始逐历元递推出后续历元的接收机位置。
由式（18）可知，若不考虑模型残差的影响，即$ \Delta {\mathit{\boldsymbol{X}}}_{2}=0 $，该式则退化为单历元的时间差分伪距及相位的位置变化估计。因此，本文模型的优势在于采用随机游走过程对时间差分后的模型残差进行了描述与滤波，并利用残差的滤波值对接收机的坐标变化进行了精确修正。此外，由于随机游走过程的使用，此滤波模型削弱了模型残差随时间的累积，在一定程度上减轻了定位结果的发散。
2 不同滤波模型对比实验
实验采用的两组观测数据是2020-11-17分别由无人机平台和移动车辆平台在湘潭市郊区的开阔环境下采集得到的，两载体均搭载了上海司南的K708板卡接收机，无人机采用了华信的HX-CH7604A小型无人机天线，汽车采用了上海司南的AT360测量天线。两载体的运动状态均包括低速近似直线的往返移动及环绕移动，速度不超过10 m/s。
受无人机续航时间限制，实验中两载体每次动态数据采集时间持续约0.5 h，并重复数次数据采集。采集的数据包含GPS和北斗卫星导航系统（BeiDou satellite navigation system，BDS）双系统观测数据，采样频率为1 Hz，截止高度角为10°。由于本文的滤波模型方法均未考虑周跳发生的情况，因此为简化处理流程，在定位解算前先利用TurboEdit算法对观测数据进行周跳探测[24]，依据探测结果选取其中无周跳的10 min左右观测时段数据用于后续滤波模型的对比分析。此外，虽然实验重复采集了多组数据，但限于文章篇幅，本文仅选取其中一组进行结果展示，其他组数据的结果与本文结果类似。实验的解算过程利用课题组内部程序完成。
利用基于伪距及多普勒观测值的常速滤波模型（模型一）、基于相位时间差分的位置变化量滤波模型（模型二）、基于等价消去原理的单点准静态滤波模型（模型三）分别对GPS的L1频率观测值和BDS的B1I频率观测值进行单系统定位解算。考虑到模型一、模型二得到的是接收机的绝对坐标，而模型三得到的是接收机的位置变化，为便于性能评估，在滤波结束后计算各历元定位坐标相对于首历元单点定位结果在东（east，E）、北（north，N）和天（up，U）方向上的相对位置。同时，以实时动态定位解算得到的移动载体厘米级固定解为基础，也计算各历元相对首历元在E、N、U方向上的相对位置作为参考值。最后，通过对比滤波得到的相对初始历元的位置与其参考值间的偏差，分析3种滤波模型在动态定位收敛后的内符合定位精度，即在不考虑定位结果系统性偏差情况下定位轨迹与真实轨迹的一致性评估。直接对比相邻历元间的坐标差结果，避免了实际轨迹经一段时间后可能发散而无法在结果中反映出来的问题。
同时，还统计了不同滤波模型定位结果的收敛时间，由于所统计的定位结果偏差中可能包含有系统性偏差，因此该收敛时间被定义为E、N、U 3个方向上的坐标偏差相对其整体变化趋于稳定时的坐标偏差的差异值均小于0.2 m所花费的时间[25-26]。
2.1 无人机实验
第一组实验数据由飞行无人机平台采集，基于GPS和BDS的定位结果在E、N、U方向上的偏差以及观测卫星数随时间的变化分别如图 1和图 2所示。受伪距和多普勒观测值精度限制，模型一的定位结果偏差存在较大波动，而模型二和模型三的定位结果在100 s内均实现了收敛。因此，选取第100秒后所有历元定位结果的均值作为其收敛值。具体来看，模型二和模型三均在1 min内实现了定位结果的收敛，且模型三的收敛速度略快，这与等价消去后其状态方程能更精准地描述参数变化有关。
![]()
图 1 GPS的无人机定位结果偏差及观测卫星数Figure 1. Positioning Deviation of Unmanned Aerial Vehicle (UAV) and the Number of Observed Satellites in GPS![]()
图 2 BDS的无人机定位结果偏差及观测卫星数Figure 2. Positioning Deviation of the UAV and the Number of Observed Satellites in BDS模型一的定位结果存在明显的系统性偏差，模型二在各方向上的定位结果亦未收敛到零。由于两者获得的定位结果均为接收机的绝对位置，而图 1和图 2反映的均是以首历元位置为起点的相对位置的偏差，因此上述系统性偏差是首历元单点定位偏差与收敛后定位偏差的综合表现。对于模型三而言，由于获得的定位结果为接收机的历元间变化量，均可归算到与接收机首历元的相对位置，因此其定位结果偏差最后均基本收敛到零附近，仅在N方向存在相对较小的系统偏差。
值得指出的是，模型二在图 1中E方向和图 2中U方向的偏差曲线在收敛后仍存在较为明显的线性变化趋势，该变化趋势是由历元间差分后各项残差的逐渐累积引起的[12]。考虑到实验中无人机观测数据较短，在长时间观测定位时模型二可能会引起结果发散。相比较而言，模型三的定位结果偏差未观察到明显的线性变化，这是因为模型三顾及了时间差分后模型残差的影响，并利用这部分残差对坐标变化量进行了修正。
为进一步对比分析3种滤波模型定位收敛后的稳定性及不同系统间的性能差异，对各滤波模型的收敛时间、收敛后定位结果在E、N、U方向及平面（horizontal，H）、三维（three-dimensional，3D）的定位误差标准差（standard deviation，STD）（即扣除系统性偏差收敛值后的均方根误差）进行了统计，其结果如表 1所示。与收敛值的定义相对应，3种滤波模型定位结果偏差的STD亦只统计100 s以后的结果，且虽然模型一的结果无明显收敛过程，但其系统偏差仍定义为100 s后定位偏差的均值。
表 1 3种滤波模型的收敛时间及收敛后的定位偏差标准差(无人机实验)Table 1. Convergence Time and STD of Position Errors After Convergence for GPS and BDS (UAV Experiment)系统 滤波
模型收敛时间/s 标准差/m E N U H 3D GPS 模型一 0.188 0.218 0.299 0.288 0.415 模型二 43 0.115 0.071 0.072 0.135 0.153 模型三 23 0.028 0.038 0.107 0.046 0.117 BDS 模型一 0.143 0.145 0.163 0.204 0.261 模型二 33 0.037 0.064 0.053 0.074 0.091 模型三 10 0.026 0.013 0.044 0.029 0.053 由表 1可以发现，无论采用哪种滤波模型，整体而言BDS的收敛速度及STD均优于GPS。产生这种现象的原因是观测到的BDS卫星数明显多于GPS，从而带来更多的冗余观测值；并且BDS混合星座使得BDS卫星的图形结构在亚太地区更优，模型强度更强。同时，模型三的收敛时间显著低于模型二，在30 s内即实现了定位结果的收敛，GPS及BDS的收敛速度分别提升了46.5%和69.7%。
表 1的STD结果反映了各模型在100 s后定位结果的稳定性。通过3种滤波模型的对比结果可以发现，模型三在收敛后的定位稳定性明显优于其他两种方法，仅在GPS的高程方向不如模型二。模型三在BDS的平面STD优于3 cm水平，点位误差的STD也仅为5.3 cm。相比模型二方法而言，GPS及BDS的平面和三维STD分别提升了66.4%、23.5%和60.8%、41.8%，表现出优异的定位性能。
2.2 车辆实验
第二组实验数据由移动车辆平台采集，最终得到GPS和BDS基于3种滤波模型的定位结果在E/N/U方向上的偏差及观测卫星数随时间的变化如图 3和图 4所示。图 3、4中各模型所表现出的定位性能与图 1和图 2基本一致，相比其他两种模型，模型三几乎无明显收敛过程且其相对位置仍最接近参考值，仅在图 4的N方向第500—600历元时段内出现了一定的偏移，但很快又收敛到参考值。
![]()
图 3 GPS的汽车采样数据定位结果偏差及观测卫星数Figure 3. Positioning Deviation of Vehicle Sampling Data and the Number of Observed Satellites in GPS![]()
图 4 BDS的汽车采样数据定位结果偏差及观测卫星数Figure 4. Positioning Deviation of Vehicle Sampling Data and the Number of Observed Satellites in BDS对本组实验中3种滤波模型的收敛时间、完成收敛后相对定位结果在E、N、U各方向及H、3D的STD进行了统计，各指标的定义与无人机实验一致，其结果如表 2所示。由收敛时间的对比可以发现，无论采用哪种滤波方法，整体而言BDS的收敛速度及STD仍然显著优于GPS，由数十秒降至10 s以内。相比无人机实验的结果，收敛速度的提升主要得益于本组实验中测量型天线的使用。同时，模型三的收敛时间仍然低于模型二，对于GPS和BDS速度分别提升了71.7%和42.9%。
表 2 3种滤波模型的收敛时间及收敛后定位误差标准差(汽车实验)Table 2. Convergence Time and STD of Position Errors After Convergence for GPS and BDS (Vehicle Experiment)系统 滤波
模型收敛时间/s 标准差/m E N U H 3D GPS 模型一 0.198 0.328 0.399 0.383 0.553 模型二 46 0.170 0.041 0.064 0.175 0.186 模型三 13 0.036 0.031 0.119 0.048 0.128 BDS 模型一 0.155 0.242 0.344 0.287 0.448 模型二 7 0.090 0.127 0.033 0.156 0.159 模型三 4 0.051 0.065 0.057 0.083 0.100 由收敛后坐标偏差的STD统计结果可以发现，模型三在收敛后的定位稳定性仍然显著优于模型一。相比于模型二而言，GPS和BDS中模型三在高程方向不及模型二，但平面及三维STD均显著优于模型二。模型三的点位误差STD在10 cm左右，不及无人机实验的定位精度，主要是由于无人机数据采集过程中环境更为开阔，观测值不易受多路径等误差的干扰。
3 结果讨论
本文利用无人机和汽车采集的两组动态观测数据对3种滤波模型分别基于GPS和BDS进行了定位实验。基于滤波结果，对3种滤波模型在不同卫星系统下的收敛时间和收敛后的偏差标准差进行了统计。经过对比发现基于等价消去原理构建的模型具有较快的收敛速度，相对于模型二而言平均提升了50%左右。
另外，在定位精度的稳定性上，模型三在两组实验中均体现出了明显的优势，表现出了该方法优异的定位性能。综合3种滤波模型在BDS和GPS下的整体表现来看，无论采用哪种滤波方法，整体而言BDS的收敛速度及STD均优于GPS。
从数学模型的角度看，模型一虽然定位结果一般，但是在定位过程中无需考虑周跳的因素，且模型较为简单，容易实现；模型二存在残余误差累积导致定位结果发散的问题，但在求定位置变化量的同时可以进行周跳的探测与修复；模型三经过等价变换后参数得到了简化，且因顾及了时间差分后残余误差的影响，极大削弱了模型二中定位结果发散的问题，但在实际应用中仍需考虑周跳的影响。
从实际应用的角度看，本文提出的基于等价消去原理的实时动态单点定位滤波模型能够为载体提供较为准确的运动轨迹，进而可以为动态多目标系统提供相对稳定的动态基准。目标之间的相对位置一般很容易获取，一旦动态位置得到确定，那么整个系统给的每一个目标都将按照预计轨迹运行。因此，本文提出的动态定位模型对动态多目标的导航具有一定的应用价值。
4 结语
为了研究如何利用有限的观测数据获取更加准确和稳定的动态定位信息，本文在回顾传统的基于伪距及多普勒观测值的常速滤波模型和基于相位时间差分的位置变化量滤波模型后，提出了一种基于等价消去原理的单点准静态滤波模型，来实现实时的动态单点定位。
基于等价消去原理的单点准静态滤波模型根据等价消去原理移除了观测方程中随时间变化的位置和钟差参数，使观测方程和对应的滤波状态方程得到极大简化。然后选用历元间差分观测量构建卡尔曼滤波，使得滤波过程充分顾及模型残差的影响。利用该滤波得到的残差对位置信息进行了有效修正，极大削弱了直接采用相位时间差分获取位置变化量带来的模型残差的累积。利用无人机和汽车采集的动态观测数据对3种滤波模型分别进行了基于GPS和BDS的单频定位实验，依据定位结果对滤波模型的定位性能进行了对比分析，得出如下结论：
1）相比基于相位时间差分的位置变化量滤波模型，本文提出的基于等价消去原理的单点准静态滤波模型能实现更快的定位收敛时间和更稳定的厘米级轨迹确定，无人机实验和车辆实验中的收敛时间均优于30 s，相对首历元位置的3D定位偏差的标准差基本在10 cm以内，平面精度可达3 cm。
2）由于观测到的BDS卫星数明显多于GPS卫星数且混合星座能够构成更优的卫星图形结构，从而带来更多的冗余观测值和更高的模型强度，因此无论采用哪种滤波方法，整体而言BDS的收敛速度及稳定性均优于GPS。
就数学模型而言，模型三经过等价变换后得到了简化，但仍需要考虑周跳的影响。同时，模型三虽然实现了相对首历元位置的高精度单点定位，但无法避免首历元单点定位带来的绝对定位误差，因此只能实现运动轨迹形状的精准描述，即对真实运动轨迹的平移。在实际应用中，可利用其他途径精确确定首历元位置，从而保证后续历元真实运动轨迹的获取。
此外，本文实验所用的观测数据时长相对较短，各模型定位精度的长时稳定性仍需进一步分析。且本文实验仅采用了单频单系统的定位模式，未对多频多系统的情况进行测试。在模型二和模型三中，本文未考虑周跳的情况，在实际的应用中需要对其进行探测和修复。这些方面的问题也将是本课题下一步深入研究的方向。
致谢: 感谢UNAVCO提供的GPS观测数据以及SNOTEL网络提供的雪深、雪水当量参考资料。 -
![]()
图 2 PBO网络测站分布情况和P351站观测环境
Figure 2. Distribution of PBO GPS Stations and Observation Environment of P351 Station
![]()
图 7 P351站雪深与SNOTEL 490站雪深对比
Figure 7. Comparison Between Snow Depths of P351 and SNOTEL 490
![]()
图 9 结合SNOTEL气候变量的雪水当量估计值与观测值对比
Figure 9. Comparison Between the SWEs Predicted with SNOTEL Climate Data and the In-situ SWEs
![]()
图 10 结合SNOTEL气候变量的雪水当量估计值与原位观测值对比差值时间序列
Figure 10. Time-Series of the Differences Between the SWEs Predicted with SNOTEL Climate Data and the In-situ SWEs
![]()
图 11 结合ClimateNA气候变量的雪水当量估计值与观测值对比结果
Figure 11. Comparison Between the SWEs Predicted with ClimateNA Climate Data and the In-situ SWEs
![]()
图 12 结合ClimateNA气候变量的雪水当量估计值与原位观测值对比差值时间序列
Figure 12. Time-Series of the Differences Between the SWEs Predicted with ClimateNA Climate Data and the In-situ SWEs
![]()
图 13 结合P351站雪深与SNOTEL气候变量的雪水当量估计值与原位观测值对比差值时间序列
Figure 13. Time-Series of the Differences Between the SWEs Predicted with Snow Depth of P351 and SNOTEL Climate Data and In-situ SWEs
![]()
图 14 结合P351站雪深与ClimateNA气候变量的雪水当量估计值与原位观测值对比差值时间序列
Figure 14. Time-Series of the Differences Between the SWEs Predicted with Snow Depth of P351 and ClimateNA Climate Data and the In-situ SWEs
表 1 P351及SNOTEL站位置信息
Table 1 Information of P351 and SNOTEL Stations
测站 位置 海拔/m 与P351站高差/m 与P351站距离/km P351 114.719°W，43.874°N 2 692.61 450 114.674°W，43.603°N 2 676.14 -126.19 30.86 489 114.672°W，43.877°N 2 731.01 -415.75 3.81 490 114.714°W，43.875°N 2 566.42 -16.46 0.41 845 114.853°W，43.799°N 1 923.29 38.40 13.39 895 114.418°W，43.769°N 2 276.86 -769.32 27.24 
下载: 导出CSV
表 2 SNOTEL站气候变量观测值
Table 2 Climate Observed Variables of SNOTEL Stations
测站 PPTWT/mm TD/℃ 450 292 21.7 489 254 21.1 490 292 21.9 845 417 21.1 895 185 22.8
下载: 导出CSV
表 3 ClimateNA气候变量预测值
Table 3 Climate Predicted Variables of ClimateNA
测站 PPTWT/mm TD/℃ P351 287 24.4 450 321 24.9 489 243 25.0 490 283 24.5 845 417 24.9 895 219 26.6
下载: 导出CSV
表 4 结合SNOTEL气候变量的雪水当量估计值与原位观测值的对比
Table 4 Comparison Between the SWEs Predicted with SNOTEL Climate Data and the In-situ SWEs
测站 R2 RMSE/cm Bias/cm 450 0.99 5.9 -4.2 489 0.98 1.6 -0.4 490 0.99 2.0 -1.2 845 0.99 3.1 -0.4 895 0.98 1.0 -0.3
下载: 导出CSV
表 5 结合ClimateNA气候变量的雪水当量估计值与原位观测值对比结果
Table 5 Comparison Between the SWEs Predicted with ClimateNA Climate Data and the In-situ SWEs
测站 R2 RMSE/cm Bias/cm 450 0.98 5.9 -4.5 489 0.98 2.6 -1.4 490 0.99 2.7 -1.8 845 0.99 4.0 -2.5 895 0.98 1.1 -0.3
下载: 导出CSV
表 6 结合P351站雪深与SNOTEL站气候变量的雪水当量估计值与原位观测值对比
Table 6 Comparison Between the SWEs Predicted with Snow Depth of P351 and SNOTEL Climate Data and the In-situ SWEs
测站 与P351站高差/m 距P351距离/km R2 RMSE/cm Bias/cm 450 -126.19 30.86 0.91 8.4 -3.6 489 -415.75 3.81 0.86 7.6 5.7 490 -16.46 0.41 0.98 3.4 -1.5 845 38.40 13.39 0.93 16.4 -14.1 895 -769.32 27.24 0.57 13.9 11.7
下载: 导出CSV
表 7 结合P351站雪深与ClimateNA站气候变量的雪水当量估计值与原位观测值对比结果
Table 7 Comparison Between the SWEs Predicted with Snow Depth of P351 and ClimateNA Climate Data and the In-situ SWEs
测站 与P351站高差/m 距P351距离/km R2 RMSE/cm Bias/cm 450 -126.19 30.86 0.91 8.6 -4.0 489 -415.75 3.81 0.89 5.8 4.2 490 -16.46 0.41 0.98 4.0 -2.2 845 38.40 13.39 0.91 17.9 -15.3 895 -769.32 27.24 0.59 13.6 11.5
下载: 导出CSV
-
[1] Walsh J E. Snow Cover and Atmospheric Variability: Changes in the Snow Covering the Earth's Surface Affect both Daily Weather and Long-Term Climate[J]. American Scientist, 1984, 72(1): 50-57
[2] Beniston M. Climatic Change in Mountain Regions: A Review of Possible Impacts[M]. Dordrecht: Springer, 2003
[3] Elder K, Rosenthal W, Davis R E. Estimating the Spatial Distribution of Snow Water Equivalence in a Montane Watershed[J]. Hydrological Processes, 1998, 12(11): 1 793-1 808 doi: 10.1002/(SICI)1099-1085(199808/09)12:10/11%3C1793::AID-HYP695%3E3.0.CO;2-K/abstract
[4] McCreight J L, Small E E. Modeling Bulk Density and Snow Water Equivalent Using Daily Snow Depth Observations[J]. The Cryosphere, 2014, 8 (2): 521-536 doi: 10.5194/tc-8-521-2014
[5] Mizukami N, Perica S. Spatiotemporal Characteristics of Snowpack Density in the Mountainous Regions of the Western United States[J]. Journal of Hydrometeorology, 2008, 9(6): 1 416-1 426 doi: 10.1175/2008JHM981.1
[6] Larson K M, Gutmann E D, Zavorotny V U, et al. Can We Measure Snow Depth with GPS Receivers? [J]. Geophysical Research Letters, 2009, 36(17): L17502 doi: 10.1029/2009GL039430
[7] Larson K M. GPS Interferometric Reflectometry: Applications to Surface Soil Moisture, Snow Depth, and Vegetation Water Content in the Western United States[J]. Wiley Interdisciplinary Reviews: Water, 2016, 3(6): 775-787 doi: 10.1002/wat2.1167
[8] Hefty J, Gerhatova L. Experience with Snow Depth Sensing Using GPS Multipath Gained at Some Permanent Stations in Central Europe[C]//The EGU General Assembly Conference Abstracts, Vienna, Austria, 2013
[9] Jin S G, Najibi N. Sensing Snow Height and Surface Temperature Variations in Greenland from GPS Reflected Signals[J]. Advances in Space Research, 2014, 53(11): 1 623-1 633 doi: 10.1016/j.asr.2014.03.005
[10] Gutmann E D, Larson K M, Williams M W, et al. Snow Measurement by GPS Interferometric Reflectometry: An Evaluation at Niwot Ridge, Colorado[J]. Hydrological Processes, 2012, 26(19): 2 951- 2 961 doi: 10.1002/hyp.8329
[11] Larson K M, Nievinski F G. GPS Snow Sensing: Results from the EarthScope Plate Boundary Observatory[J]. GPS Solutions, 2013, 17(1): 41-52 doi: 10.1007/s10291-012-0259-7
[12] Gleason S, Gebre-Egziabher D. GNSS Applications and Methods[M]. London: Artech House, 2009
[13] De Michele C, Avanzi F, Ghezzi A, et al. Investigating the Dynamics of Bulk Snow Density in Dry and Wet Conditions Using a One-Dimensional Model[J]. The Cryosphere, 2013, 7(2): 433-444 doi: 10.5194/tc-7-433-2013
[14] Meløysund V, Leira B, Høiseth K V, et al. Predicting Snow Density Using Meteorological Data[J]. Meteorological Applications, 2007, 14(4): 413-423 doi: 10.1002/met.40
[15] Hill D F, Burakowski E A, Crumley R L, et al. Converting Snow Depth to Snow Water Equivalent Using Climatological Variables[J]. The Cryosphere, 2019, 13(7): 1 767-1 784 doi: 10.5194/tc-13-1767-2019
[16] Jacobson M D. Estimating Snow Water Equivalent for a Slightly Tilted Snow-Covered Prairie Grass Field by GPS Interferometric Reflectometry[J]. EURASIP Journal on Advances in Signal Processing, 2014, 2 014(1): 1-9 doi: 10.1186/1687-6180-2014-61?site=asp.eurasipjournals.springeropen.com
[17] McCreight J L, Small E E, Larson K M. Snow Depth, Density, and SWE Estimates Derived from GPS Reflection Data: Validation in the Western US[J]. Water Resources Research, 2014, 50(8): 6 892- 6 909 doi: 10.1002/2014WR015561
[18] Arogundade A B, Qualls R J. Is Climate Really Changing? Visible Signs of Climate Change in Idaho [C]//The 38th Annual Convention of the National Society of Black Engineers, Pittsburgh, Pennsylvania, America, 2012
[19] Herring T A, Melbourne T I, Murray M H, et al. Plate Boundary Observatory and Related Networks: GPS Data Analysis Methods and Geodetic Products[J]. Reviews of Geophysics, 2016, 54(4): 759-808 doi: 10.1002/2016RG000529
[20] Wang T L, Hamann A, Spittlehouse D L, et al. ClimateWNA -High-Resolution Spatial Climate Data for Western North America[J]. Journal of Applied Meteorology and Climatology, 2012, 51 (1): 16-29 doi: 10.1175/JAMC-D-11-043.1
[21] Wang T L, Hamann A, Spittlehouse D, et al. Locally Downscaled and Spatially Customizable Climate Data for Historical and Future Periods for North America[J]. PLoS One, 2016, 11(6) : e0156720 doi: 10.1371/journal.pone.0156720
[22] 吴继忠, 王天, 吴玮. 利用GPS-IR监测土壤含水量的反演模型[J]. 武汉大学学报·信息科学版, 2018, 43(6): 887-892 doi: 10.13203/j.whugis20160088 Wu Jizhong, Wang Tian, Wu Wei. Retrieval Model for Soil Moisture Content Using GPS-Interferometric Reflectometry[J]. Geomatics and Information Science of Wuhan University, 2018, 43(6): 887-892 doi: 10.13203/j.whugis20160088
[23] Roesler C, Larson K M. Software Tools for GNSS Interferometric Reflectometry (GNSS-IR)[J]. GPS Solutions, 2018, 22(3): 1-10 http://smartsearch.nstl.gov.cn/paper_detail.html?id=086f05c89c7e7a22cc60a6a86ad6711a
[24] Larson K M, Small E E, Gutmann E, et al. Using GPS Multipath to Measure Soil Moisture Fluctuations: Initial Results[J]. GPS Solutions, 2008, 12 (3): 173-177 doi: 10.1007/s10291-007-0076-6
[25] Lundberg A, Richardson-Näslund C, Andersson C. Snow Density Variations: Consequences for Ground-Penetrating Radar[J]. Hydrological Processes, 2006, 20(7): 1 483-1 495 doi: 10.1002/hyp.5944
[26] Jonas T, Marty C, Magnusson J. Estimating the Snow Water Equivalent from Snow Depth Measurements in the Swiss Alps[J]. Journal of Hydrology, 2009, 378(1/2): 161-167 http://www.researchgate.net/profile/Christoph_Marty/publication/222813786_Estimating_the_snow_water_equivalent_from_snow_depth_measurements_in_the_Swiss_Alps._J_Hydrol/links/0f31752ebb793cb8d1000000.pdf
[27] Hill D F, Burakowski E A, Crumley R L, et al. Converting Snow Depth to Snow Water Equivalent Using Climatological Variables[J]. The Cryosphere, 2019, 13(7): 1 767-1 784 doi: 10.5194/tc-13-1767-2019
[28] Gutmann E D, Larson K M, Williams M W, et al. Snow Measurement by GPS Interferometric Reflectometry: An Evaluation at Niwot Ridge, Colorado[J]. Hydrological Processes, 2012, 26(19): 2 951- 2 961 doi: 10.1002/hyp.8329
[29] Nievinski F G, Larson K M. Inverse Modeling of GPS Multipath for Snow Depth Estimation-Part I: Formulation and Simulations[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52 (10): 6 555-6 563 doi: 10.1109/TGRS.2013.2297681
[30] Sturm M, Taras B, Liston G E, et al. Estimating Snow Water Equivalent Using Snow Depth Data and Climate Classes[J]. Journal of Hydrometeorology, 2010, 11(6): 1 380-1 394 doi: 10.1175/2010JHM1202.1
[31] Dyer J L, Mote T L. Spatial Variability and Trends in Observed Snow Depth over North America[J]. Geophysical Research Letters, 2006, 33 (16) : L16503 doi: 10.1029/2006GL027258
-
期刊类型引用(1)
1. 赵智辉，杨金虎. 基于GNSS技术对滑坡体稳定性试验研究. 山西建筑. 2024(03): 86-89 . 
百度学术
其他类型引用(0)
计量
- 文章访问数: 1488
- HTML全文浏览量: 476
- PDF下载量: 106
- 被引次数: 1
