留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

构建阈值模型改善TurboEdit实时周跳探测

张小红 曾琪 何俊 康超

张小红, 曾琪, 何俊, 康超. 构建阈值模型改善TurboEdit实时周跳探测[J]. 武汉大学学报 ● 信息科学版, 2017, 42(3): 285-292. doi: 10.13203/j.whugis20150045
引用本文: 张小红, 曾琪, 何俊, 康超. 构建阈值模型改善TurboEdit实时周跳探测[J]. 武汉大学学报 ● 信息科学版, 2017, 42(3): 285-292. doi: 10.13203/j.whugis20150045
ZHANG Xiaohong, ZENG Qi, HE Jun, KANG Chao. Improving TurboEdit Real-time Cycle Slip Detection by the Construction of Threshold Model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(3): 285-292. doi: 10.13203/j.whugis20150045
Citation: ZHANG Xiaohong, ZENG Qi, HE Jun, KANG Chao. Improving TurboEdit Real-time Cycle Slip Detection by the Construction of Threshold Model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(3): 285-292. doi: 10.13203/j.whugis20150045

构建阈值模型改善TurboEdit实时周跳探测

doi: 10.13203/j.whugis20150045
基金项目: 

中电集团54所高校合作项目 KX132600031

预研基金 51324040103

预研基金 9140A24020713JB11342

详细信息
    作者简介:

    张小红, 博士, 教授, 主要研究方向为GNSS精密单点定位及其地学应用。xhzhang@sgg.whu.edu.cn

  • 中图分类号: P228

Improving TurboEdit Real-time Cycle Slip Detection by the Construction of Threshold Model

Funds: 

The CET 54 Universities Coorperation Project KX132600031

the Pre-research Fund 51324040103

the Pre-research Fund 9140A24020713JB11342

More Information
    Author Bio:

    ZHANG Xiaohong, PhD, professor, speciallizes in GNSS precise point positioning and its application in Geoscience. E-mail:xhzhang@sgg.whu.edu.cn

图(9) / 表(1)
计量
  • 文章访问数:  1812
  • HTML全文浏览量:  90
  • PDF下载量:  659
  • 被引次数: 0
出版历程
  • 收稿日期:  2015-09-14
  • 刊出日期:  2017-03-05

构建阈值模型改善TurboEdit实时周跳探测

doi: 10.13203/j.whugis20150045
    基金项目:

    中电集团54所高校合作项目 KX132600031

    预研基金 51324040103

    预研基金 9140A24020713JB11342

    作者简介:

    张小红, 博士, 教授, 主要研究方向为GNSS精密单点定位及其地学应用。xhzhang@sgg.whu.edu.cn

  • 中图分类号: P228

摘要: 在比较现有各类TurboEdit改进算法的基础上,选取了一套适合实时周跳探测的组合方法,并针对其缺陷进行优化改进,加入阈值模型。在滑动平均滤波MW(Melborne-Wubbena)组合周跳探测方法中加入了自适应观测序列均方根变化的阈值模型,有效减少了周跳探测的误判和漏判;在相邻历元求差的GF(geometry-free)组合周跳探测中加入了随采样率及卫星高度角变化的加权阈值模型,显著提升了采样率较低时周跳探测的可靠性。大量数据测试检验表明,该组合方法有效可行,相比其他改进方法不仅在探测精度上有一定提高,在避免周跳误判、漏判上也有较大改善。

English Abstract

张小红, 曾琪, 何俊, 康超. 构建阈值模型改善TurboEdit实时周跳探测[J]. 武汉大学学报 ● 信息科学版, 2017, 42(3): 285-292. doi: 10.13203/j.whugis20150045
引用本文: 张小红, 曾琪, 何俊, 康超. 构建阈值模型改善TurboEdit实时周跳探测[J]. 武汉大学学报 ● 信息科学版, 2017, 42(3): 285-292. doi: 10.13203/j.whugis20150045
ZHANG Xiaohong, ZENG Qi, HE Jun, KANG Chao. Improving TurboEdit Real-time Cycle Slip Detection by the Construction of Threshold Model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(3): 285-292. doi: 10.13203/j.whugis20150045
Citation: ZHANG Xiaohong, ZENG Qi, HE Jun, KANG Chao. Improving TurboEdit Real-time Cycle Slip Detection by the Construction of Threshold Model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(3): 285-292. doi: 10.13203/j.whugis20150045
  • 随着用户对导航定位的精度和实时性要求越来越高,精密单点定位 (precise point positioning,PPP) 技术和实时动态差分测量 (real-time kinematic,RTK) 技术已成为GNSS领域的主要研究对象。这两种技术都是以载波相位观测值作为主要观测量,而周跳的发生会造成观测值序列的中断,影响模糊度的固定,最终导致定位性能的降低。

    自20世纪80年代以来,许多周跳探测方法陆续被提出,包括电离层残差法、高次差法、卡尔曼滤波法[1-3]等。但这些周跳探测方法大都基于事后处理,且均存在一定的局限性,TurboEdit方法[4]联合MW (Melborne-Wubbena) 组合和GF (geometry-free) 组合,因其具有单站探测、探测效率高、灵活性好等特点而被广泛采用,然而该方法也存在一些缺陷,如采用了噪声较大的伪距观测值、阈值固定以及拟合方式引入人为误差等。近年来,不少学者针对TurboEdit方法的缺陷,提出了一些改进的方法。方荣新、袁玉斌等采用多项式拟合对GF组合进行滤波[5, 6],该方法对拟合历元个数及拟合阶数要求较高,同时也会带来拟合误差和舍入误差;王振杰等提出了基于MW组合滑动前后向窗口和修正 (phase ionospheric residual, PIR) 组合历元一阶差商滑动窗口的“两步法”探测周跳[7],然而,该方法不仅引入了新的检验量,而且后向窗口只能用于事后处理;吴继忠等引入了差分无电离层观测值,且每个历元要计算卫星位置及站星距离, 增加了计算量[8];Liu提出了联合使用MW组合和总电子含量变化率 (total electron contents rate, TECR) 来探测周跳,但该方法对观测数据采样率要求较高[9];文献[10-11]提出了联合前后向滑动窗口的MW组合和电离层残差组合二次历元间差分 (second-order time-difference phase ionopheric residual,STPIR) 的周跳探测和修复方法,均只能用于事后处理。

    本文在总结比较各类TurboEdit改进算法的基础上,选取了一套适合于实时周跳探测的组合算法,针对MW组合和GF组合周跳探测构建自适应的阈值模型,以期提升周跳探测的性能,并通过大量数据测试比较验证了本文方法的优越性。

    • TurboEdit方法是利用双频载波相位观测值和伪距观测值构造的MW组合和GF组合观测值来联合探测周跳。

    • MW组合观测模型为:

      (1)

      式中,φ1φ2P1P2分别为f1f2频段上的相位观测值和伪距观测值;λWL=c/(f1-f2) 和NWL=N1-N2分别表示宽巷波长和宽巷模糊度。MW周跳检验量则为:

      (2)

      使用MW方法进行周跳探测时,还需要不断递推计算前i个历元的平均宽巷模糊度以及均方根,递推计算为:

      (3)
      (4)

      其中,NWL为前i个历元宽巷模糊度的平均值;σ(i) 表示前i个历元的标准差,若满足如下两个条件,则认为当前历元存在周跳[4]

      (5)
      (6)
    • 采用载波相位观测值和伪距观测值的GF组合模型为:

      (7)
      (8)

      式中,LGFPGF分别表示载波和伪距的GF观测值;ΔI表示两个频率上的电离层延迟之差。

      由于伪距观测值中含有约30 cm的噪声影响,因此, 为了保证周探测的准确性与稳定性,需要对PGF(i) 进行多项式拟合生成Qi来滤去伪距噪声的影响,多项式拟合阶数必须满足条件m=min[(N/100+1), 6],N为历元个数,令式 (7) 减去式 (8),得到GF组合周跳探测模型:

      (9)

      对于GF组合法,若满足如下两个条件则认为当前历元存在周跳[4]

      (10)
      (11)

      上述方法存在以下几点弊端:①卫星在刚进入和即将离开视场时,卫星高度角较低导致大气延迟误差和多路径误差的影响较大,如果把这部分数据加入到中间弧段宽巷模糊度的均值和均方根的计算中,会导致计算出的均方根较大,从而可能出现周跳漏判;同理也会造成后半弧段出现周跳误判;②在GF组合周跳探测中对伪距电离层残差组合进行了多项式拟合,而多项式拟合涉及多项式类型、拟合次数等人为因素,降低了数据处理的准确性和客观性;③该算法的探测性能依赖于电离层的变化,即对采样率要求较高,若接收机失锁时间较长也可能造成该方法失效。

    • 针对TurboEdit周跳探测方法的缺陷,近几年已有许多关于TurboEdit改进算法的研究,下面简要介绍几种具有代表性的改进方法。

    • 1) 滑动平均滤波模型

      在MW方法中采用滑动平均的方法替代宽巷模糊度中对每一历元之前的所有数据整体求平均的方法[6]。即将式 (3) 和式 (4) 替换成式 (12):

      (12)
      (13)

      其中,m表示参加平均值计算的滑动窗口宽度,本文中m=min (N/10,25),N表示弧段内历元个数。上式的前提是i≥m,若i < m,则只计算前i个历元的宽巷模糊度均值和均方根。

      2) 前后向滑动窗口模型

      前后向滑动窗口模型是对当前探测历元前后一定宽度窗口内的宽巷模糊度进行平滑滤波处理[7, 10],从而大大削弱了观测噪声对周跳误判的影响。前后向滑动窗口内宽巷模糊度均值以及周跳检验量的计算公式为:

      (14)
      (15)
      (16)

      其中, NWL, Bwd(i) 和NWL, Fwd(i) 分别表示待检测历元后向m个历元和前向n个历元的宽巷模糊度均值。在探测之前要确保后向窗口内无周跳发生, 因此, 要对后向窗口作逆向数据剔除处理。周跳探测阈值一般根据观测条件及采样率选为一常数,也可以像原TurboEdit方法那样求出均方根再进行判断,此处不再赘述。

      上述两种改进方法都是基于滑动窗口滤波的思想,有效避免了原TurboEdit方法中因为弧段内数据质量前后不同对周跳探测的影响。其中,前后向滑动窗口模型采用双向平滑滤波,理论较为严密,探测性能较好,但算法稍复杂,而且需要探测历元前后的数据,适用于事后处理;滑动平均滤波模型只对探测历元之前的部分数据进行平滑滤波,适用于实时周跳探测,但NWL(i) 是由当前历元的MW组合求得,相比平滑过的NWLσ(i) 而言,对多路径效应和观测噪声的影响更为敏感,容易造成周跳误判。

    • 1) 相邻历元求差法

      相邻历元求差法去除了噪声较大的伪距观测值,仅将相邻两个历元的相位GF组合观测值求差来探测周跳,表达式如下:

      (17)

      式中,ΔI′=ΔI(i)-ΔI(i-1) 表示相邻历元间电离层延迟之差。一般情况下,电离层变化缓慢,ΔI′近似为0,可以根据ΔLGF(t) 来判断周跳。

      2) 滑动窗口拟合模型

      在无周跳发生时,MW组合近似为一常数,而GF组合则为光滑的曲线。在GF组合周跳探测时可以采用一种基于滑动窗口的曲线拟合方法[6],即从第i个历元开始取该历元之前n个历元的数据进行拟合,求出第i个历元的拟合值QGF(i) 和拟合标准差σ(i),如果满足:

      (18)

      则认为当前历元i异常,出现了周跳或粗差。

      上述两种方法相对于原TurboEdit方法而言,去除了噪声较大的伪距观测值,探测精度得到极大提升。其中,相邻历元求差法仅采用相邻两个历元的观测信息来探测周跳,简单易行,探测精度高,但其有效性依赖于电离层的变化;滑动窗口拟合模型可以看作是原方法和相邻历元求差法在滤波长度上的扩展,该方法从理论上来说较严密,但该方法的性能受到拟合长度和拟合阶数的限制,相位GF组合的探测精度本身较高,若拟合长度和拟合阶数选取得不合适,将会影响周跳的探测。

    • 本文基于MW组合的滑动滤波模型和GF组合的相邻历元求差法来联合探测周跳。不管是MW组合的滑动滤波周跳探测还是GF组合的历元差探测周跳,阈值的选择对周跳探测的准确性影响至关重要。已有方法通常选取经验阈值,算法的适应性较差。为此,本文提出构建自适应的阈值模型,来提升周跳探测的整体性能。

    • 在改进的TurboEdit方法中,令宽巷模糊度变化量与均方根的比值为:

      (19)

      则MW组合探测周跳的关键条件|D|≥k,其中,k一般默认为4。采用该方法进行实测数据的周跳探测时,由于只对一侧的数据进行平滑滤波,因此可能会将波动较大的噪声误判为周跳。图 1给出了一组采样率为1 s的静态观测数据的MW组合周跳探测情况,原始观测数据无周跳,图中曲线为截取一小时内比值D的变化情况,两侧的直线表示原阈值k0=±4,若曲线波动超过阈值,则该算法误判为当前历元发生周跳。

      图  1  D值变化曲线

      Figure 1.  Change of D Value

      图 1可以看出,大约有10个左右的历元会做出周跳误判,结合按式 (4) 推求的均方根发现,周跳误判历元的均方根值均较小,此时曲线上可能由噪声引起的不到半周的波动都会使式 (19) 成立,从而引起周跳误判。因此, 需要重新构建阈值模型,使之能够有效减少周跳误判的发生。为了得到具有统计意义、准确可靠的阈值模型,本文测试了几十组不同接收机、不同采样率下的观测数据,其中已事先在每颗卫星的数据中添加了3~5组两周周跳,计算每组数据中的|D|值与均方根,并将所有数据中的周跳误判历元和添加周跳历元的|D|值和均方根合并一起统计。不同采样率下的统计结果如图 2~图 4所示,“均方根”是指当前历元计算的窗口内宽巷模糊度均方差。其中红色方块表示周跳误判历元的探测情况;绿色圆点表示添加两周周跳历元的探测情况。

      图  2  1 s采样率下|D|值变化曲线

      Figure 2.  Change of|D| Value at 1 s Interval

      图  3  15 s采样率下|D|值变化曲线

      Figure 3.  Change of |D|Value at 15 s Interval

      图  4  30 s采样率下|D|值变化曲线

      Figure 4.  Change of |D| Value at 30 s Interval

      加入阈值模型之前,阈值设为定值4,即|D|值大于4时认为发生周跳,小于4时则认为无周跳。图 2~图 4中的红色方块表示无周跳但|D|值超过4的历元,即周跳误判。由图中可以看出,周跳误判主要集中在均方根较小的历元,且|D|值随均方根的增大而减小;当均方根大于0.4时,周跳误判率明显降低;当均方根大于0.6时,误判率基本为0。图中绿色圆点表示添加周跳的历元,可以看出,当均方根较小时,周跳历元|D|值较大,容易探测;当均方根位于0.4~0.6之间时,添加周跳历元与周跳误判历元较接近;当均方根大于0.6时,周跳历元的|D|值逐渐减小,少数历元的|D|值甚至小于4,可能会发生漏判。

      为了解决上述问题,有必要重新构建阈值模型,即在图 2~图 4中找到一条合适的分界线,能最大程度地分离周跳误判历元和添加周跳历元的|D|值。当均方根较小时,两者较易分离,综合图 2~图 4并考虑不同采样率下模型的通用性和精简性,此处选择线性函数k=6-5×Vsig(0 < Vsig≤0.4) 分离两者,其中, Vsig为均方根值;当均方根位于0.4~0.6之间时,两者较难分离,为了避免周跳漏判,k值设为定值4;当均方根大于0.6时,基本只存在周跳漏判,此时可以适当降低阈值,但为了遵循3σ准则,即k≥3,选择线性函数k=4-2.5×(Vsig-0.6)(0.6 < Vsig≤1) 来避免一定的周跳漏判;当均方根大于1时,k值设为3,此时相位序列波动很大,很可能由较大噪声引起,该种情况出现的几率较小,周跳也较难探测,需联合其他方法进行探测。最终确定的分界线,即阈值模型,如图 2~图 4中的黑线所示,对应的函数模型为:

      (20)

      上述随均方根自适应变化的阈值模型比一些通过采样率和高度角来约束阈值等方法[12]更为准确直接。虽然该模型不能分离所有的周跳误判和周跳历元,但对于单一数据来说能有效减少周跳误判,另外该方法对均方根较大时的周跳漏判也有一定改善,探测性能相比改进之前的方法有明显提高。

    • 在GF组合相邻历元求差法探测周跳时,由于载波观测值的精度很高,在无周跳发生时,GF组合序列的波动范围约在1~2 cm内,GF组合周跳探测的精度很高,通常可达到半周以上,而常见较难探测的整周跳组合为 (1,1),它对GF组合序列的影响为λ1-λ2=5.4 cm,因此,GF阈值一般设为5 cm,以便能探测出较小整周跳。在实际数据周跳探测时,该阈值略为严格,特别是在电离层较活跃或采样间隔较大时,容易造成周跳误判。为了得出合理的阈值模型,图 5~图 7分别给出了无周跳情况下,不同采样率、不同卫星数据的GF组合历元差的时间序列图及对应的卫星高度角变化曲线,其中30 s采样率数据为赤道附近电离层较活跃时的数据。

      图  5  1 s采样率下GPS PRN4

      Figure 5.  GPS PRN4 at 1 s Interval

      图  6  15 s采样率下GPS PRN25

      Figure 6.  GPS PRN25 at 15 s Interval

      图  7  30 s采样率下GPS PRN12

      Figure 7.  GPS PRN12 at 30 s Interval

      图 5~图 7可以看出,GPS卫星相邻历元间的GF组合变化量总体上随着采样间隔的增加而增加,这主要是由于电离层的变化量增大引起的;对于不同采样率的数据,GF组合历元差与卫星高度角呈现高度的相关性。由此可见,周跳探测阈值应随着不同采样率数据和不同卫星高度角适当地变化。由图 5~图 7,当采样间隔较小 (R≤5 s) 时,GF历元差序列波动较小,此时阈值可设为常规的5 cm;当采样间隔较大 (R≥30 s) 时,阈值可以设为15 cm;当采样间隔位于两者之间时,可以仿照文献[12]中的随采样间隔R变化的模型选择线性函数进行过渡。基于上述分析并结合平时大量数据处理经验,得到随采样间隔R变化的经验阈值模型 (单位为m):

      (21)

      对于每一组数据,GF组合历元差与高度角成反比,主要原因是卫星高度角低时,多路径效应和观测噪声的影响较大,因此,需要在上述初始阈值中引入一个与卫星高度角有关的加权因子l[13],加权因子的表达式如下:

      (22)

      式中,e表示当前历元的卫星高度角;E表示临界卫星高度角。只有当e小于E时,加权因子才起作用,对初始阈值T0进行约束。E的取值视实际观测环境而定,通常情况下取15°~30°之间的某个值,本文取30°。综合式 (21) 和式 (22),可以得出GF组合周跳探测的加权阈值模型如下:

      (23)

      上述加权阈值模型是在对大量不同地区、不同采样率、不同卫星高度角数据进行测试的基础上得来的。该模型能适用于不同观测环境下数据的周跳探测,与原来的固定阈值相比有较大的优势。它不仅减少了低卫星高度角时的周跳误判,而且在一定程度上削弱了采样率过低对GF组合周跳探测性能的影响。

    • 本文选取5组不同采样率、不同观测环境下的数据进行实验分析,其中1 s采样率数据一组,为司南K505接收机采集的动态船载数据;15 s和30 s采样率数据各两组,分别为NovAtel ProPark6和Trimble NetR9接收机采集的静态数据。其中为了突出GF组合加权阈值模型的优势,选取的15 s和30 s采样率数据均是电离层较活跃时的观测数据。原始观测数据经MW组合和GF组合检验均无周跳,为了测试构建阈值模型后的MW组合和GF组合的周跳探测性能,每组数据任选一颗GPS卫星,在其L1L2载波相位观测值序列上的等间隔位置加入常见较难探测的周跳,添加的周跳组合如表 1所示。

      表 1  周跳模拟情况

      Table 1.  Situation of Cycle-Slips Simulation

      序号1234567
      L15190002
      L2417-120.52

      图 8表示采样率为1 s、15 s和30 s时的宽巷模糊度变化量与均方根的比值D的变化曲线,其中红线k表示自适应均方根变化阈值模型,蓝线k0表示原阈值±4。三组数据加入周跳的历元间隔分别为800、180和80。由图 8(a)~8(c)中的D值和阈值k的位置关系可以看出,MW组合可以探测出所有L1L2载波上周跳之差为两周的周跳组合;对于周跳之差为1周及以内的周跳组合,该方法仅在某些观测数据质量较好的弧段才能探测出来;另外,MW组合探测不出两个载波上发生的相同周跳。由此可知,MW组合周跳探测的精度约为2周,在观测数据较好时可以达到1周,且其探测精度随采样率的增加而降低。

      图  8  1 s、15 s和30 s采样率下D值变化曲线

      Figure 8.  Change of D Value at 1s, 15 s and 30 s Interval

      对比原阈值模型和自适应均方根变化阈值模型可以看出,图 8(b)中的第19和853历元以及图 8(c)中的第362历元的D值均超出蓝线范围而未超出红线,图 8(a)中的观测数据表现得更为明显,多处历元的D值超出蓝线范围。此时若仍采用原阈值±4,则会发生多处周跳误判,而改进后的自适应均方根变化阈值模型则极大地避免了周跳误判。另一方面,在图 8(a)的第800历元和图 8(b)的第180历元的D值在原阈值范围以内,本来会发生周跳漏判,采用改进后的阈值模型则探测出了此处的周跳。另一方面,图 8(a)中的数据由于采集于湖面,多路径误差较为严重,导致D值曲线波动较大,有几处发生了周跳误判,但加入阈值模型后的周跳探测性能总体上提高了。综上分析,自适应均方根变化阈值模型能够有效减少周跳误判和漏判。

      图 9表示采样率为1 s、15 s和30 s下的GF组合历元差变化曲线。其中,红线表示加权阈值模型;蓝线表示卫星高度角。三组数据加入周跳的历元间隔分别为800、80和70。由图 9(a)可以看出,当采样率为1 s时,GF组合历元差序列的波动较小,仅在高度角很低时才出现2 cm左右的波动,此时GF组合周跳探测的精度较高,达到0.5周左右,只有发生在第1和第3组的近似波长之比的周跳组合未能探测出来;由图 9(b)图 9(c)可以看出,当采样率较低,电离层较活跃时,GF组合历元差序列存在一定的趋势项,且该趋势项与采样率和卫星高度角有较强的相关性。此时若采用固定阈值0.05 m,则在高度角较低时会发生周跳误判,而加权阈值模型考虑了采样率和卫星高度角两个因素的影响,其变化趋势与GF组合历元差序列的变化趋势基本一致,该阈值模型在采样率较低、电离层活跃时可以有效减少周跳误判。但同时也会漏判一些对GF组合序列影响较小的周跳组合,如图 9(c)中的 (1,1) 和 (2,2) 等,另外对于L1L2载波上近似波长之比的周跳组合会探测失效,如图 9(b)图 9(c)中的 (9,7) 和 (5,4),此时可根据GF组合历元差序列适当地调整式 (21) 中的参数或通过事后质量检验来探测剩余的小周跳。因此,在实际周跳探测过程中,应当根据观测环境和条件合理地选择阈值模型。

      图  9  1 s、15 s和30 s采样率下GF组合历元差变化曲线

      Figure 9.  Change of the Difference of GF Combination Between Epochs at 1 s, 15 s and 30 s Interval

    • 本文在滑动平均滤波MW组合周跳探测方法中加入了自适应观测序列均方根变化的阈值模型,该模型通过反映观测数据质量的均方根来约束阈值,从而能够有效减少周跳误判和漏判的发生,提高MW组合周跳探测的性能;另一方面在相邻历元求差的GF组合周跳探测方法中加入了随采样率及卫星高度角变化的加权阈值模型,该模型削弱了采样率和卫星高度角过低对该方法探测性能的影响,显著提升了GF组合周跳探测的可靠性。添加阈值模型后的TurboEdit方法与原来的方法相比不仅在探测精度上有一定的提高,在减少周跳误判和漏判上也有较大改善,并且两者能探测出对方探测失效的周跳盲区,因此,两者组合,可以探测出绝大多数的周跳。

参考文献 (13)

目录

    /

    返回文章
    返回