GNSS real-time cycle slip repair algorithm and PPP test analysis
-
摘要: 全球卫星导航系统(GNSS)能够为用户提供定位、导航和授时(PNT)服务，被广泛应用于国防安全保障和国民经济建设中. 实时精密单点定位(RT-PPP)是一种高精度卫星导航定位方法，针对信号中断导致的重新初始化时间长的问题，提出一种基于历元间差分伪距和载波相位观测量的周跳修复算法，设计了一种RT-PPP算法，并介绍其实现流程. 采用国际GNSS服务(IGS)观测站数据进行周跳修复实验，成功率在99%以上，缩短了重新收敛时间. 在楼顶进行推车实验，RT-PPP在水平方向定位精度优于1 cm，高程定位方向精度为2~3 cm.
-
关键词:
- 全球卫星导航系统(GNSS) /
- 周跳修复 /
- 实时精密单点定位(RT-PPP) /
- 定位实验 /
- 精度评定
Abstract: Global Navigation Satellite System (GNSS) can provide users with positioning, navigation and timing services, and is widely used in national defense security and national economic construction. Real time precise point positioning (PPP) is a high precise satellite navigation positioning method, this paper discusses the current problems of the algorithm that long reinitialization time caused by signal interruption , and proposes a cycle slip repair algorithm based on the differential pseudorange and carrier phase observations. It designs a real-time precise point positioning algorithm and introduces its implementation process. The success rate of the cycle slip repair experiment using the data of the IGS observatory is more than 99%, and shorten the reconvergence time. The cart experiment was carrier out on the roof of the building, and the real-time precise point positioning accuracy was better than 1 cm in the horizontal direction, and the accuracy was 2~3 cm in the vertical direction. -
0. 引 言
全球卫星导航系统(GNSS)是一种利用卫星测距信号为用户接收机提供全天候、高精度的定位、导航和授时(PNT)服务的现代化系统，主要包括美国的GPS，俄罗斯的GLONASS，欧盟的Galileo和中国的北斗卫星导航系统(BDS)[1]. 目前，GNSS高精度定位算法包括差分定位算法和精密单点定位(PPP)算法[2]. 由于PPP算法不依赖地面参考站，能够独立解算接收机坐标、接收机钟差和大气延迟参数，被广泛应用于低轨卫星定轨、时间同步、大气监测等领域[3-6].
由于最终的精密星历通常在2周之后才能获得，传统的PPP一般采用事后处理模式. 为了满足工程领域和科研领域对高精度位置等信息的实时性需求，国内外相关学者和机构基于地面观测数据对卫星轨道和钟差进行实时估计，并开展实时PPP (RT-PPP)研究工作[7-9]. 目前，RT-PPP存在的主要问题有：1)由于PPP数学模型参数个数和误差项较多，以及实时轨道钟差精度有限等因素，一般需要几十分钟的初始化时间[10]；2)受观测值中的周跳影响，不连续的载波相位模糊度会导致定位结果重新初始化；3)采用方差阵和观测量残差对定位结果进行内符合精度评定，不能准确描述参数估计精度；4)收敛时间一般是将坐标估值与真实值对比得到的统计结果，在RT-PPP尤其是动态PPP应用中，真实坐标未知导致收敛时间的不确定性. 以上几种问题都在一定程度上降低了RT-PPP的实用性、稳定性和可靠性，这也是当前RT-PPP的应用瓶颈和GNSS数据处理领域的研究重点[11].
针对目前RT-PPP存在的重新收敛问题，本文基于国际GNSS服务(IGS)实时产品开展相关研究工作，首先介绍RT-PPP观测方程，并提出一种周跳修复算法，用以解决信号中断后重新收敛时间过长的问题. 结合实时产品和观测数据流，论述RT-PPP算法及其实现流程. 根据IGS观测站进行周跳修复实验，验证其对信号中断后的定位结果的改善程度. 同时，在楼顶进行轨道推车实验，进一步评定本文提出的RT-PPP算法所能达到的定位效果.
1. RT-PPP算法
1.1 观测方程
PPP是一种基于单个站观测数据，通过固定导航卫星轨道和钟差，使用模型改正固体潮、天线相位中心(APC)偏差等系统误差，独立解算高精度接收机坐标、接收机钟差、大气延迟和模糊度参数的卫星定位算法[12]. 根据使用的观测方程形式不同，可分为无电离层组合PPP和非组合PPP. 两者均采用双频伪距和载波相位观测值，最终的参数估值精度相当. 不同之处在于非组合PPP观测量精度高，能够同时估计电离层延迟，适用于多频数据处理. 无电离层组合PPP参数个数少，模型强度更高，并且与IGS发布的精密卫星轨道和钟差产品具有更好的一致性[13]. 本文采用无电离层组合PPP观测方程
$$ \left\{\begin{aligned} &P_r^s = \rho _r^s + c \cdot {\rm{d}}{t_r} - c \cdot {\rm{d}}t_p^s + m_r^s \cdot zt{d_r} + {\varepsilon _P} \\ &L_r^s = \rho _r^s + c \cdot {\rm{d}}{t_r} - c \cdot {\rm{d}}t_p^s + m_r^s \cdot zt{d_r} + B_r^s + {\varepsilon _L} \end{aligned}. \right. $$ (1) 式中：
$ P_r^s $ 和$ L_r^s $ 分别为无电离层伪距和载波相位观测量；$\; \rho _r^s$ 为接收机到卫星的几何距离，卫星位置为实时卫星轨道；${\rm{d}}{t_r}$ 为接收机钟差；${\rm{d}}t_p^s$ 为卫星钟差；$ c $ 为光速；$ m_r^s $ 为对流层延迟投影函数；$ zt{d_r} $ 为天顶(U)方向对流层延迟；$ B_r^s $ 为无电离层组合相位模糊度，不具备整数特性；$ {\varepsilon _P} $ 和$ {\varepsilon _L} $ 则分别为伪距和载波相位观测噪声.1.2 周跳修复算法
由于信号失锁等原因造成的周跳会使载波相位模糊度发生未知的整数跳变，从而造成坐标参数的重新收敛. 周跳在GNSS实际观测时难以避免，将导致坐标参数需要长时间重新收敛，这对RT-PPP应用尤为不利. 用于PPP周跳快速修复的方法包括联合惯性导航辅助周跳修复、使用多普勒观测数据的周跳修复方法[14-15]. 本文采用历元间差分伪距和载波相位观测值修复周跳的方法，来避免RT-PPP中的重新收敛问题. 按照数据处理的顺序，周跳修复算法的步骤可以分为：建立电离层延迟变化模型、固定宽巷模糊度、固定窄巷模糊度和修复周跳.
1.2.1 建立电离层延迟变化模型
电离层延迟变化模型主要用于辅助宽巷模糊度固定，因为电离层延迟是历元间差分宽巷相位观测方程中最主要的误差源，间隔30 s时的影响约为几个厘米，时差为2 min时能够达到分米级，影响宽巷模糊度固定. 因此，需要对每个卫星的电离层延迟变化进行建模和改正.
$$ {\begin{bmatrix} {{L_I}({t_1})} \\ {{L_I}({t_2})} \\ \vdots \\ {{L_I}({t_n})} \end{bmatrix}} = {\begin{bmatrix} {{t_1}}&1 \\ {{t_2}}&1 \\ \vdots & \vdots \\ {{t_n}}&1 \end{bmatrix}} {\begin{bmatrix} { - \gamma \dot I} \\ {Am{b_I}} \end{bmatrix}} . $$ (2) 无几何距离相位观测值
${L_I}$ 中主要包含电离层变化$ - \gamma \dot I$ 和模糊度$Am{b_I}$ . 在连续弧段和较短时间范围内，模糊度及相位延迟等误差为常数，不影响${L_I}$ 时间序列中的电离层变化信息，可用于建立电离层变化模型. 使用${L_I}$ 序列线性拟合电离层延迟模型$\dot I$ ，并对一定时间间隔内的电离层延迟变化进行预测，当电离层环境复杂时，可采用二阶多项式拟合.1.2.2 固定宽项模糊度
基于电离层延迟变化模型，使用历元间差分窄巷伪距
${\dot P_N}$ 和宽巷相位${\dot L_W}$ 固定宽巷模糊度. 根据误差源的短期稳定性，式(3)中消除了伪距和相位硬件延迟以及潮汐改正等影响，包括相对坐标、接收机钟速和宽巷模糊度参数. 电离层延迟变化是主要误差项，通过第一步建立的电离层延迟变化模型改正，对流层延迟使用估计值或者模型改正，不再设置参数. 根据模糊度浮点解和方差阵，使用LAMBDA算法固定整周宽巷模糊度${N_1} - {N_2}$ .$$\left\{\begin{aligned} &{{\dot P}_N} - {{\dot I}_N} = \dot \rho + c{\rm{d}}{{\dot t}_r} - c{\rm{d}}{{\dot t}^s}\\ &{{\dot L}_W} + {{\dot I}_W} = \dot \rho + c{\rm{d}}{{\dot t}_r} - c{\rm{d}}{{\dot t}^s} + {\lambda _W}({N_1} - {N_2}) \end{aligned}.\right. $$ (3) 由于方程中使用了窄巷伪距组合，观测噪声为分米级，计算得到的坐标参数精度约也为分米级. 为了避免伪距噪声对模糊度浮点解的影响，理论上可以使用周跳前所有连续弧段内的数据，组成时间差分方程对模糊度参数进行卡尔曼滤波. 由于方程中系统误差项与时间间隔
$\Delta T$ 相关，时差越大，误差相关性越小，采用时间间隔权函数对不同时刻的观测值定权$$ W(\Delta T) = \left\{ {\begin{array}{*{20}{c}} {1.0},&{\Delta T \leqslant 30\;{\rm{s}}} \\ {\tan \left(\displaystyle\frac{\text{π} }{4} \cdot \frac{{30}}{{\Delta T}}\right)},&{\Delta T \gt 30\;{\rm{s}}} \end{array}} \right. . $$ (4) 1.2.3 固定窄项模糊度
基于电离层延迟变化模型以及宽巷模糊度，使用历元间差分宽巷相位
$ {\dot L_W} $ 和无电离层相位${\dot L_C}$ 固定窄巷模糊度. 此时，无模糊度的宽巷相位组合转换为厘米级的伪距观测值，避免了使用伪距观测值对波长较短的窄巷模糊度固定的不利影响.$$ \left\{\begin{aligned} &{{\dot L}_W} + {{\dot I}_W} - {\lambda _W}({N_1} - {N_2}) = \dot \rho + c{\rm{d}}{{\dot t}_r} - c{\rm{d}}{{\dot t}^s} \\ &{{\dot L}_C} - \frac{{{\lambda _N}}}{{{\lambda _2}}}{\lambda _W}({N_1} - {N_2}) = \dot \rho + c{\rm{d}}{{\dot t}_r} - c{\rm{d}}{{\dot t}^s} + {\lambda _N}{N_1} \end{aligned} .\right. $$ (5) 与宽巷模糊度固定类似，式(5)中消除了相关系统误差，仅包含相对坐标、钟速和模糊度参数，同样采用LAMBDA算法固定窄巷模糊度
${N_1}$ . 需要注意的是，在动态PPP模式下当接收机天线水平旋转且所有卫星都发生周跳时，相位缠绕误差、模糊度参数和钟差参数耦合. 为保持模糊度整数特性，需要引入一颗模糊度为整数的参考卫星. 此时，相位缠绕误差和模糊度基准误差会被钟差参数吸收，所有模糊度保持整数特性.1.2.4 修复周跳
将宽巷模糊度和窄巷模糊度转换为原始频率的周跳
${N_1}$ 和${N_2}$ ，对观测值进行修复，实现PPP周跳后的快速重新初始化. 为避免错误的周跳修复对参数估计造成不良影响，应该在修复后重新使用TurboEdit方法探测周跳并标记弧段.1.3 RT-PPP流程
与常规PPP相比，RT-PPP的不同之处在于实时精密星历和观测数据获取，以及对实时数据的处理方法. 本文设计了一种RT-PPP数据处理流程，如图1所示. RT-PPP数据处理流程主要包括：实时精密轨道和钟差获取、实时观测数据获取和PPP数据处理.
1)实时精密轨道和钟差获取
通过IGGNtrip获取实时数据分析中心播发的轨道钟差改正数和广播星历. IGGNtrip模块包括网络数据接收、RTCM解码、数据存储和统计功能，采用共享内存机制将数据接收和数据处理分离[16]. RT-PPP主程序调用IGGNtrip API返回实时产品和广播星历数据流，并将这两种数据恢复成PPP可用的实时精密轨道和钟差.
2)实时观测数据获取
实时观测数据按照来源可以分为网络数据流和串口数据流. 网络数据流是由IGGNtrip提供的IGS站实时观测数据，串口数据流是通过串口工具获取的GNSS接收机实时观测数据. 串口在传输数据时采用特定比特率发送和接收一定数据长度的二进制数据，并不是一个时刻完整的观测数据，这就需要建立缓存机制和循环冗余校验来确保观测数据完整性. 目前RT-PPP软件能够接收10 Hz采样的高频数据，支持的数据格式包括BINEX (BINary EXchange)和RTCM.
3) PPP数据处理
PPP数据处理是RT-PPP的核心，主要包括实时观测数据和实时精密星历时间同步、观测值定权、PPP数据处理和收敛时间判断. 其中，时间同步是针对当前实时产品相对于实时观测数据5~20 s的时间延迟，采用分层式实时产品时延处理方法，实现精密星历与观测值同步. 观测值定权是在当前实时产品精度有限的情况下，不使用观测值的标称精度，而采用事先统计的实时卫星钟差精度确定观测值先验精度和不同系统之间的权重关系. PPP数据处理中，采用卡尔曼滤波算法逐历元估计测站坐标、对流层延迟和浮点解模糊度. 基于GPT2w经验对流层模型约束天顶方向对流层湿分量延迟参数，加速RT-PPP初始化过程，并采用周跳修复算法来解决重新初始化问题. 另外，由于RT-PPP应用中不能通过与真实坐标对比来确定收敛时间，而统计的收敛时间跟观测条件有关，具有一定不确定性. 根据模糊度参数是否收敛判断坐标收敛时间，解决RT-PPP尤其是动态PPP中收敛时间不明确的问题.
2. 实时周跳修复实验
为了验证RT-PPP周跳修复算法的可行性，本节采用2020年第250 天JFNG站和GMSD站30 s采样间隔的BDS/GPS/GLONASS数据进行周跳修复实验. 宽巷模糊度和窄巷模糊度固定成功的ratio阈值都设为3.0，周跳修复成功率定义为宽巷和窄巷模糊度固定成功率. 采用CNES 的CLK92实时产品作为精密星历，进行实时数据处理. 对当天13时至23时所有时刻的观测数据都设置一个数值为0的周跳，对这段时间的数据尝试修复周跳，统计模糊度固定ratio值和周跳修复成功率.
图2中，蓝色点代表宽巷模糊度固定ratio阈值，红色点代表窄巷模糊度固定ratio阈值，绿色线为阈值，都以倒数的形式来表示. 当蓝色点和红色点低于绿色的分界线时，代表模糊度固定成果，反之则代表模糊度固定失败. 从图2中可以看出，JFNG站的宽巷模糊度固定成功率略低于GMSD站，这主要是因为JFNG站附近高大的树木较多，伪距多路径效应严重. 除此之外，两个站点的模糊度固定成功率都在99%以上，ratio阈值超过10.0的观测时刻占约为97%，窄巷伪距残差为0.26 m，与先验值相当.
采用JFNG站当天BDS/GPS/GLONASS观测数据，以及CLK92产品模拟PPP-RTK周跳修复. 在14时至22时每个整点时刻人为设置一个数值为卫星号乘以小时数的周跳，并对GLONASS卫星加50，BDS卫星加100. 图3为周跳修复前和周跳修复后的动态PPP效果对比.
从图3可以看出，在整点时刻人为引入周跳后，如果只对周跳进行探测并不修复，定位结果需要不断地重新初始化. 经统计，未采用周跳修复算法的RT-PPP在信号中断后，平均收敛时间为8.5 min. 而如果使用本文所论述的方法对周跳进行修复，RT-PPP定位结果则不需要重新初始化，只需要在初始化阶段完成对坐标参数的滤波收敛. 因此，本文所提出的周跳修复算法能够有效解决RT-PPP在信号发生中断后的重新收敛问题.
3. RTK定位实验
为了验证RT-PPP在天线运动状态下的定位效果，本文在某楼顶进行轨道推车实验. 位于楼顶的轨道为闭合环线，周长约50 m，轨道面上有一个可以人为拖动的推车. 在轨道推车实验中，人为拖动的推车在直线处的平均移动速度约为1.0 m/s，在转弯处的平均速度约为0.2 m/s. 实验所采用的天线类型为抗多路径的TRM59900，安装在位于推车平台中一个三脚架上. GNSS接收机类型为大地测量型TRIMBLE NETR9，同时记录BDS/GPS/GLONASS数据，采样频率为1 Hz. GNSS接收机通过串口与运行在笔记本电脑上的RT-PPP程序连接,数据传输格式为BINEX.
RT-PPP使用串口数据流的GNSS观测数据，并采用CLK90产品. 由于卫星轨道参考中心位于APC，不需要进行APC改正. 同时，在流动站附近的一个观测墩上安装同样型号的GNSS接收机和天线作为基准站，与推车上的流动站进行实时动态(RTK)处理，作为RT-PPP坐标结果的参考值. 图4为轨道推车实验现场，图5为RTK和RT-PPP处理的流动站运动轨迹对比.
轨道推车实验时间为北京时间2020年5月16日14时至16时，在15时28分左右人为拖动推车平台顺时针移动，15时30分左右停止，其余时间平台保持静止. 流动站PPP-RTK在18.5 min后收敛，单频RTK在首个历元就能够成功固定模糊度，ratio阈值最小为3.2. 将流动站RTK结果转化为ITRF下的绝对坐标，与RT-PPP定位结果位于统一的坐标系统. RT-PPP和RTK的定位结果均显示流动站在15时38分29秒开始移动，并于30分14秒停止. 图6为15时20分至35分时间段内，RT-PPP与RTK相比在三维方向上的坐标偏差. 可以看出在静止阶段和运动阶段，RT-PPP坐标偏差变化不明显. 轨道推车PPP-RTK定位精度如表1所示，静止阶段RT-PPP在东(N)、北(E)、U三维方向上坐标偏差均方根(RMS)分别为0.69 cm、0.65 cm、2.39 cm，运动阶段的坐标偏差RMS分别为0.79 cm、0.74 cm、3.01 cm.
表 1 不同阶段RT-PPP定位精度统计表cm 运动状态 N E U 静止阶段 0.69 0.65 2.39 运动阶段 0.79 0.74 3.01 4. 结束语
本文首先介绍了RT-PPP的研究意义和现状，及其在实际应用中存在的问题. 针对信号中断后，传统PPP算法需要进行重新初始化的问题，提出一种实时周跳修复算法. 基于IGS观测站实测数据模拟周跳修复效果，JFNG站和GMSD站周跳修复成功率都在99%以上，周跳修复成功后，能够避免长达约10 min的重新初始化时间. 基于IGGNtrip模块获取实时精密星历，并通过串口获取接收机二进制数据，进行楼顶推车实验. 结果表明，将RTK结果作为参考，RT-PPP在水平方向定位精度优于1 cm，高程方向为2~3 cm. 综上，本文提出的一种周跳修复的RT-PPP算法能够有效解决信号中断后的重新初始化问题，并且能够达到厘米级动态定位效果，在相关科研和工程应用中具有推广价值.
随着卫星导航系统的逐步完善和更新，以及用户对定位性能需求的不断提升，PPP算法也在不断发展. 相关研究主要聚焦在以下几个方面：1) PPP-RTK算法能够显著提高PPP定位精度，并缩短其收敛时间. 该算法首先基于卫星端相位偏差产品实现用户端模糊度固定，并采用参考站网络内插的大气延迟来缩短模糊度首次固定时间，使PPP达到与RTK相同的定位效果，是目前GNSS精密定位算法领域的研究热点，正处于实验研究到工程应用的转换阶段[17]. 2)低轨导航卫星增强PPP是当前论证的热点，美国的“铱星二代”、中国的“鸿雁”、等低轨卫星具有轨道高度低、运动快的特点. 卫星几何图形变化快，历元间方程相关性低，使得低轨卫星联合导航卫星能够将PPP的收敛时间从数十分钟缩短到几秒钟，从而真正发挥PPP技术的优越性[18]. 3) PPP与惯性导航、激光雷达、摄像头、伪卫星、超宽带(UWB)等定位手段进行组合，能够有效解决GNSS信号拒止条件下，如城市峡谷、隧道、室内等环境中连续定位问题，提高定位服务的连续性、可靠性和完好性，在无人驾驶、人工智能领域具有重要应用价值[19].
-
[1] MONTENBRUCK O, STIGENBERGER P, PRANGE L, et al. The Multi-GNSS Experiment (MGEX) of the International GNSS Service (IGS) – achievements, prospects and challenges[J]. Advances in space research, 2017, 59(7): 1671-1697. DOI: 10.1016/j.asr.2017.01.011
[2] ZUMBERGE J F, HEFLIN M B, JEFFERSOND C, et al. Precise point positioning for the efficient and robust analysis of GPS data from large networks[J]. Journal of geophysical research solid earth, 1997, 102(B3): 5005-5017. DOI: 10.1029/96JB03860
[3] 李星星, 张小红, 李盼. 固定非差整数模糊度的PPP快速精密定位定轨[J]. 地球物理学报, 2012, 55(3): 833-840. [4] 施闯, 张东, 宋伟, 等. 北斗广域高精度时间服务原型系统[J]. 测绘学报, 2020, 49(3): 269-277. [5] 郭敏, 张捍卫, 夏朋飞. GNSS天顶对流层延迟的短时天气预报分析[J]. 测绘科学, 2021, 46(4): 28-36. DOI: 10.16251/j.cnki.1009-2307.2021.04.005 [6] 王伟, 臧文驰, 彭竞, 等. 基于RT-PPP的低轨卫星实时高精度时间同步方法[J]. 全球定位系统, 2021, 46(5): 26-32. DOI: 10.12265/j.gnss.2021050301 [7] 楼益栋. 导航卫星实时精密轨道与钟差确定[D]. 武汉: 武汉大学, 2008. [8] 刘腾. 多模GNSS非组合精密单点定位算法及其电离层应用研究[D]. 北京: 中国科学院大学, 2017. [9] XIAO G W, LIU G Y, OU J K, et al. MG-APP: an open-source software for multi-GNSS precise point positioning and application analysis[J]. GPS solutions, 2020, 24(3). DOI: 10.1007/s10291-020-00976-1
[10] 卢福康, 余学祥, 肖星星, 等. BDS-2/BDS-3/GPS/Galileo双频精密单点定位精度分析与评价[J]. 全球定位系统, 2022, 47(5): 35-44. [11] 张小红, 胡家欢, 任晓东. PPP/PPP-RTK新进展与北斗/GNSS PPP定位性能比较[J]. 测绘学报, 2020, 49(9): 1084-1100. [12] KOUBA J, P HEROUX. Precise point positioning using IGS orbit and clock products[J]. GPS solutions, 2001, 5(2): 12-28. DOI: 10.1007/PL00012883
[13] 张小红, 左翔, 李盼. 非组合与组合PPP模型比较及定位性能分析[J]. 武汉大学学报(信息科学版), 2013, 38(5): 561-565. [14] 刘帅, 孙付平, 张伦东, 等. INS辅助周跳修复以实现精密单点定位瞬时重新收敛[J]. 中国惯性技术学报, 2015, 23(5): 607-614. DOI: 10.13695/j.cnki.12-1222/o3.2015.05.010 [15] DING W W, OU J K. Instantaneous re-initialization of real time kinematic PPP by adding doppler observation[J]. Journal of astronautics, 2013, 34(6): 750-800. DOI: 10.3873/j.issn.1000-1328.2013.06.008
[16] 贺伟欣. GNSS实时数据获取与精密单点定位研究[D]. 北京: 中国科学院大学, 2015. [17] KHODABANDEH A. Single-station PPP-RTK: correction latency and ambiguity resolution performance[J]. Journal of geodesy, 2021, 95(4): 42. DOI: 10.1007/s00190-021-01490-z
[18] 王磊, 李德仁, 陈锐志, 等. 低轨卫星导航增强技术——机遇与挑战[J]. 中国工程科学, 2020, 22(2): 144-152. [19] 朱锋. GNSS/SINS/视觉多传感器融合的精密定位定姿方法与关键技术[D]. 武汉: 武汉大学, 2019. -
期刊类型引用(1)
1. 王群，吴俊烽，肖名豪，黄亮星. 基于GNSS的安全授时系统研究. 电子元器件与信息技术. 2024(01): 30-33 . 
百度学术
其他类型引用(1)
下载:
计量
- 文章访问数: 408
- HTML全文浏览量: 93
- PDF下载量: 41
- 被引次数: 2
