文章快速检索  
  高级检索
GNSS-IR双频数据融合的土壤湿度反演方法
荆丽丽1, 杨磊1, 汉牟田2, 洪学宝2, 孙波1, 梁勇1     
1. 山东农业大学 信息科学与工程学院, 泰安 271019;
2. 北京航空航天大学 电子信息工程学院, 北京 100083
摘要: 目前全球导航卫星系统反射信号干涉测量(GNSS-IR)土壤湿度反演研究仅针对单一频点展开,提出用熵值法将2个频点数据进行融合以改进土壤湿度反演精度。首先,利用频谱分析法分别解析出各频点的信噪比(SNR)序列的振荡频率,计算出对应的等效天线高度,并利用最小二乘法求解各频点信噪比序列相位;然后,通过熵值法进行2个频点的相位观测量融合;最后,利用融合结果与实测土壤湿度建立经验模型,实现土壤湿度反演。利用地基观测实验获得的全球定位系统(GPS)L1和L2信噪比数据对该方法进行了验证,结果表明:L1和L2双频融合反演结果平均标准差为0.6%,比L1单频反演结果提高64.73%,比L2单频反演结果提高32.12%;均方根误差为0.37%,比L1频点降低72.8%,比L2频点降低73.4%。
关键词: 全球导航卫星系统反射信号干涉测量(GNSS-IR)     土壤湿度反演     双频数据融合     熵值法     全球定位系统(GPS)    
Soil moisture inversion method based on GNSS-IR dual frequency data fusion
JING Lili1, YANG Lei1, HAN Moutian2, HONG Xuebao2, SUN Bo1, LIANG Yong1     
1. College of Information Science and Engineering, Shandong Agricultural University, Tai'an 271019, China;
2. School of Electronic and Information Engineering, Beihang University, Beijing 100083, China
Received: 2018-09-18; Accepted: 2019-01-18; Published online: 2019-01-28 08:31
Foundation item: National Natural Science Foundation of China (41171351); National Key R & D Program of China (2016YFC0803104); National High-tech Research and Development Program of China (2013AA102301); Open Project of National Engineering Research Center for Information Technology in Agriculture (KF2015W003); The Basic Public Welfare Research Project in Zhejiang Province (LGN19D040001); Shandong Agricultural University Top Disciplines Foundation (xxxy201703)
Corresponding author. YANG Lei, E-mail: yanglei_sdau@163.com
Abstract: At present, the study of soil moisture inversion in the field of global navigation satellite signal-interferometer and reflectometry (GNSS-IR) is only for single frequency deployment. In the paper, we propose a method that uses the entropy method to fuse two frequency to improve the accuracy of soil moisture inversion. First, the spectrum analysis method is used to analyze the oscillation frequency of the signal-to-noise ratio (SNR) sequence of each frequency point, and calculate the corresponding equivalent antenna height. The different frequency phase of SNR sequence can be solved by least square method. Then, the phase observation of two frequencies is fused by the entropy method. Finally, an empirical model was established by using the fusion results and the measured soil moisture to achieve soil moisture inversion. The method was verified by global positioning system (GPS) SNR ratio data obtained in frequency L1 and L2 by ground-based observation experiments. The results show that the average standard deviation of the L1 and L2 inversion results after dual frequency fusion is 0.6%, which is 64.73% higher than the L1 frequency inversion results and 32.12% higher than the L2 frequency inversion results. And the RMSE is 0.37%, 72.8% lower than L1 frequency and 73.4% lower than L2 frequency.
Keywords: global navigation satellite signal-interferometer and reflectometry (GNSS-IR)     soil moisture inversion     dual frequency data fusion     entropy method     global positioning system (GPS)    

土壤水分又称土壤湿度或土壤含水量,是联系陆-气相互作用的关键物理量之一[1]。准确测定土壤水分可以有效管理灌溉、保证作物稳产高产等[2-3],因此研究快速、实时、准确、高效、廉价的大面积农田土壤湿度监测方法是精准农业领域一个至关重要的研究方向。全球导航卫星系统反射信号干涉测量(Global Navigation Satellite Signal-Interferometer and Reflectometry,GNSS-IR)是一种新兴的地基遥感技术[4],其利用GNSS直射与地表反射信号的干涉效应,通过获取干涉信号所携带的反射面的相关特性,进行地物参数的遥感,如土壤湿度[5-7]、积雪厚度[8-11]等。

基于GNSS-IR的土壤湿度测量于2008年由Larson等[12]首次提出,该技术使用接收机记录的信噪比(Signal-to-Noise Ratio, SNR)数据进行土壤湿度测量,其物理本质是利用直射信号和反射信号之间的干涉效应,而接收机记录的SNR数据可以看作是对干涉信号的度量,Larson等[12]通过实验证明了SNR的振荡幅度、频率和相位3个观测量与土壤湿度相关。经过近年的发展,GNSS-IR技术已得到国内外专家学者的广泛关注。国外方面,2014年, Alonso-Arroyo等[13]改变接收天线的极化方式,提出通过跟踪垂直极化(V-Pol)和水平极化(H-Pol)SNR数据间的相位差进行土壤湿度测量的方法,该方法提高布鲁斯特角度估计的准确性,从而提高土壤水分探测的准确性,并进行了实验验证。2016年,Roussel等[14]使用参考站接收机,对高度角为2°~70°的全球定位系统(Globa Positioning System, GPS)和GLONASS的SNR数据进行了处理分析,并对从低仰角(2°~30°)和高仰角(30°~70°)SNR数据中提取出的观测量进行了融合,结果表明融合后观测量与土壤湿度间的相关性高于融合前观测量与土壤湿度间的相关性。2017年, Yang等[15]提出干涉信号的解析模型,并证明了从SNR数据中提取介电常数进而反演土壤湿度的可行性,结果表明北斗信号在土壤湿度反演方面效果明显。

国内方面,2015年,敖敏思等[16]结合仿真和实测土壤湿度数据、GPS观测值开展对比试验,发现SNR能有效跟踪土壤湿度变化,最大有效测量范围为45 m,并指出利用指数函数能较好地描述SNR多径延迟相位与土壤湿度之间的关系。同年,徐晓悦等[17]通过Lomb-Scargle变换分析了SNR数据中的多径反射分量,分析结果显示, 幅度与土壤湿度的相关性较强。汉牟田等[18]于2016年根据GNSS的SNR和干涉效应估计方法,推导出利用GNSS干涉信号幅度反演土壤湿度的半经验模型,并对此进行了仿真验证。2016年,严颂华等[19]探讨并通过实验验证了北斗B1波段信号干涉测量土壤水分的可行性。2017年,金双根等[20]对GNSS-R干涉测量的应用和植被覆盖问题进行了综合分析。2018年, 段睿等[21]针对单天线模式下土壤湿度均方根误差较差等问题,提出SVRM辅助的土壤湿度反演方法,并通过实验证明该方法能有效提高土壤湿度均方根误差。

以上研究仅利用单一频点GNSS干涉信号反演土壤湿度,忽视了不同频点信号间的差异性及潜在的更丰富的观测信息。本文通过双频数据融合,综合双频观测信息,提升土壤湿度反演性能。

1 GNSS-IR土壤湿度反演原理

GNSS-IR是一种通过土壤反射的GNSS信号与直射信号的干涉效应实现土壤湿度反演的遥感技术[22]。GNSS接收天线在接收直射信号的同时,因周围环境(海面、土壤表面等)的反射也会接收到反射信号,当天线架设高度较低时,反射信号同直射信号相比具有相同的频率,因此可在天线处产生较为稳定的干涉效应。干涉场景如图 1所示。

图 1 GNSS-R干涉测量原理图 Fig. 1 Schematic diagram of GNSS-R interferometry

这一干涉效应体现在接收机记录的SNR中。SNR可以用直射和反射信号表示为[12]

(1)

式中:AdAm分别为直射、反射信号的幅度;ψ为两信号的相位差;θ为卫星高度角。当GNSS卫星高度角θ发生变化时,AdAmψ均会发生变化,造成SNR的波动。

根据式(1),将二阶多项式拟合法的结果作为直射分量,并将其从SNR序列中分离得到SNR的反射分量。通过频谱分析法得到SNR序列的频率,进而求得等效天线高度,选取值最大的等效天线高度值作为等效天线高度的估计值。

去除直射分量后,SNR的反射分量可表示为

(2)

式中:h为等效天线高度;λ为波长;ϕ为相位。

根据式(2)所述的形式,对SNR的反射分量序列进行最小二乘拟合,得到相位观测量,建立其与实测土壤湿度之间的映射关系,即可用来反演土壤湿度。

2 基于熵值法的双频数据融合

根据香农信息论,信息熵可以量化地衡量信息的不确定性。在GNSS-IR土壤湿度反演中,不同频点SNR序列的相位观测量所含的土壤湿度信息不同,可以根据文献[23]所述的熵值法计算各频点SNR序列相位观测量的信息熵来确定其权值,加权作差后得到融合相位观测量。

具体地,对于频点k(k=1, 2),第i颗卫星的第j天相位观测量为xij(k)(i=1, 2, …, M; j=1, 2, …, N),其在整个观测过程中的相位观测量序列为Xi(k)

首先计算相位观测量序列Xi(k)与实测土壤湿度SM的正负相关性ri(k)

(3)

并根据计算结果选择式(4)求得第i颗卫星第j天相位观测量的指标:

(4)

将第i颗卫星第j天相位观测量代入式(5),共N天的数据,计算得到第i颗卫星第j天相位观测量的比重zij(k)

(5)

分别将第i颗卫星第j天相位观测量的比重zij(k)代入式(6),计算得到第i颗卫星相位观测量的指标熵:

(6)

将第i颗卫星相位观测量的指标熵代入式(7),求得第i颗卫星相位观测量的熵值法冗余度:

(7)

根据第i颗卫星相位观测量熵值法冗余度和所有卫星冗余度和的比值,计算得到第i颗卫星相位观测量的指标熵的权重:

(8)

将第i颗卫星相位观测量的指标熵的权重代入式(9)计算指标评价得分,得到第i颗卫星第j天的相位观测量的指标评价得分:

(9)

按上述方法求得第i颗卫星第j天各个频点相位观测量的指标评价得分后,分别与相应的原始相位观测量相乘并代入式(10)进行数据融合,得到融合后的相位观测量:

(10)

式中:Xij为融合后的相位观测量。

根据式(10)进行2个频点相位观测量的融合,将数据分为训练集与测试集,用训练集建立相位观测量与土壤湿度间的经验模型,并进行土壤湿度反演。

3 实验数据处理及结果分析 3.1 实验概况

本文使用2014-02-04—2014-03-21在法国图卢兹(Toulouse)市的Lamasquere(43°29′14.45″N,1°13′44.11″E)采集的GPS L1和L2 SNR数据。实验场地四周空旷无遮蔽,整个实验期间为裸土,现场如图 2所示。

图 2 GNSS-IR土壤湿度地基实验场地信息 Fig. 2 Information of GNSS-IR soil moisture ground-based experiment site

实验使用Leica GR25接收机和AR10基准站天线,天线指向天顶,天线架设高度为1.70 m,采样频率为1 Hz。距离天线相位中心在地表投影约2 m的位置布置2枚ML3 Theta Probe土壤湿度传感器采集原位土壤湿度数据,精度为±1%,深度分别为2 cm和5 cm,采样间隔为10 min。

3.2 实验数据处理

数据处理流程如图 3所示。

图 3 数据处理过程 Fig. 3 Flow of data processing

具体过程如下:①获取实验数据。从采集的数据中提取SNR、时间、卫星高度角、方位角等信息,将L1和L2频点的SNR定义为SNR(1)、SNR(2)。②数据预处理。截取卫星高度角在2°~30°范围的SNR数据,同时删除SNR为零的数据。③数据分离。将数据分为训练集和测试集。④去除趋势项。对SNR(1)、SNR(2)的卫星高度角正弦值进行二阶多项式拟合,并以此作为直射分量的近似值,去除直射分量,得到反射分量SNRm(1)、SNRm(2)。⑤频谱分析。利用Lomb-Scargle方法对反射分量进行频谱分析,求得频率f(1)f(2),进而求得等效天线高度h(1)h(2)。⑥求相位观测量。根据式(2)对SNRm进行最小二乘拟合求得相位观测量P(1)P(2)。⑦数据融合。根据式(10)进行融合处理,得到融合后的相位观测量P。⑧建模并验证。选取融合前后的相位观测量分别与对应的土壤湿度观测值进行模型建立,将测试数据代入模型进行验证。

3.3 实验结果分析

根据3.2节处理过程,将L1和L2的SNR序列相位观测量进行融合,分别建立L1、L2及融合的相位观测量与同比土壤湿度间的一元线性回归模型,结果如图 4~图 7所示。

图 4 PRN4建模结果 Fig. 4 Modeling results of PRN4
图 5 PRN7建模结果 Fig. 5 Modeling results of PRN7
图 6 PRN12建模结果 Fig. 6 Modeling results of PRN12
图 7 PRN24建模结果 Fig. 7 Modeling results of PRN24

图 4~图 7分别为PRN4、PRN7、PRN12及PRN24的建模结果,其中图 4(a)图 5(a)图 6(a)图 7(a)为融合后的相位观测量与同比土壤湿度的建模结果,图 4(b)图 5(b)图 6(b)图 7(b)为单频点(以L1频点为例)的相位观测量与同比土壤湿度的建模结果。从图中可以直观地看出融合后的相位观测量与土壤湿度间的相关性更高。

模型建立完成后,将测试数据代入模型得到反演结果。同时,对单频SNR序列做同样处理得到单频反演结果。图 8给出了PRN4、PRN7、PRN12、PRN24 4颗卫星的处理结果。

图 8 不同卫星反演结果对比 Fig. 8 Comparison of different satellite inversion results

图 8为各颗卫星反演结果对比,其中融合后PRN4的标准差为0.68 %,比L1提高42.93%,比L2提高43.23%;均方根误差(RMSE)为0.49%,比L1降低81.58%,比L2降低77.52%;PRN7的标准差为0.62 %,比L1提高43.85%,比L2提高47.01%;RMSE为0.98%,比L1降低58.82%,比L2降低58.47%;PRN12的标准差为0.81%,比L1提高51.30%,比L2提高41.03%;RMSE为1.28%,比L1降低13.51%,比L2降低46.67%;PRN24的标准差为0.92 %,比L1提高3.14%,比L2降低3.84%;RMSE为0.95%,比L1降低64.15%。

对所有卫星同一时间段(13:00—15:00)的反演结果取平均,结果如图 9所示。

图 9 反演结果平均对比分析 Fig. 9 Inversion results average comparison analysis

图 9(a)为融合前L1频点、L2频点和融合后反演结果以及FDR采集的同比数据随时间变化的对比结果,从对比中可以清晰显示融合后的反演结果与FDR的匹配度更高。图 9(b)为反演结果线性对比,图中计算了融合前两频点的反演结果、融合后的反演结果与同比数据间的决定系数RL12RL22R2及均方根误差RMSEL1、RMSEL2、RMSE。通过计算得到如下结果:融合后的决定系数R2为0.97,RMSE为0.37%,R2比L1频点结果提高73.2%,比L2频点提高38.5%;RMSE比L1频点降低72.8%,比L2频点降低73.4%。

通过统计所有卫星的标准差,得到平均结果如下:L1和L2双频融合反演结果平均标准差为0.6%, 精度比L1单频反演结果提高64.73%, 比L2单频反演结果提高32.12%。

4 结论

本文提出了一种基于GNSS-IR双频数据融合的土壤湿度反演方法, 该方法以不同频点SNR序列相位观测量所对应的指标评价得分作为融合的权重标准,通过熵值法对2个频点SNR相位观测量进行融合处理。利用地基观测获取的实验数据对该方法的性能进行了验证,结果表明该方法可显著提高土壤湿度反演精度。从处理结果中可得到如下结论:

1) 不同频点间所含的土壤信息存在差异。熵值法不仅能够克服主观确定权值的随机性和臆断性,而且作为多指标评价方法能够避免不同频点间信息的重复性,综合考虑信息间的差异性,两频点中质量好的数据能得以保留,从而使反演结果得到提高。

2) 综合反演结果分析,熵值法在双频数据融合处理中能够有效提高反演结果的准确性,不仅反演精度有所提高,融合后数据的变化趋势也更加趋近于同比数据。最终结果也表明融合后的反演结果对比单一频点的反演结果有显著提高,且融合后的平均标准差为0.6%。

致谢   感谢北京航空航天大学杨东凯教授在学术上给予的指导和帮助,感谢F. Baup和K. Boniface收集气象数据,以及Roussel.N和F. Frappart收集GNSS观测数据,感谢J. Darrozes为本文提供的数据。

参考文献
[1]
AUBER J C, BIBAUT A, RIGAL J M.Characterization of multipath on land and sea at GPS frequencies[C]//Proceedings of the 7th International Technical Meeting of the Satellite Division of the Institute of Navigation(ION GPS 1994), 1994: 1155-1171.
[2]
GARRISON J L, KATZBERG S, HILL M. Effect of sea roughness on bistatically scattered range coded singals from the global positioning system[J]. Geophysical Research Letters, 1998, 25(13): 2257-2260. DOI:10.1029/98GL51615
[3]
KOMJATHY A, ZAVOROTNY V, AXELRAD P, et al.GPS signal scattering from sea surface: Comparison between experimental data and theoretical model[C]//Proceeding of the 5th International Conference on Marine and Coastal Environments, 1998: 1-10.
[4]
丁金才. GPS气象学及其应用[M]. 北京: 气象出版社, 2009.
DING J C. GPS meteorology and its application[M]. Beijing: Meteorological Publishing Press, 2009. (in Chinese)
[5]
RODRIGUEZ-ALVAREZ N, BOSCH-LLUIS X, CAMPS A, et al. Review of crop growth and soil moisture monitoring from a ground-based instrument implementing the interference pattern GNSS-R technique[J]. Radio Science, 2016, 46(6): 1-11.
[6]
CAMPS A, PARK H, PABLOS M, et al.Soil moisture and vegetation impact in GNSS-R TechDemoSat-1 observations[C]//36th IEEE International Geoscience and Remote Sensing Symposium.NJ: IEEE Press, 2016: 1982-1984. https://www.researchgate.net/publication/309773062_Soil_moisture_and_vegetation_impact_in_GNSS-R_TechDemosat-1_observations
[7]
PALOSCIA S, MACELLONI G, SANTI E, et al. A multifrequency algorithm for the retrieval of soil moisture on a large scale using microwave data from SMMR and SSM/I satellites[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(8): 1655-1661. DOI:10.1109/36.942543
[8]
NIEVINSKI F G.Forward and inverse modeling of GPS multipath for snow monitoring[D].Boulder: University of Colorado Boulder, 2013. https://www.researchgate.net/publication/258848060_Forward_and_Inverse_Modeling_of_GPS_Multipath_for_Snow_Monitoring
[9]
RODRIGUEZ-ALVAREZ N, AGUASCA A, VALENCIA E, et al. Snow thickness monitoring using GNSS measurements[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(6): 1109-1113. DOI:10.1109/LGRS.2012.2190379
[10]
LARSON K M, NIEVINSKI F G. GPS snow sensing:Results from the earth scope plate boundary observatory[J]. GPS Solutions, 2013, 17(1): 41-52. DOI:10.1007/s10291-012-0259-7
[11]
BINBIN L I, ZHANG Y, YANG S, et al. Snow depth altimetry using GNSS signal with single antenna[J]. GNSS World of China, 2016, 41(6): 37-41.
[12]
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
[13]
ALONSO-ARROYO A, CAMPS A, AUGUASCA A, et al. Improving the accuracy of soil moisture retrievals using the phase difference of the dual-polarization GNSS-R interference patterns[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(12): 2090-2094. DOI:10.1109/LGRS.2014.2320052
[14]
ROUSSEL N, FRAPPART F, RAMILLIEN G, et al. Detection of soil moisture variations using GPS and GLONASS SNR data for elevation angles ranging from 2° to 70°[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2016, 9(10): 4781-4794. DOI:10.1109/JSTARS.2016.2537847
[15]
YANG T, WAN W, CHEN X W, et al. Using BDS SNR observations to measure near-surface soil moisture fluctuations:Results from low vegetated surface[J]. IEEE Geoscience and Remote Sensing Letters, 2017, 14(8): 1308-1312. DOI:10.1109/LGRS.2017.2710083
[16]
敖敏思, 朱建军, 胡友健, 等. 利用SNR观测值进行GPS土壤湿度监测[J]. 武汉大学学报(信息科学版), 2015, 40(1): 117-120.
AO M S, ZHU J J, HU Y J, et al. Comparative experiments on soil moisture monitoring with GPS SNR observations[J]. Geomatics and Information Science of Wuhan University, 2015, 40(1): 117-120. (in Chinese)
[17]
徐晓悦, 郑南山, 谭兴龙. GPS-R技术辅助的土壤水含量变化监测[J]. 测绘科学技术学报, 2015, 32(5): 465-468.
XU X Y, ZHENG N S, TAN X L. Monitoring of soil moisture fluctuation in mining areas based on GPS-R[J]. Journal of Geomatics Science and Technology, 2015, 32(5): 465-468. DOI:10.3969/j.issn.1673-6338.2015.05.006 (in Chinese)
[18]
汉牟田, 张波, 杨东凯, 等. 利用GNSS干涉信号振荡幅度反演土壤湿度[J]. 测绘学报, 2016, 45(11): 1293-1300.
HAN M T, ZHANG B, YANG D K, et al. Soil moisture retrieval utilizing GNSS interference signal amplitude[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(11): 1293-1300. DOI:10.11947/j.AGCS.2016.20160145 (in Chinese)
[19]
YAN S H, ZHAO F, CHEN N C, et al. Soil moisture estimation based on BeiDou B1 interference signal analysis[J]. Science China Earth Sciences, 2016, 59(12): 2427-2440. DOI:10.1007/s11430-015-0013-7
[20]
金双根, 张勤耘, 钱晓东. 全球导航卫星系统反射测量(GNSS+R)最新进展与应用前景[J]. 测绘学报, 2017, 46(10): 1389-1398.
JIN S G, ZHANG Q Y, QIAN X D. New progress and application prospects of global navigation satellite system reflectometry(GNSS+R)[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10): 1389-1398. DOI:10.11947/j.AGCS.2017.20170282 (in Chinese)
[21]
段睿, 张波, 汉牟田, 等. SVRM方法的单天线GNSS-R土壤湿度反演[J]. 导航定位学报, 2018, 6(1): 34-39.
DUAN R, ZHANG B, HAN M T, et al. Inversion of soil moisture with singal-antenna GNSS-R signal using SVRM[J]. Journal of Navigation and Positioning, 2018, 6(1): 34-39. (in Chinese)
[22]
杨磊.GNSS-R农田土壤湿度反演方法研究[D].泰安: 泰安农业大学, 2017.
YANG L.Study of GNSS-R cropland soil moisture retrieval method[D].Tai'an: Shandong Agricultural University, 2017(in Chiense). http://www.cqvip.com/QK/90069X/201801/674427807.html
[23]
王富喜, 毛爱华, 李赫龙, 等. 基于熵值法的山东省城镇化质量测度及空间差异分析[J]. 地理科学, 2013, 33(11): 1323-1329.
WANG F X, MAO A H, LI H L, et al. Quality measurement and regional difference of urbanization in shandong province based on the entropy method[J]. Scientia Geographica Sinica, 2013, 33(11): 1323-1329. (in Chinese)
http://dx.doi.org/10.13700/j.bh.1001-5965.2018.0555
北京航空航天大学主办。
0

文章信息

荆丽丽, 杨磊, 汉牟田, 洪学宝, 孙波, 梁勇
JING Lili, YANG Lei, HAN Moutian, HONG Xuebao, SUN Bo, LIANG Yong
GNSS-IR双频数据融合的土壤湿度反演方法
Soil moisture inversion method based on GNSS-IR dual frequency data fusion
北京航空航天大学学报, 2019, 45(6): 1248-1255
Journal of Beijing University of Aeronautics and Astronsutics, 2019, 45(6): 1248-1255
http://dx.doi.org/10.13700/j.bh.1001-5965.2018.0555

文章历史

收稿日期: 2018-09-18
录用日期: 2019-01-18
网络出版时间: 2019-01-28 08:31

相关文章

工作空间