边缘计算的大规模传感器(65周年施闯)

边缘计算的大规模传感器(65周年施闯)(1)

本文内容来源于《测绘学报》2022年第7期(审图号GS京(2022)0495号)

广域实时精密定位与时间服务系统

施闯1

边缘计算的大规模传感器(65周年施闯)(2)

,辜声峰2,楼益栋2,郑福1,宋伟1,张东1,毛飞宇2

1. 北京航空航天大学卫星导航与移动通信融合技术工业和信息化部重点实验室, 北京 100191;2. 武汉大学卫星导航定位技术研究中心, 湖北 武汉 430079

基金项目:国家自然科学基金(41931075;42174029)

摘要:广域实时精密定位与时间服务已成为GNSS应用领域研究热点, 目前国内外学者围绕其模型算法已展开大量的研究。本文重点论述广域实时精密定位与时间服务数据的处理方法和服务系统, 给出了基于不同基准约束的卫星钟差解算数学模型, 提出通过引入外接原子钟测站、标准时间源(UTC/BDT)等不同时间基准, 构建卫星拟稳基准、外接原子钟跟踪站拟稳基准及标准时间源等约束下的钟差解算模型, 分析了时间基准对精密单点定位和精密单点授时的影响。本文采用实时卫星轨道、钟差、相位偏差、电离层延迟等服务产品及跟踪站实时数据, 验证了系统产品可靠性及终端定位与时间服务性能。实测结果表明: GPS轨道径向精度1.8 cm, 钟差STD精度约0.05 ns; BDS-3轨道径向精度6.7 cm, 钟差STD精度优于0.1 ns; GPS和BDS-2电离层改正精度分别为0.74 TECU与1.03 TECU。基于该产品实现了用户端PPP、PPP-RTK及PPT、PPT-RTK服务, 满足了用户实时厘米级定位和优于0.5 ns的单站时间传递服务, 当采用GPS BDS-2 PPP-RTK解算时, 平面收敛至5 cm约需要12 min。

关键词:广域实时 精密单点定位 精密时间服务 PPP-RTK PPT-RTK

边缘计算的大规模传感器(65周年施闯)(3)

边缘计算的大规模传感器(65周年施闯)(4)

引文格式:施闯, 辜声峰, 楼益栋, 等. 广域实时精密定位与时间服务系统[J]. 测绘学报,2022,51(7):1206-1214. DOI: 10.11947/j.AGCS.2022.20220153

边缘计算的大规模传感器(65周年施闯)(5)

SHI Chuang, GU Shengfeng, LOU Yidong, et al. Real-time wide-area precise positioning and precise timing service system[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(7): 1206-1214. DOI: 10.11947/j.AGCS.2022.20220153

阅读全文:http://xb.chinasmp.com/article/2022/1001-1595/20220711.htm

引 言

得益于全天时、全天候、全球覆盖、高时效、低成本等优势,包括中国的BDS、美国的GPS、俄罗斯的GLONASS及欧洲的Galileo在内的四大全球导航卫星系统(GNSS), 已成为各国定位导航与授时(positioning, navigation, and timing, PNT)的关键基础设施[1]。随着数字信息时代中导航与位置服务等新兴产业的发展,卫星导航系统基本服务性能已难以满足日益增长的精密定位与授时需求。以北斗卫星导航系统为例,其基本定位与授时服务精度分别为10 m与20 ns[2-4],然而自动驾驶、新一代通信技术、精密测控等进一步提出了实时厘米级定位,以及亚纳秒级授时需求。

为实现广域、甚至全球范围内厘米级定位及亚纳秒级授时服务,美国喷气推进实验室(Jet Propulsion Laboratory, JPL)建立了全球差分系统(global differential GPS, GDGPS),并采用RTG软件实现了GPS双频高程20 cm,水平10 cm的实时精密定位。值得注意的是,虽然该系统称为差分系统,但实际上是基于精密单点定位技术(precise point positioning, PPP)实现[5]。随后,国际GNSS服务组织(IGS)于2007年启动了实时试验计划(real-time pilot project, RTPP, https://igs.org/rts/monitoring/),展开广域高精度实时GNSS服务(real-time service, RTS)应用研究,以满足实时地震灾害预警、空间天气监测等应用需求[6-7]。RTS服务采用状态域表达(state space representation, SSR),将GNSS观测误差分解为轨道、钟差等,实现用户端实时PPP处理。文献[8]基于一周数据评估了RTPP组合产品IGS03服务性能,指出实时GPS轨道和钟差精度分别为5、8 cm,GLONASS轨道和钟差精度分别为13、24 cm,同时GPS和GLONASS产品完整率分别高于95%和90%。近年来,北斗和Galileo系统的建设为卫星导航授时服务提供了新的契机,部分IGS实时分析中心开始处理并播发多系统实时轨道、钟差产品。文献[9]以欧洲定轨中心(CODE)事后产品为参考,对法国宇航局(CNES)四系统实时产品进行了比较分析,并通过全球分布多系统跟踪站PPP处理,指出其定位精度随着参与定位的系统数增加而提高。随着非差模糊度固定技术的发展[10-11],文献[12]比较了基于整数钟模型以及未校正相位偏差模型的PPP模糊度固定,并推荐一种兼容当前IGS钟差产品的相位钟模型,以提升PPP定位性能(ftp://igs.gnsswhu.cn/pub/whu/phasebias)。除实时轨道钟差产品,IGS进一步提供了实时全球电离层产品服务(RT-GIM),随着非差非组合模型的发展,高精度电离层增强PPP也成为加快广域实时精密定位用户收敛速度的重要手段之一[13-14]。基于非差模糊度固定及大气增强研究进展,相关学者提出PPP-RTK技术,在局域范围采用非差模型实现了与网络RTK相当的定位效果[15-16]。除IGS提供的RTS服务外,美国Navcom公司的StarFire、荷兰Fugro公司的StarFix/SeaStar、瑞典Hexagon公司的VeriPos等系统相继建成,提供全球实时精密定位服务。然而不论是IGS或商业机构,其实时服务系统主要聚焦于高精度导航与定位服务。

GPS技术在时频传递领域的应用可追溯至文献[17]在20世纪80年代提出的共视法(common view, CV)时间传递。在该方法中两个站通过对同一颗卫星观测,消除卫星端共同的轨道钟差,提高时间传递精度。为提高数据利用率及数据处理可靠性,相关学者提出全视法(all-in-view, AV)时间传递,利用整个导航卫星星座,将全球任意位置测站的时间归算至IGS时间尺度IGST[18]。CV法和AV法都采用伪距观测值实现,考虑相位观测值精度比伪距观测值高2个数量级,相关学者采用相位观测值与IGS事后精密星历产品,通过PPP技术实现了稳定度高于10-15量级的时频传递[19-20]。鉴于PPP技术在时频传递领域的优势,国际计量局(Bureau International des Poids et Mesures, BIPM)提出采用PPP技术取代伪距单点定位技术实现国际原子时(Temps Atomique International, TAI)时频传递[21]。与导航定位中的PPP模型算法发展相对应,国内外学者围绕多系统PPP时间传递[22]、基于RTS的实时PPP时间传递[23]、非差非组合PPP时间传递[24]、非差模糊度固定PPP时间传递[25],以及大气增强时间传递[26]等展开了系统深入的研究。考虑GNSS接收机晶振,尤其原子钟短期稳定度高,通过对接收机钟差建模并将其引入估计模型能显著提升GNSS时间服务精度[27-28]

随着精密单点定位、非差模糊度固定及大气增强技术的发展,以PPP、PPP-RTK算法模型不断完善,为厘米级定位及亚纳秒级授时提供了重要的技术支撑。考虑PPP、PPP-RTK对高精度轨道、钟差、相位偏差等产品的依赖,融合大地测量、卫星导航及互联网技术,自主设计、自主运行的广域实时精密定位与时间服务系统可以进一步提升北斗/GNSS服务性能,拓展北斗/GNSS高水平应用。上述研究主要基于RTS产品展开终端算法研究,且主要局限于导航定位。针对广域实时精密定位和时间服务系统,国内学者也已分别展开有益的探索:文献[29]首先构建了GPS/GLONASS广域实时精密定位原型系统(wide-area precise positioning, WPP)。随着北斗卫星导航系统的发展[30],尤其是北斗精密定轨等技术的成熟[31],基于北斗地基增强系统,有关学者进一步构建了北斗/GNSS广域实时精密定位服务系统[32-33]。文献[34]扩展了广域实时精密定位服务,提出并构建了广域高精度时间服务系统(wide-area precise time, WPT)。针对星基PPP-RTK服务需求,文献[35]设计并采用仿真数据验证了广域星基电离层增强的可行性。然而上述系统时间服务中并未考虑时间溯源需求,且实际服务中仍采用PPP浮点解展开。为此,本文在前期研究的基础上,给出了基于不同基准约束的卫星钟差解算数学模型,通过将基于PPP的时间传递溯源至UTC(k),实现了精密单点授时服务(PPT);同时将模糊度固定及电离层增强引入该实时系统,进一步实现了PPP-RTK与PPT-RTK服务;此外,本文对系统架构以及算法的比较分析,指出了广域实时精密定位与时间服务的联系与区别,在此基础上,评估了本系统实时服务产品精度,以及PPP/PPP-RTK与PPT/PPT-RTK的性能。

1 广域实时精密定位与时间服务

1.1 服务系统

北斗/GNSS广域实时精密定位与时间服务系统(WPP/WPT service system)提供精密定位和精密时间服务,主要由观测数据源、实时精密数据处理平台及用户终端等部分构成,如图 1所示。其中,黑色部分为定位与时间服务公有部分,蓝色部分为时间服务特有部分。

边缘计算的大规模传感器(65周年施闯)(6)

图 1 广域实时精密定位与时间服务系统流程Fig. 1 Diagram of WPP/WPT service system

图选项

外部数据收集IGS全球跟踪站及中国区域跟踪站实时(事后)GNSS观测数据与广播星历,同时从IGS/IERS服务器下载闰秒、天线表、行星表等表文件。广域实时精密数据处理平台解算卫星精密轨道与钟差、信号偏差,以及大气延迟模型等服务产品,在对生成的产品进行质量检测后进行编码播发。用户终端接收精密产品,结合本地GNSS观测数据,采用PPP方法解算获得接收机位置与钟差参数。对于蓝色所示时间服务部分,部分跟踪站外接高精度高稳定原子钟实现时间参考基准引入。另外通过实时引入UTC(k)或北斗时BDT,可进一步实现时间服务溯源。在本文WPP/WPT系统中,服务端通过引入中国计量院时频参考,实现产品的UTC(NIM)溯源。用户终端针对本地晶振长期频率稳定度差、易受外界环境与老化影响,结合实时处理得到的接收机钟差结果,利用精密时钟调控算法将北斗长稳与终端内部晶振的短稳进行最优化组合,实现终端时钟的时频特性统一至服务系统端的时频参考。

系统不同部分通过通信链路实现数据传输,其中分布于全球的GNSS观测站实时数据流以RTCM3.X或接收机原始二进制格式,采用NTRIP协议同步至数据处理平台。数据处理平台生成精密服务产品后采用RTCM SSR/IGS SSR格式编码,通过NTRIP协议播发至定位授时用户。

1.2 数学模型

广域实时精密定位及时间服务关键技术采用非差PPP方法为基础,国内外学者围绕PPP定位数学模型已展开大量卓有成效的研究,然而鲜有文献深入分析PPP定位与时间服务在数学模型方面的异同。为实现多频率多系统观测统一处理,本文以非差非组合数学模型为基础,分别给出服务端卫星钟差估计数学模型和终端PPP解算数学模型,并重点讨论时间基准偏差对定位与时间服务的影响。

非差非组合模型如式(1)所示[35]

边缘计算的大规模传感器(65周年施闯)(7)

(1)

式中,Pr,fsiΦr,fsi分别为以距离为单位,频率f∈{1, 2, 3}上接收机r至卫星si∈{s1,s2, …,sn}的伪距和相位观测值;ρrsi=|rsi-rr|为对应几何距离;tsitr分别为卫星与接收机钟面时与“真实”时间之差;TZ为天顶对流层延迟,αrsi为对应投影函数;IZsi为卫星si天顶电离层延迟,以TECU为单位,

边缘计算的大规模传感器(65周年施闯)(8)

f和高度角e的转换函数;bfsibr,f分别为卫星端与接收机端伪距偏差;Nr,fsi为浮点模糊度参数,λf为其对应波长;εpεΦ分别为包含多路径误差的伪距和相位测量噪声。

目前IGS钟差处理以无电离层组合(ionosphere-free, IF)为基础,结合式(1)及IF组合

边缘计算的大规模传感器(65周年施闯)(9)

边缘计算的大规模传感器(65周年施闯)(10)

(2)

为消除钟差与伪距偏差之间的相关性,式(2)中卫星与接收机钟差参数分别定义为tsi=tsi bIFsitr=tr br, IF。伪距偏差的引入意味着有GNSS技术确定的接收机钟差与真实接收机钟差之间存在一个系统性偏差,因此在授时应用中一般需要进行接收机硬件延迟校准[34]

进一步由于“真实”时间无法获得,需要人为定义时间基准

边缘计算的大规模传感器(65周年施闯)(11)

(3)

式中,t=(tstr)T为卫星钟差与接收机钟差组成的参数向量;Qtt为钟差参数向量协方差向量,其不同取值对应不同的钟差基准

边缘计算的大规模传感器(65周年施闯)(12)

(4)

式中,Ukk维单位阵,此时认为对应的钟稳定度一致。考虑对不同频标建模时,还可根据各时频源稳定度确定矩阵Uk[37]。结合式(4),当没有外接原子钟跟踪站或标准时间源引入时,|t|2Qtt=min对应卫星钟差重心基准;当引入外接原子钟跟踪站时,|t|2Qtt=min对应原子钟跟踪站拟稳基准;当引入标准时间源时,|t|2Qtt=min将系统时间基准溯源至该标准时间。此时卫星钟差不再是“真实”钟差,而受时间基准t0影响

边缘计算的大规模传感器(65周年施闯)(13)

(5)

将卫星钟差式(5)引入式(2)即得到用户端PPP模型,此时接收机钟差为

边缘计算的大规模传感器(65周年施闯)(14)

(6)

由式(6)可以看出,定位用户不受时间基准t0影响,然而时间服务用户却受到t0影响。当用户rArB同时接收该系统时间服务时,其时间同步结果为

边缘计算的大规模传感器(65周年施闯)(15)

(7)

即时间同步不受所选时间基准影响。

2 试验分析

根据上述分析,时间基准的选择并不影响广域实时系统定位及时间同步服务,然而PPP、PPP-RTK单向授时精度却受限于时间基准的稳定性。为验证本文广域实时精密定位与时间服务系统产品精度与终端服务性能,本文采用FUSING软件实现上述时间基准约束算法,并引入模糊度固定及大气增强技术,用于系统数据处理服务。目前FUSING已发展为集多系统实时滤波定轨、实时钟差与信号偏差估计、大气延迟建模与监测等于一体的综合性业务化软件平台[11, 36-37]。本文采用系统2022年1月25日至29日实时输出轨道、钟差产品,再结合实时数据流存储观测数据,通过事后仿实时处理得到相位偏差及电离层延迟产品,并基于上述产品实现定位与钟差解算,分析系统产品精度以及系统服务性能。其中实时轨道产品基于全球IGS实时数据流解算,测站分布如图 2所示;相位偏差及电离层延迟基于中国区域分布部分陆态网与北斗地基增强网实时数据流解算,测站分布如图 3蓝色部分所示;图 3红色所示站点为34个用户定位站;图 3绿色所示站点为1个用户时间服务站,该站时钟类型为外接UTC(NIM),具有稳定的时频参考。

边缘计算的大规模传感器(65周年施闯)(16)

图 2 IGS实时全球跟踪站分布Fig. 2 Global distribution of IGS real-time tracking stations

图选项

边缘计算的大规模传感器(65周年施闯)(17)

图 3 中国区域实时跟踪站分布Fig. 3 Regional distribution of real-time tracking stations in China

图选项

2.1 精密服务产品精度分析

本文以德国地学中心事后多系统精密轨道与钟差产品为参考,评估了北斗、GPS、Galileo和GLONASS四系统卫星实时轨道钟差产品。表 2给出了不同系统轨道产品切向、法向、径向及三维平均RMS值,其中GPS与Galileo卫星轨道精度较优,三维RMS平均值分别为0.05、0.04 m。表 3是不同系统卫星钟差RMS与STD值,GPS与Galileo卫星钟差同样精度最优,STD平均值均为0.05 ns。图 4给出了2022年DOY 25—DOY 29共5 d的多系统卫星钟差时间序列(不同系统纵坐标刻度范围不同)。图 5则是相应时段GPS与北斗二号宽巷与窄巷UPD的时间序列。值得说明的是,相位偏差解算虽采用非差非组合观测值,但考虑宽巷非核准相位延迟(uncalibrated phase delay, UPD)更为稳定,有利于实时播发,且目前绝大部分产品以宽巷窄巷方式表达与改正。因此本文在获得非差非组合模糊度后,将其转为宽巷和窄巷模糊度,并提取UPD产品。

表 1 试验说明Tab. 1 Experiment details

类别服务端产品解算终端定位时间服务
卫星系统BDS-2 BDS-3(仅IGS跟踪站) GPS GLONASS
观测时段2022年DOY 25—DOY 29
采样率/s30
观测量轨道钟差:无电离层组合相位偏差电离层:非差非组合非差非组合
观测噪声/m伪距:0.3;相位:0.003
解算模式轨道钟差:实时解算相位偏差电离层:实时数据流存储,事后仿实时解算实时数据流存储,事后仿实时解算
相位中心IGS_14模型改正
相位缠绕模型改正
潮汐固体潮、极移潮汐、海洋潮汐,IERS协议模型改正
相对论IERS协议模型改正
电离层/服务端模型改正
对流层GPT2模型 过程噪声估计
测站分布轨道钟差:图 2(90)相位偏差电离层:图 3蓝色(200)定位:图 3红色(34)时间服务:图 3绿色(1)

表选项

表 2 2022年DOY 25—DOY 29轨道误差RMS平均值Tab. 2 RMS values of orbit error for DOY 25—DOY 29 of 2022

系统切向法向径向3D
G0.0430.0270.0180.054
C20.5320.1530.1260.568
C30.1220.0490.0670.148
E0.0360.0190.0120.042
R0.0900.0630.0360.115

表选项

表 3 2022年DOY 25—DOY 29钟差RMS和STD平均值Tab. 3 RMS and STD values of clock error for DOY 25—DOY 29 of 2022

系统RMSSTD
G0.4030.054
C20.5460.095
C30.3830.096
E0.2020.052
R1.6640.297

表选项

边缘计算的大规模传感器(65周年施闯)(18)

图 4 2022年DOY 25—DOY 29钟差时间序列Fig. 4 Time series of clock error for DOY 25—DOY 29 of 2022

图选项

边缘计算的大规模传感器(65周年施闯)(19)

图 5 2022年DOY 25—DOY 29 UPD时间序列Fig. 5 Time series of UPD for DOY 25—DOY 29 of 2022

图选项

为评估实时电离层模型精度,通过34个用户站固定坐标,采用后处理模式反算电离层为参考,截止高度角设定为15°。图 6为电离层内插残差RMS值,测站按纬度从高到低排列,可以看出,由于低纬度地区电离层更为活跃,其改正精度低于1TECU。对比GPS和BDS-2,其电离层改正精度分别为0.74、1.03 TECU。

边缘计算的大规模传感器(65周年施闯)(20)

图 6 2022年DOY 25—DOY 29用户站固定坐标反算电离层评估精度Fig. 6 RMS values of ionospheric delay residuals for DOY 25—DOY 29 of 2022

图选项

2.2 用户定位与时间服务性能分析

对图 3红色所示站点,采用PPP和PPP-RTK两种策略进行定位分析,定位每小时重启一次。图 7是68%水平定位误差的时间序列,可以看出平面方向收敛速度PPP-RTK相较于PPP有显著的优势,收敛至5 cm所用时间从半小时缩短至10 min。图 8进一步根据图 3中红色测站离最近基准站(图 3中蓝色测站)距离,采用不同颜色给出了不同站PPP-RTK定位68%收敛时序图。由图 8可以看出,由于本系统中国区域电离层延迟采用广域星基建模实现[35],因此不同跟踪站收敛序列与其离最近基准站之间距离没有明显相关性。

边缘计算的大规模传感器(65周年施闯)(21)

图 7 2022年DOY 25—DOY 29 68%水平定位误差序列Fig. 7 Time series of positioning error in 68% level of all stations for DOY 25—DOY 29 of 2022

图选项

边缘计算的大规模传感器(65周年施闯)(22)

图 8 2022年DOY 25—DOY 29不同跟踪站68%水平定位误差序列Fig. 8 Time series of positioning error in 68% level of individual station for DOY 25—DOY 29 of 2022

图选项

表 4是相应的后半小时定位RMS统计值,无论是平面还是高程上PPP-RTK的定位精度都高于PPP,其中平面甚至提升了71%。

表 4 2022年DOY 25—DOY 29定位RMS值Tab. 4 RMS values of positioning error for DOY 25—DOY 29 of 2022cm

类别PPPPPP-RTK提升/(%)
平面4.701.3471.4
高程5.114.845.3

表选项

考虑单向时间传递服务,由于BJC2站外接UTC(NIM),其接收机钟本身相对稳定,因此可用于PPT、PPT-RTK钟差估值性能评估,测试场景如图 9所示。终端钟差解算单天重启,其接收机钟差平均STD见表 5,可以看出,本文广域实时精密定位与时间服务系统能满足优于0.5 ns单向时间传递。同时GPS单向时间传递精度优于BDS-2,而对比PPT和PPT-RTK结果可以看出,PPT-RTK对接收机钟差稳定性的提升并不显著。

边缘计算的大规模传感器(65周年施闯)(23)

图 9 单向时间传递服务测试场景Fig. 9 Test scenario of one-way time transfer service

图选项

表 5 2022年DOY 25—DOY 29钟差序列STD值Tab. 5 STD values of receiver clock error for DOY 25—DOY 29 of 2022

类别PPTPPT-RTK提升/(%)
GPS0.280.290
BDS-20.450.3815.6

表选项

图 10进一步采用重叠阿伦方差评估了接收机钟差稳定性,对于GPS而言,PPT-RTK相对于PPT估计结果略有提升;对BDS-2而言,PPT-RTK模型的稳定度明显优于PPT模型,其提升约为27.8%。对比GPS和BDS-2,GPS静态精密授时精度要高于BDS-2,两者的万秒稳定度分别达2.14E-14和5.56E-14量级,短期项稳定度GPS比BDS-2高1~2个量级。

边缘计算的大规模传感器(65周年施闯)(24)

图 10 2022年DOY 25—DOY 29 GPS/BDS-2 PPT和PPT-RTK解算的接收机钟差稳定度Fig. 10 Overlapping Allan deviation of receiver clock estimated from PPT and PPT-RTK for DOY 25—DOY 29 of 2022

图选项

3 结论

本文介绍了一种融合GNSS定位与时间服务的广域实时高精度服务系统,重点基于GNSS广域实时精密定位与时间服务数学模型异同,给出了基于不同基准约束的卫星钟差解算数学模型,比较了卫星重心基准、外接原子钟跟踪站拟稳基准及标准时间源等不同时间基准的引入钟差解算法,并分析了其对PPP、PPT的影响。在此基础上结合系统实测轨道、钟差产品,以及实时数据流观测数据,采用事后仿实时方式生成相位偏差与电离层延迟模型等,验证了系统端产品精度及终端服务性能。结果表明,系统能生成包括BDS-3在内的全球4大卫星导航系统轨道、钟差、相位偏差等产品,其中GPS轨道径向精度1.8 cm,钟差STD精度约0.05 ns;BDS-3轨道径向精度6.7 cm,钟差STD精度优于0.1 ns。中国区域覆盖电离层模型在低纬度地区改正精度略低,总体而言,GPS和北斗二电离层改正精度分别为0.74、1.03 TECU。基于该产品,用户可实现实时厘米级精密定位和优于0.5 ns秒的单站时间传递服务,PPP-RTK相比PPP能显著提升定位收敛速度和精度,而PPT-RTK对PPT时间服务精度提升较为有限,稳定度有一定提升。

作者简介

第一作者简介:施闯(1968-), 男, 博士, 教授, 博士生导师, 研究方向为北斗/GNSS高精度定位导航授时及其应用。E-mail: shichuang@buaa.edu.cn

初审:张艳玲

复审:宋启凡

终审:金 君

资讯

,

免责声明:本文仅代表文章作者的个人观点,与本站无关。其原创性、真实性以及文中陈述文字和内容未经本站证实,对本文以及其中全部或者部分内容文字的真实性、完整性和原创性本站不作任何保证或承诺,请读者仅作参考,并自行核实相关内容。文章投诉邮箱:anhduc.ph@yahoo.com

    分享
    投诉
    首页