2018, Vol. 38
1 mm的天顶对流层延迟(ZTD)可导致单点定位产生0.1 mm水平误差和2~6 mm垂直误差[1]，是GNSS高精度导航定位中的主要误差源之一。当前，IGS(international GNSS service)分析中心及数据中心为用户提供了不同时间分辨率的ZTD产品(JPL、CDDIS提供5 min产品，NGS、EUR、EMR、GFZ提供1 h产品，CODE、ESA提供2 h产品)，其精度约为2~5 mm[2]。很多学者利用上述ZTD产品进行对流层延迟的时间特征研究[3-4]，建立了适用于区域和全球的ZTD模型[5-7]，并试图建立高时间分辨率的全球ZTD模型[8]。同时，ZTD产品及所建区域或全球模型在其他应用中也获得良好效果。例如，利用全球ZTD模型作为先验约束显著地提高了精密单点定位的收敛速度和定位精度[10]，利用高时间分辨率ZTD产品(时间分辨率5 min)进行大气可降水量反演[11]，利用高精度ZTD模型或数据改正SAR影像大气误差[12]等。
综上可知，ZTD产品及基于此所建模型的时间分辨率不断提高。但是，ZTD产品本身作为一类数据资源，高时间分辨率在满足实时性的同时也会增加数据存储、传输负担并降低数据处理响应速度。当ZTD产品建模精度与产品标称精度水平基本一致时，若进一步提高时间分辨率则弊(降低数据处理响应速度)大于利(增强应用效果)。针对目前ZTD产品发布时间分辨率缺乏确定依据的问题，本文兼顾产品数据存储量大小、ZTD时变特征(建模精度)与产品标称精度对产品发布时间分辨率进行规范化研究。首先，利用历元差分法和区间拟合法分析不同时间分辨率的历元差分和区间拟合误差的统计特征；其次，利用Lomb-Scargle变换[13](简称L-S变换)分析ZTD产品(残余延迟)的短周期性特征，并利用奇异值分解法[14]研究ZTD产品特征；最后，以2014全年12个IGS站点ZTD产品(5 min时间间隔)进行综合研究，分析并确定了ZTD产品发布的规范时间分辨率。
1 ZTD时序分析方法 1.1 区间拟合在区间拟合之前，对不同时间分辨率的对流层产品进行一阶历元差分，通过历元差分ZTD的直方图统计大致确定误差特征突变的关键时间分辨率。进而，对各关键时间分辨率的ZTD时序产品进行区间拟合误差分析，则能进一步研究规范或合理的时间分辨率。
本文利用移动多项式模型对各关键时间分辨率的ZTD时间序列进行区间拟合，并以区间拟合误差进行ZTD时间特征分析。其中，移动拟合区间时长为m=k+2(多项式阶数取k=3、k=2、k=1)。根据区间拟合模型参数(加密)计算5 min时间间隔的ZTD时序，与发布的5 min原始ZTD产品时序同步作差获得区间拟合误差，进而对拟合误差进行概率统计分析。
1.2 L-S频谱分析研究表明，ZTD产品时序中含有明显的年周期和半年周期的趋势项，对原始ZTD时序的频谱分析应先去除趋势项的影响。带有半年和年周期项的拟合模型为：
| $ \begin{align} &y({{t}_{i}})=a+\sum\limits_{\text{ }j=1}^{2}{{{b}_{j}}\text{sin}\ \left( \frac{2\text{ }\!\!\pi\!\!\text{ }j}{365.25}{{t}_{i}} \right)+} \\ &\ \ \ \ \ \ \ \ \sum\limits_{j=1}^{2}{{{c}_{j}}\text{cos}\left( \frac{2\text{ }\!\!\pi\!\!\text{ }j}{365.25}{{t}_{i}} \right)~+x({{t}_{i}})} \\ \end{align} $ | (1) |
式中，y(ti)为ti时刻的原始ZTD，a为常系数，b、c为周期项系数，x(ti)为去除趋势项的ZTD残余延迟。
实际的ZTD时序由于存在数据缺失，往往形成非均匀时间序列。为有效消除因数据缺失而导致的虚假频率，本文采用Lomb-Scargle变换对ZTD残余延迟x(t)进行频谱分析。L-S变换的归一化功率谱Px(f)计算式为[13]：
| $ \begin{align} &{{P}_{x}}\left( f \right)=\frac{1}{2{{\sigma }^{2}}}\left\{ \frac{{{\left[\sum{\Delta {{x}_{t}}\text{cos}2\text{ }\!\!\pi\!\!\text{ }f\left( t-\tau \right)~} \right]}^{2}}}{\sum{\text{co}{{\text{s}}^{2}}2\text{ }\!\!\pi\!\!\text{ }f\left( t-\tau \right)}} \right.+\text{ } \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. \frac{{{\left[\sum{\Delta {{x}_{t}}\text{sin}2\text{ }\!\!\pi\!\!\text{ }f\left( t-\tau \right)~} \right]}^{2}}}{\sum{\text{si}{{\text{n}}^{2}}2\text{ }\!\!\pi\!\!\text{ }f\left( t-\tau \right)}} \right\} \\ \end{align} $ | (2) |
式中，x、Δxt(=x(t)－x)分别为残余延迟均值、零均值残余延迟，f、τ分别表示频率和相位。
1.3 奇异值分解设IGS站点ZTD时序构成的数据矩阵为A(m×n，m个观测历元，n个IGS站点)，且满足rank(A)=n，则存在m阶正交矩阵U和n阶正交矩阵V，使得下式成立[14]：
| $ \mathit{\boldsymbol{A}}=\mathit{\boldsymbol{US}}{{\mathit{\boldsymbol{V}}}^{\text{T}}} $ | (3) |
式中，U和V称为核矩阵，且均为正交矩阵；S=diag(λ1, λ2, …, λr, 0)，λi(i=1, 2, …, r)为矩阵A的全部非零奇异值，且满足λ1≥λ2≥…≥λr > 0。
非零奇异值λi可反映数据矩阵A的分布特征，且奇异值越大，与之对应的行(列)包含的特征越多。SVD分解多用来进行特征提取[14]和模式识别[15]，本文引入SVD分解进行ZTD产品时间分辨率的研究。计算各时间分辨率与原始分辨率的ZTD数据矩阵奇异值之比(λr(i)/λr(0)×100%)，通过比值突变分析确定ZTD产品发布的规范或合理的时间分辨率。
2 ZTD产品时间分辨率的确定 2.1 ZTD数据源目前IGS分析中心提供各基准站自2007年起至今时间分辨率为5 min的ZTD产品，其精度约为2~5 mm。本文顾及ZTD纬度、高程差异，在全球范围内选取了12个IGS站点，并将这12个IGS站点2014年全年的ZTD产品进行稀疏化，得到时间分辨率为15 min、30 min、45 min、60 min、90 min、120 min、150 min、180 min的ZTD时序产品，将5 min原始产品和上述数据作为本文实验数据进行ZTD产品时间分辨率的研究。其中，IGS站点的纬度、高程以及2014年收集的ZTD产品实际天数整理如表 1所示。
|
表 1 IGS站点与ZTD信息 Tab. 1 Information of IGS stations and zenith tropospheric delay |
由§2.1可知，IGS分析中心提供标称精度2~5 mm、时间分辨率高达5 min的ZTD产品，用户可根据需要自行下载使用。但在实际应用中，需进一步基于ZTD产品采用插值算法获得观测历元的天顶对流层延迟。本节利用§2.1全部数据对不同时间分辨率的ZTD产品进行一阶历元差分计算，考察ZTD差分绝对值与标称精度的关系，从而研究不同时间分辨率的ZTD历元差绝对值的分组概率特征。需要说明的是，ZTD数据产品以d为单位，在单天内按各时间分辨率进行一阶差分计算。当缺失数据长度超过分辨率对应的时间间隔时，则放弃一阶差分计算。限于篇幅，图 1仅显示了ABMF站点的统计结果。
|
图 1 ZTD历元差分值概率统计 Fig. 1 The probability histogram of the epoch-differenced ZTD |
从图 1可以发现，对采用的9种时间分辨率，ZTD时序的历元差绝对值均存在着“小集中、大分散”的统计特征。而且，当时间间隔在30 min以内时，ZTD历元差绝对值的主要分布范围随时间间隔的增加而显著增大，且最大分组概率与时间间隔呈明显的负相关；当时间间隔大于30 min时，历元差绝对值主要分布范围和最大分组概率趋于稳定。此外，考察历元差值与ZTD标称精度的关系可知，5 min间隔的历元差绝对值与ZTD标称精度(按平均3 mm误差分析)基本处于相同水平，表明5 min的ZTD产品相邻离散点值变化可以忽略不计。当时间分辨率降低至15 min和30 min，历元差绝对值主要分布范围显著增大，但90%以上不超过10 mm(约3倍标称误差)，表明采用线性模型拟合基本可达ZTD产品标称级精度；当分辨率进一步降低至45 min以上时，10 mm以上误差比例显著增加，表明需要采用高次多项式模型才能提高插值精度。分析其原因，当时间间隔从5 min逐渐增加时，相邻历元ZTD时间相关性及其变化梯度均较高，从而导致分布范围的迅速增大和最大分组概率的迅速减小。当时间间隔增大到45 min以后，相邻历元ZTD时间相关性及其变化梯度均较低，故ZTD历元差主要分布范围和最大分组概率稳定。因此，通过9种时间分辨率的分组概率统计结果与分析，初步判定ZTD产品发布的时间分辨率为30 min，该时间间隔内的ZTD建模精度可达标称精度水平。同时，选定5 min、15 min、30 min、45 min、60 min、90 min、120 min共7种为后续研究的关键时间分辨率。
2.3 区间拟合误差分析利用§1.1区间拟合方法对§2.2中7种关键时间分辨率的ZTD时序进行区间拟合，获得ZTD拟合误差并进行概率统计分析。详细方案为：1)拟合策略分别选用区间长度为3的一阶线性模型(记为方案A)、区间长度为4的二阶多项式模型(记为方案B)和区间长度为5的三阶多项式模型(记为方案C)，对7种关键时间分辨率的ZTD数据进行区间拟合，根据拟合模型重建时间间隔5 min的ZTD时序。2)以IGS分析中心的时间分辨率5 min的ZTD产品作为标准值，与区间拟合重建的ZTD时序作差并统计不同关键时间分辨率下的拟合误差，分组概率统计结果如表 2(单位%)所示。
|
|
表 2 ZTD区间拟合误差概率统计 Tab. 2 Probability statistics of the interval fitting errors |
从表 2方案A可看出，各关键时间分辨率下的拟合误差分布存在显著差别。各区间拟合误差小于ZTD产品标称误差3 mm的比例随时间分辨率降低而降低，分别为99.84%、99.23%、92.82%、81.20%、70.48%、55.27%、46.78%。由此可知，时间间隔30 min是拟合误差分布的时间分辨率拐点。换言之，时间间隔不超过30 min时，拟合误差绝大部分小于3 mm(90%以上)；时间间隔超过30 min时，拟合误差小于3 mm的比例急剧下降且分布范围更广，反映了低时间分辨率ZTD产品对拟合模型的复杂性有更高要求。
综合方案A、方案B和方案C可得，当时间间隔在30 min以内时，相邻历元之间的时间相关性较强，拟合误差小于ZTD产品标称误差3 mm的比例分别达到92.82%(线性拟合)、95.17%(二阶多项式)、96.41%(三阶多项式)。当时间间隔大于30 min时，相邻历元的时间相关性明显降低，模型拟合误差受时间间隔的影响较大，导致拟合误差超过3 mm的比例显著增大。因此，ZTD产品发布的合理时间分辨率可基本确定为30 min。此外，对于所采用的7种关键时间分辨率数据，方案B和方案C的拟合误差概率统计较为一致，略优于方案A(最大差值约为4%)，表明采用二阶多项式是30 min时间间隔ZTD产品的合理建模策略。
2.4 残余延迟频谱分析为检验天顶对流层延迟的高频短期周期性，本文基于§2.1原始ZTD产品(5 min时间分辨率)利用§1.2方法获得ZTD残余延迟并进行L-S频谱分析。具体步骤为：1)根据式(1)对原始ZTD时序进行周期趋势项拟合，并去除趋势项得到残余延迟x(t)；2)利用L-S法对残余延迟进行频谱分析，获得功率谱并分析ZTD产品的高频短周期性特征。为节省篇幅，图 2~3仅显示了ABMF和AJAC站点的原始ZTD及残余延迟，相应站点的残余延迟频谱分析见图 4。
|
图 2 原始ZTD及残余延迟时序(ABMF站) Fig. 2 Original and residual ZTD of ABMF |
|
图 3 原始ZTD及残余延迟时序(AJAC站) Fig. 3 Original and residual ZTD of AJAC |
|
图 4 ABMF站和AJAC站残余延迟频谱分析结果 Fig. 4 Spectrum analysis for the ZTD of ABMF and ABMF and AJAC |
从图 2(a)、图 3(a)可以看出，ZTD数据具有较为明显的长周期项(半年、年)和较为复杂的高频短周期项，若直接对原始ZTD时序进行频谱分析会导致错误的频率结果。图 2(b)、图 3(b)为去除长周期趋势项的残余延迟序列，可以发现，残余延迟序列中含有大量的高频短周期项。图 4显示了ZTD残余延迟序列的频谱分析结果。可以看出，ABMF站与AJAC站的ZTD时序中含有相同的短周期特征，其周期项约为12 h，这与已有文献[4, 7]关于ZTD周期特征的结果基本一致。此外，两个IGS站均含有更高频率的短周期项，如ABMF含约3.5 h的短周期项、AJAC含约2.5 h的短周期项。根据纳斯奎斯采样定理，如果重采样数据完整保留ZTD的全部特征，则ZTD产品发布的时间间隔必须在1 h范围以内。
2.5 奇异值分解为进一步研究ZTD时序的时间特征，本节利用§1.3奇异值分解方法对12个站点不同分辨率的ZTD时序进行特征提取，从而讨论分析ZTD产品发布的规范时间分辨率。具体的实验步骤为：1)选用§2.1中12个IGS站年积日85(春季)和年积日263(秋季)全天的ZTD时序分别构成数据矩阵(§2.2中7种关键时间分辨率)；2)对上述12个IGS站各关键时间分辨率的ZTD数据矩阵进行SVD分解，并统计各数据矩阵的奇异值之比，结果如表 3(单位%)所示。
|
|
表 3 不同分辨率ZTD数据矩阵奇异值比 Tab. 3 Singular value ratio of ZTD matrix with time different resolutions |
从表 3年积日85结果可知，不同时间采样率的数据矩阵所得的奇异值比有较大差别:15 min时间间隔的奇异值比均约为57%, 30 min时间间隔的奇异值比均约为40%，45 min时间间隔及以上的奇异值比显著下降而且12个奇异值比离散程度增加。其中，45 min时间间隔的最大最小奇异值比之差约7.5%，60 min时间间隔的最大最小奇异值比之差约8.7%，90 min时间间隔的最大最小奇异值比之差约14.6%，120 min时间间隔的最大最小奇异值比之差约13.2%。由此可判断，时间间隔30 min是数据矩阵奇异值比特征的时间拐点。换言之，时间间隔30 min的ZTD数据矩阵可较好地保留ZTD产品的时间特征，建议IGS分析中心将ZTD产品发布的时间分辨率规范为30 min。从年积日263结果可得相同结论，在此不再赘述。
3 结语1) 分析天顶对流层产品发布的时间分辨率现状，综合产品数据存储量、建模和标称精度水平等，指出天顶对流层产品发布时间分辨率规范化的现实意义，进而探讨基于历元间差分、区间拟合、L-S频谱分析和奇异值分解算法的ZTD时序特征分析方法。
2) 利用12个IGS站点2014年天顶对流层产品，综合原始5 min和重采样15 min、30 min、45 min、60 min等时间分辨率ZTD产品进行历元差分概率统计、区间拟合误差分析、L-S频谱分析和奇异值分解研究，各方法研究结果及分析结论较为一致，即ZTD产品发布的合理或规范时间分辨率为30 min。
| [1] |
Li W, Yuan Y, Ou J, et al. New Versions of the BDS/GNSS Zenith Tropospheric Delay Model IGGtrop[J]. Journal of Geodesy, 2015, 89(1): 73-80 DOI:10.1007/s00190-014-0761-5
(  0) |
| [2] |
Byun S H, Bar-Sever Y E. A New Type of Troposphere Zenith Path Delay Product of the International GNSS Service[J]. Journal of Geodesy, 2009, 83(3): 1-7
( 0) |
| [3] |
吕慧珠, 黄文德, 闻德保. 一种基于频谱分析和AR补偿的对流层延迟预报模型[J]. 大地测量与地球动力学, 2015, 35(2): 283-286 (Lü Huizhu, Huang Wende, Wen Debao. A Tropospheric Delay Prediction Model Based on Spectrum Analysis and the AR Compensation[J]. Journal of Geodesy and Geodynamics, 2015, 35(2): 283-286)
( 0) |
| [4] |
Jin S, Luo O F, Gleason S. Characterization of Diurnal Cycles in ZTD from a Decade of Global GPS Observations[J]. Journal of Geodesy, 2009, 83(6): 537-545 DOI:10.1007/s00190-008-0264-3
( 0) |
| [5] |
Takeichi N, Sakai T, Fukushima S, et al. TroposphericDelay Correction with Dense GPS Network in L1-SAIF Augmentation[J]. GPS Solutions, 2010, 14(2): 185-192 DOI:10.1007/s10291-009-0133-4
( 0) |
| [6] |
毛健, 朱长青, 郭继发. 一种新的全球天顶对流层延迟模型[J]. 武汉大学学报:信息科学版, 2013, 38(6): 684-688 (Mao Jian, Zhu Changqing, Guo Jifa. A New Global Zenith Tropospheric Delay Model[J]. Geomatics and Information Science of Wuhan University, 2013, 38(6): 684-688)
( 0) |
| [7] |
姚宜斌, 胡羽丰, 张豹. 利用多源数据构建全球天顶对流层延迟模型[J]. 科学通报, 2016(24): 2 730-2 741 (Yao Yibin, Hu Yufeng, Zhang Bao. Establishment of a Global Zenith Tropospheric Delay Model Using Multi-Source Data[J]. Chinese Science Bulletin, 2016(24): 2 730-2 741)
( 0) |
| [8] |
姚宜斌, 胡羽丰, 余琛. 一种改进的全球对流层天顶延迟模型[J]. 测绘学报, 2015, 44(3): 242-249 (Yao Yibin, Hu Yufeng, Yu Chen. An Improved Global Zenith Tropospheric Delay Model[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(3): 242-249 DOI:10.11947/j.AGCS.2015.20140089)
( 0) |
| [9] |
姚宜斌, 余琛, 胡羽丰, 等. 利用非气象参数对流层延迟估计模型加速PPP收敛[J]. 武汉大学学报:信息科学版, 2015, 40(2): 188-192 (Yao Yibin, Yu Chen, Hu Yufeng, et al. Using Non-Meteorological Parameters Tropospheric Delay Estimation Model for Accelerating Convergence of PPP[J]. Geomatics and Information Science of Wuhan University, 2015, 40(2): 188-192)
( 0) |
| [10] |
Shi J, Xu C, Guo J, et al. Local Troposphere Augmentation for Real-Time Precise Point Positioning[J]. Earth, Planets and Space, 2014, 66(1): 1-13 DOI:10.1186/1880-5981-66-1
( 0) |
| [11] |
刘立龙, 姚朝龙, 熊思, 等. 基于插值气压的GPS反演大气可降水量研究[J]. 大地测量与地球动力学, 2013, 33(3): 72-78 (Liu Lilong, Yao Chaolong, Xiong Si, et al. Estimation of Precipitable Water Vapor Based on Interpolated Pressure Data[J]. Journal of Geodesy and Geodynamics, 2013, 33(3): 72-78)
( 0) |
| [12] |
Bekaert D P S, Hooper A, Wright T J. A Spatially-Variable Power-Law Tropospheric Correction Technique for InSAR Data[J]. Journal of Geophysical Research Solid Earth, 2015, 120(2): 1 345-1 356 DOI:10.1002/2014JB011558
( 0) |
| [13] |
Mortier A, Faria J P, Correia C M, et al. BGLS: A Bayesian Formalism for the Generalised Lomb-Scargle Periodogram[J]. Astronomy & Astrophysics, 2014, 573(2): 675-690
( 0) |
| [14] |
Fioranelli F, Ritchie M, Griffiths H. Classification of Unarmed/Armed Personnel Using the Net RAD Multi-Static Radar for Micro-Doppler and Singular Value Decomposition Features[J]. IEEE Geoscience & Remote Sensing Letters, 2015, 12(9): 1 933-1 937
( 0) |
| [15] |
王文胜, 陈伏兵, 杨静宇. 一种基于奇异值分解的特征抽取方法[J]. 电子与信息学报, 2005, 27(2): 294-297 (Wang Wensheng, Chen Fubing, Yang Jingyu. A Method of Feature Extraction Based on SVD[J]. Journal of Electronics & Information Technology, 2005, 27(2): 294-297)
( 0) |