顾及映射函数误差的对流层延迟两步估计法

范昊鹏，孙中苗，张丽萍，刘晓刚

1. 信息工程大学地理空间信息学院，河南 郑州 450052； 2. 西安测绘研究所，陕西 西安 710054； 3. 地理信息工程国家重点实验室，陕西 西安 710054

对于甚长基线干涉测量(very long baseline interferometry，VLBI)、全球导航卫星系统(global navigation satellite system，GNSS)等空间大地测量技术，对流层延迟是其主要误差源之一。在利用这些测量技术求解大地测量参数(包括测站坐标、速度、地球定向参数等)的过程中，一般做法是：将对流层延迟中的主项，即天顶干延迟部分作为延迟的先验值；因干延迟计算精度较高，可达毫米级，远小于天顶湿延迟值，所以可将天顶湿延迟作为干延迟的改正值估计得出；最终，干、湿天顶延迟分别与对应的映射函数相乘并求和得到斜路径延迟[1-2]。以下如不作特殊说明，统一将天顶干、湿延迟简称为干、湿延迟。

这种估计策略并未考虑映射函数的误差，如果选用GMF(global mapping function)、NMF(Niell mapping function)等经验映射函数模型，则干、湿分量有可能存在10-2、10-1量级的误差[3-4]，此时与干、湿延迟相乘会产生厘米级甚至分米级的斜路径延迟误差，进而对测站位置等大地测量参数的解算造成不可忽略的影响[5]。虽然VMF(Vienna mapping function)系列模型[6-7]相比GMF和NMF具有更高的精度和时间分辨率，但VMF模型需要利用气象数据事后反算得到，因而在实际观测中无法实时获取。这对未来大地测量产品的快速、高精度发布较为不利[8]。

本文以GMF经验模型为例，分析了该映射函数的误差特点，利用统计方法构建了该误差随高度角变化的函数模型，进而提出了顾及映射函数误差的对流层延迟两步估计法，并用试验证明了该方法的有效性。

1 传统估计方法

传统估计方法中，对流层延迟的数学模型如下

L=ZHD·mfh(el)+ZWD·mfw(el)+

maz(el)·[Gnscos(α)+Gewsin(α)]

(1)

式中，L表示斜路径对流层延迟；ZHD、ZWD分别表示天顶对流层干、湿延迟；mfh(el)、mfw(el)表示高度角为el时干、湿分量的映射函数；maz(el)表示高度角为el时对流层水平梯度映射函数，具体形式为[9]

(2)

式中，α为方位角；Gns、Gew分别被定义为水平梯度的南北、东西分量。

此处需说明的是：①一般情况下，干延迟可按照Saastamonien模型计算[10]，当已知地表气压时，认为干延迟误差可忽略；其中，地表气压可通过气象仪器观测获取，或利用数值气象模型等方法得到。本文采用的是后者，具体是欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts，ECMWF)提供的分层格网气象数据。除此之外，后文用于评定精度的映射函数、湿延迟等参数的参考值，也根据该数据计算得到。ECMWF数据根据多源气象信息综合得到，在国际气象学研究中具有十分权威的地位，常被学者用作评判对流层延迟精度的参考数据。②映射函数通常是将待求地区的位置和年积日代入GMF等全球经验模型计算得到(具体公式可参照文献[11—12])，并忽略其误差。此时观测方程为

Li=[mfw(eli)maz(eli)cos(αi)maz(eli)sin(αi)]

(3)

式中，下标i表示第i次观测的量。

但事实上，映射函数并非准确。以GMF映射函数模型为例，绘出西安地区(大地坐标为34°10′41.74″N，108°59′8.16″E，504.427 m)2016年5°高度角时，GMF与映射函数参考值如图1(此时暂不考虑方位角对映射函数的影响)。其中映射函数参考值采用ECMWF分层格网气象数据并根据射线追踪方法反算得到[13-14](以下试验的参考值计算方法与之一致)。

图1 西安地区2016年5°高度角时干、湿分量映射函数Fig.1 Hydrostatic and wet mapping function of Xi’an area during 2016 at elevation of 5°

统计映射函数各分量误差见表1。

表1西安地区5°高度角的映射函数各分量误差

Tab.1ErrorsofeachmappingfunctioncomponentsofXi’anareaatelevationof5°

映射函数分量最小值最大值平均值标准差干分量-0.00730.03560.01980.0061湿分量-0.72211.87320.2800.3541

统计一年当中误差分布频数直方图见图2。

图2 映射函数误差频数直方图Fig.2 Histogram of mapping function errors

从图1、表1可知，高度角为5°时，西安地区干、湿分量映射函数均存在明显的系统偏差，湿分量的RMS超过0.3，且最大误差超过1.8；另外由图2可知，当高度角一定时，一年当中的干、湿分量映射函数误差大体呈正态分布。

2 顾及映射函数误差的估计方法

顾及映射函数误差的对流层延迟数学模型如下

maz(el)·[Gnscos(α)+Gewsin(α)]

(4)

Li=[maz(eli)cosαimaz(eli)sinαiZHD

(5)

此处需说明两点：①为避免观测方程中出现两种不同待估参数相乘的现象，式(4)中引入湿延迟先验值，并在式(5)中忽略湿延迟改正与映射函数改正相乘的项；②观测方程待估参数中的映射函数改正项，即δmfh(eli)、δmfw(eli)包含了随高度角而变的量，需将高度角相关项剥离，方可按照最小二乘原理求解。故下文首先对δmfh(eli)和δmfw(eli)(也可将其视为映射函数误差)的变化特性进行分析。

2.1 映射函数误差特性分析

在我国境内选取西安、武夷、乌鲁木齐、日喀则、海口5个地区，分别代表温带季风、亚热带季风、温带大陆、高原和热带季风5种主要的气候类型区域进行研究。各地具体位置见表2。

表2 各地区大地坐标

统计2016年各地区映射函数误差的年均值和标准差随高度角变化的情况见图3。

从图3可知，该误差的平均值和标准差随高度角增大而迅速减小，当高度角大于40°时，映射函数误差对斜路径延迟的影响将降至亚毫米量级，基本可忽略不计；同时，该误差呈现出较强的地域性，即随地理位置不同而变，尤其是湿分量误差平均值的变化最明显。以下以西安地区为例进一步研究。

统计西安2013—2016年各年映射函数误差的平均值和标准差随高度角的变化情况见图4。

图3 各地区干、湿分量映射函数误差的平均值和标准差随高度角的变化Fig.3 Variation of average and standard deviation of hydrostatic and wet mapping function errors with elevation in different areas

从图4可知，相比地域性特点，映射函数误差的时域性表现出较强的一致性，即该误差的平均值和标准差各年纵向相比变化趋势十分接近；同时观察各曲线的形状特征，可考虑利用指数函数对其拟合，具体形式如下，其中a1、a2、b1、b2为待拟合系数

y=a1exp(a2x)+b1exp(b2x)

(6)

统计各年份、各分量的平均值、标准差曲线的拟合系数见表3。可看出，各年拟合系数十分接近，与上文结论吻合；故考虑利用前3年拟合系数的平均值代表后续年份的情况，即对2013—2015年的a1、a2、b1、b2分别取平均，以其平均值代表2016年的系数，并统计这种方法的拟合精度，结果见表4。

表3 映射函数误差拟合系数及残差标准差统计

续表3

图4 2013—2016年西安地区干、湿分量映射函数误差的平均值和标准差随高度角的变化Fig.4 Variation of average and standard deviation of hydrostatic and wet mapping function errors with elevation in Xi’an during 2013—2016

表4利用2013—2015年拟合系数平均值代替2016年系数后的拟合精度

Tab.4Statisticsoffittingaccuracywhenthecoefficientsof2016arechangedwith3-yearaverageonesfrom2013to2015

映射函数分量精度最小值最大值平均值标准差干分量平均值-7.1397e-48.6899e-5-9.8463e-51.9445e-4标准差-1.1015e-43.5432e-44.0411e-51.3045e-4湿分量平均值-0.0044 0.0017 -5.2254e-40.0014 标准差-0.00929.0738e-4-0.00240.0025

(7)

(8)

(9)

(10)

(11)

表5 统计参数各项系数

2.2 观测方程改化

根据上文分析，可将观测方程式(5)写作如下形式

Li=[maz(eli)cos(αi)maz(eli)sin(αi)

(12)

此处需说明的是，在利用式(12)求解之前，应首先根据式(3)，即传统估计方法解算得到湿延迟先验值，也即ZWD0，同时得ZWD0的权逆阵；进而采用部分参数加权平差[15-16]进行解算。最终参数平差值的计算公式为

(13)

(14)

式中，V为残差向量；n为观测向量维数；t1为非随机参数维数，此处为4，即除去具有先验信息的湿延迟后的待估参数个数。

此即为顾及映射函数误差的两步估计法，即第一步利用传统估计方法(式(3))对湿延迟进行估计；随后根据部分参数加权平差，利用式(12)得到最终结果。为叙述方便，后文统一将其简称为两步法。

3 在新一代VLBI网观测中的应用效果仿真

新一代VLBI全球观测系统(VLBI global observing system，VGOS)建成后，未来大地测量参数的解算精度有望达到1 mm量级[8,17]，将成为毫米级地球科学领域的重要技术基础，然而大气延迟误差依旧会成为制约VGOS测量精度的重要因素[18-19]。故下文通过仿真，研究传统方法与两步法的估计效果。

根据未来VGOS网的指标[20-21]，每台站单天的观测将达到上千次。故本文首先模拟生成1440个观测方向(平均每分钟1次)，即1440组方位角(0°～360°)、高度角(5°～90°)序列，各观测方向具体在测站上空的分布及其在不同高度角、方位角区间的占比可参见图5、图6；随后以西安为例，根据该地区2016年年积日185 d的分层格网气象数据(数据来自ECMWF，原始分辨率为6 h，此处通过线性插值生成时间分辨率为1 min的数据)，利用射线追踪方法计算各方向斜路径延迟值，并将其作为真实值；接着在该值中加入标准差为4 ps的白噪声和方差为0.3 ps2/sec的随机漫步噪声(分别代表时延测量误差和时钟误差[22]，该误差量级与VGOS设备精度一致)，以此作为虚拟观测量；最后分别采用传统方法和两步法每60历元估计一次结果，统计两种方法得到的斜路径延迟残差(见图7，精度统计见表6)和湿延迟估计值(见图8，精度统计见表7)。

表6 西安地区斜路径延迟残差

表7 西安地区湿延迟估计误差

图5 模拟观测方向在测站上空的分布外圈为方位角轴，向径为高度角轴Fig.5 Distribution of simulated observing orientations over the station outer circle represents the azimuthal axis, while inside diameter represents the elevation angle axis

图6 模拟观测方向在各高度角、方位角区间的占比Fig.6 Percentages of simulated observing orientations over various elevation angle and azimuth ranges

从图7、表6可知，相比传统法，两步法的斜路径延迟残差被大大削减，RMS精度提升至6 mm以内，改善幅度达46.22%，同时最大误差降低54.25%；从图8、表7知，湿延迟估计精度有小幅提升，但提升效果十分有限，改善程度在亚毫米级，这主要与第二次平差中湿延迟初值较多依赖第一次平差结果有关。

图7 两种方法计算的斜路径延迟误差Fig.7 Slant path delay errors of two different methods

图8 湿延迟参考值与两种方法的估计值Fig.8 Reference wet delay values and estimated ones by two methods

同理，将上文方法用于其他地区，得到的结果见表8、表9。可见该方法在不同地区均有相近的改善效果。

4 结 语

本文首先对对流层映射函数的误差进行研究，利用以往若干年的分层气象格网数据构建了映射函数误差随高度角变化的函数模型，并提出了顾及映射函数误差的对流层延迟两步估计法。试验表明：①本文构建的映射函数误差模型，能精确逼近研究地区映射函数误差在各高度角的情况，可为分析各地区映射函数误差提供依据；②两步法可有效削减斜路径延迟残差，改善幅度达50%左右，这意味着采用该方法可将对流层延迟对VLBI观测结果的影响削弱50%，同时对对流层湿延迟的估计结果也有一定程度的改善。另外，虽然目前VMF系列映射函数精度较高，但模型系数的发布往往滞后期较长，而本文提出的方法可将映射函数改正值和湿延迟一并估计，有利于近实时性地获取映射函数结果。然而本文方法尚存在不足之处，如构建两步法的观测方程时，认定湿延迟改正的量级足够小，故忽略了该改正与映射函数改正相乘的项；但如果先验湿延迟精度不够高(如大于1 cm)，则可能会引入毫米级甚至1 cm以上的误差。因此如何选用更合理的湿延迟先验值是进一步提高估计精度的关键。

表8 各地区斜路径延迟残差

表9 各地区湿延迟估计误差