CMONOC观测约束下的中国大陆地区MODIS PWV校正

刘 备，王 勇，娄泽生，占 伟

1. 天津城建大学地质与测绘学院，天津 300384； 2. 中国地震局第一监测中心，天津 300180

水汽(precipitable water vapor，PWV)是大气的主要组成部分之一，其对天气和气候的变化有着重要的影响。PWV分布对于研究中小尺度灾害天气，以及短时天气预报中具有重要的意义[1-4]。PWV时空变化导致雷达信号的延迟，降低了合成孔径雷达差分干涉测量(InSAR)的精度[5-7]。常规PWV获取手段具有局限性，如无线电探测法精度高但时间分辨率较低[8-9]。目前，利用GNSS获取PWV技术以及GNSS PWV应用的研究已经趋于成熟[10-14]。对于GNSS反演的PWV精度，国内外多位学者开展了GNSS PWV与无线电探空、水汽辐射计(WVR)PWV比较，获得了两者偏差为1～2 mm的结论[4,15-16]，因此可认为GNSS反演的PWV精度高，可满足气象应用的精度要求。受测站间距数十千米的限制，GNSS PWV存在空间上不连续呈离散分布的特征。利用MODIS影像反演的PWV在空间上是面状连续的，可对大区域PWV进行分析研究，且MODIS PWV的变化趋势与GNSS PWV接近，MODIS PWV可以弥补GNSS PWV在空间上的不连续性。受降水、云层、地表反射光谱不确定等影响，MODIS PWV的相对误差超过10%，如果可以消除地表反射光谱、温度廓线等的不确定性影响，其相对误差可优于5%，因此认为MODIS PWV的精度较低[17-19]。目前对于GNSS PWV与MODIS PWV的相关研究主要集中于小区域以及短时段[20-24]。

融合GNSS PWV与MODIS PWV各自的优势可用于天气预报和InSAR大气校正。为了推广MODIS PWV产品的应用，需要提高MODIS PWV精度。中国大陆构造环境监测网络(CMONOC)积累了全国260个GNSS站点观测资料，本文拟利用CMONOC观测资料开展中国大陆地区MODIS PWV校正研究。首先开展GNSS PWV与MODIS PWV的相关性比较；然后依据不同气候类型的测站分布，构建基于GNSS的不同气候类型MODIS PWV校正模型；最后进行模型的可靠性分析，并进行精度评定。

1 研究数据

本文研究数据包括GNSS PWV和MODIS PWV，研究数据时间为2016年1月—12月。

1.1 GNSS PWV反演

GNSS PWV由CMONOC观测数据经高精度定位定轨软件GAMIT计算获得。CMONOC是一个由260个连续观测和2000个不定期观测站点构成的观测网络。GNSS PWV解算方案如下：解算软件为GAMIT10.6，星历为IGS精密星历，Relax解算模式，卫星高度角10°，测站天顶对流层延迟(ZTD)的解算为每两小时估算一个值，ZTD由静力学延迟(ZHD)与对流层湿延迟(ZWD)组成，ZHD采用Saastamoinen模型结合测站纬度、大地高和气压观测数据精确计算获得，ZWD由ZTD减去ZHD获得[4]。ZWD转化为PWV，可采用公式PWV=∏×ZWD计算获得，转换系数∏可由Bevis提出的公式计算获得[25]。对于由GNSS观测和气象数据反演获得的PWV序列，针对PWV明显异常值，认为是粗差，给予删除预处理。GNSS PWV单位为毫米(mm)，解算精度为1～2 mm，可满足气象应用的要求。

1.2 MODIS PWV反演

MODIS PWV反演方案如下：解算软件为ENVI5，利用MOD05 PWV产品的近红外波段可反演PWV，MODIS PWV单位为mm。MODIS影像可通过NASA网站下载MODIS的MOD05 PWV产品获得论文研究所需MODIS PWV数据(https:∥ladsweb.modaps.eosdis.nasa.gov/search/)。

GNSS PWV的获取时间分辨率为2 h，而MODIS PWV的时间分辨率为每天2个观测值，为了便于两者比较及计算相关性，本文选择将GNSS PWV的时间分辨率抽稀成每天2个观测值，2次观测值时间与MODIS PWV时间最为接近。

2 MODIS PWV校正

利用GNSS PWV进行MODIS PWV，首要任务是确定两者之间的关系，进而选择合适的模型校正方法。本节首先开展GNSS PWV与MODIS PWV的相关性分析，计算两者的偏差；然后选择回归分析方法构建基于GNSS PWV的MODIS PWV校正模型。

2.1 GNSS PWV与MODIS PWV相关性分析

由于CMONOC站点多，若每一个站点都开展GNSS PWV与MODIS PWV的相关性分析和建模工作，模型工作量大且使用也不方便。因此，本文依据中国大陆地区的气候分布类型，将CMONOC站点分别划入到各气候类型(图1)，并依据气候类型分别开展GNSS PWV与MODIS PWV的相关性比较。图2为随机选择两个站点的MODIS PWV与GNSS PWV对比的散点图。

图1 CMONOC站点气候带划分图Fig.1 Map of climate zone stations of CMONOC

图2 MODIS PWV与GNSS PWV对比Fig.2 Comparison between MODIS PWV and GNSS PWV

从图2可以看出MODIS PWV与GNSS PWV之间具有较好的相关性，且两者的演变趋势较为一致，这些特点为利用GNSS PWV对MODIS PWV进行校正提供了基础。

由于我国具有多种气候类型，每种气候类型之间降水差异较大，另外，我国的降水具有明显的季节性，因此，需要对数据进行分类。表1为按照气候类型以及季节划分之后的MODIS PWV与GNSS PWV的相关性统计，表2为两者之间的平均偏差和均方根误差统计。表1中的各气候类型的样本数计算如下：首先统计2016年各站点GNSS PWV与MODIS PWV的对应个数获得站点样本数，然后将各气候类型包含站点的样本数相加获得各气候类型的样本数。对于4个季节的样本数，就是分别求出各季节GNSS PWV与其对应的MODIS PWV序列的对应个数，然后将各气候类型包含站点的季节样本数相加获得各气候类型的季节样本数。

表1 MODIS PWV与GNSS PWV相关性统计

Tab.1 Statistics of correlation between MODIS PWV and GNSS PWV

数据时间高原山地热带季风温带大陆温带季风亚热带季风全年相关性0.654 0.1430.6620.7760.530样本数5045760218341197615355春季相关性0.373 0.3660.4520.6720.402样本数1306176561730343558夏季相关性0.279-0.4180.7470.7990.307样本数1227203555930363808秋季相关性0.569-0.2210.4840.740.381样本数1340193558530754007冬季相关性0.512 0.1680.1760.2400.316样本数1172188507328313682

由表1可以看出MODIS PWV与GNSS PWV具有较好的相关性，除热带季风气候区域外，在其他气候区域两者的全年相关性超过了0.5。高原山地全年的相关性优于4个季节的相关性，其原因在于：MODIS PWV与GNSS PWV分季节的相关性低的原因在于短期小起伏不匹配影响了两者的相关性，而全年的相关性更好的原因是短期小起伏不匹配对于一年的时间尺度来说其影响较小，因此两者的全年相关性优于分季节相关性。高原山地气候类型GNSS PWV与MODIS PWV在夏季的相关性明显低于其他3个季节，热带季风气候类型的GNSS PWV与MODIS PWV在夏、秋两个季节的相关性不同于全年和其他两个季节的相关性，相关性为负，这是因为该区域降水主要集中于夏季和秋季，而降水发生前的云层较厚影响了MODIS PWV精度。

由表2可知，MODIS PWV与GNSS PWV之间存在一些偏差，因高原山地、温带大陆和温带季风气候类型区域的PWV值较小，MODIS PWV与GNSS PWV的偏差与均方根误差相对较小，为10 mm以内。热带季风和亚热带季风区域全年的PWV值较大且波动较大，降水过程多，因此这两个区域MODIS PWV与GNSS PWV的偏差与均方根误差较大，超过10 mm，部分季节超过了20 mm。

通过MODIS PWV与GNSS PWV的误差统计分析，说明MODIS PWV精度不高，有必要进行校正；而MODIS PWV与GNSS PWV的相关性分析，说明可利用GNSS PWV开展MODIS PWV校正研究。

表2 MODIS PWV与GNSS PWV平均偏差和均方根误差统计

2.2 MODIS PWV校正模型

对MODIS PWV校正时，以年、季节两个尺度进行模型构建。其中季节划分是以3—5月为春季；6—8月为夏季；9—11月为秋季；12至次年2月为冬季。模型建立好之后对模型进行精度评定，分别对5个气候类型的MODIS PWV进行两个尺度的校正与建模，模型的一般表达式为

PWVGNSS=A×PWVMODIS+C

(1)

式中，PWVGNSS为GNSS PWV；PWVMODIS为MODIS PWV；A为模型系数；C为常数。为了兼顾模型构建所需样本数较多以及检验样本充分，本文选择80%样本数据用于模型构建，20%样本数据用于模型验证，数据选择为随机选取。表3为MODIS PWV校正模型中各项参数的统计。

MODIS PWV校正模型构建后需要进行模型的可靠性检验，本文采用预留的20%样本进行模型的可靠性验证。表4为对模型验证结果的平均偏差和均方根误差的统计。图3为用直方图绘制的校正前后MODIS PWV与GNSS PWV差值的均方根误差。

表3 MODIS PWV校正模型统计

表4 模型校正后MODIS PWV与GNSS PWV平均偏差均方根误差统计

图3 MODIS PWV(校正前后)与GNSS PWV差值的均方根误差Fig.3 The RMS of the difference between MODIS PWV and GNSS PWV, and the difference between corrected MODIS PWV and GNSS PWV

对比表2、表4和图3，可知：无论哪个气候类型，全年模型的精度都低于春季、秋季，以及冬季这3个季节的模型精度，但是优于夏季模型。此外，热带季风区的夏季和秋季MODIS校正模型的模型系数表现为负数。其原因主要是夏季降水多，MODIS影像受降水以及云层影响较大，而GNSS PWV不受季节、降水等天气因素的影响，所以导致夏季的MODIS PWV精度低于其他季节，并且MODIS PWV与GNSS PWV的相关性受这些因素的影响，表现出异常性，也可以说降水的多少对模型精度有影响，降水越少，模型精度越高。此外，由于GNSS PWV的变化具有年周期、半年周期，以及更小尺度的周期，对全年数据进行建模，由于周期性质复杂，所以导致全年模型的精度低于季节模型的精度。另外，对比分析表2、表4和图3，还可知：热带季风气候类型地区的MODIS PWV校正效果优于其他区域，如热带季风气候区域的全年模型均方根误差校正精度由原来的20.9 mm，达到8.59 mm，校正效果明显。

3 模型推广及验证

对我国进行分区构建MODIS PWV校正模型，这极大地减少了模型数量，但是区域模型是否就能代替站点模型，还需验证。以下对区域模型用于区域内站点MODIS PWV校正的精度进行验证。如果对所有站点开展单站点模型构建并进行检验，由于CMONOC站点达260个，工作量甚大，本文篇幅限制也不容许把所有站点模型及检验结果列出，本文区域模型的检验随机选择两个站点进行。

3.1 模型推广——区域模型用于单站点检验

对以上构建的区域模型进行推广，验证区域模型用于单站点的精度。对各气候类型选择2个站点进行区域模型的检验，检验之前需要构建单站点MODIS PWV校正模型，GNSS站建模的样本数和验证的样本数选取的原则参照与气候类型类似的方法，建模数据约80%，检验数据为20%，模型如表5所示。

表5 单站点MODIS校正模型

选用单站点检验数据用于单站点模型和区域模型的MODIS PWV校正，并与同时段的GNSS PWV和MODIS PWV相比较，图4为部分单站点模型和区域模型的MODIS PWV校正值与GNSS PWV、MODIS PWV的比较。统计区域模型用于单站点校正与单站点模型直接校正的平均偏差和均方根误差(表6)。

图4 MODIS PWV、GNSS PWV、区域模型PWV和站点模型PWV的对比Fig.4 Comparison among MODIS PWV, GNSS PWV, regional model PWV, and single station model PWV

根据表6和图4可以得出结论：利用区域模型对单站点MODIS PWV进行校正，其模型精度与单站点模型的精度相差不大，大部分站点区域模型用于单站点MODIS PWV校正的精度可达到毫米级，某些情况下可以用区域模型代替站点模型对MODIS PWV进行校正。当对精度要求不是特别高的时候，以及站点数据不能满足要求时，可以利用区域模型对站点进行MODIS PWV的校正工作。另外，区域模型对区域内MODIS PWV校正，不需要对每个站点都单独进行校正工作，这大大减少了人力与物力。

3.2 MODIS PWV校正模型精度评定

MODIS PWV是面状连续的，上文对MODIS PWV的校正是基于点状数据进行的，无法将MODIS PWV面状连续的优点利用起来，而ArcGIS软件能够通过一定的插值方法将点状数据以面状的形式显示出来。

表6 区域模型用于单站点校正与单站点模型直接校正的误差统计

插值方法有多种，中国地形起伏差异较大，直接采用克里金或反距离加权方法插值效果不佳。Delaunay三角网的优点是结构良好，数据结构简单，与不规则的地面特征和谐一致，可表示线性特征和叠加任意形状的区域边界。本文选择区域气候模型在GNSS站点校正，然后通过采用构建Delaunay三角形，再在三角形内进行水汽内插的方式获得中国大陆地区的MODIS PWV插值结果。利用ArcGIS中的空间分析功能对插值后的图像进行作差，这样就可以得到面状的差值图，弥补前面仅对点状数据进行分析的不足，充分利用MODIS PWV的面状优势。图5为2016年11月24日UTC 03:40时刻的GNSS PWV、MODIS PWV，以及校正后MODIS PWV之间的差值图。

图5 GNSS PWV、MODIS PWV，以及校正后MODIS PWV之间的差值Fig.5 Difference diagram among GNSS PWV, MODIS PWV and the calibration-PWV of MODIS

由图5可知：亚热带季风气候、温带季风气候的MODIS PWV校正后效果优于校正前精度，温带大陆性气候、高原山地气候和热带季风气候校正前后的差异不显著。为了精确获得各类型气候的校正效果，对各气候区所有GNSS站点校正前后的MODIS PWV与GNSS PWV之差，分别统计校正前后的平均偏差和均方根误差(表7)。

表7 模型校正前与校正后MODIS PWV与GNSS PWV平均偏差和均方根误差统计

Tab.7 Mean deviation and root mean square error statistics between, corrected and uncorrected MODIS PWV and GNSS PWVmm

高原山地热带季风温带大陆温带季风亚热带季风均方根误差校正前0.13-5.49-0.71-2.53-6.63校正后3.558.511.09-0.252.81平均偏差校正前3.5411.882.113.2711.46校正后1.348.511.631.114.79

由表7可知：从平均偏差来看，各气候类型校正后的MODIS PWV优于校正前MODIS PWV；就均方根误差而言，对于亚热带季风和温带季风气候类型来说，校正后效果优于校正前，温带大陆气候校正前后相当，高原山地和热带季风气候校正后效果低于校正前。原因分析如下：因选择的数据时间为2016年11月24日，处于秋末冬初季节，高原山地和温带大陆气候区域的水汽值较低，水汽值低于10 mm，水汽变化小，校正效果较差。而温带季风气候和亚热带季风气候水汽值相对较高，因此校正效果明显。热带季风气候水汽值虽高，由于该区域站点较少，改善效果欠佳。

4 结 论

本文通过对CMONOC的GNSS资料处理、分析，然后利用GNSS PWV对MODIS PWV进行校正，并结合气象学等多学科，开展了MODIS PWV校正研究，得出了以下结论：利用GNSS PWV对MODIS PWV进行校正，区域模型和站点模型的校正后精度均能够达到毫米级，达到了气象业务应用的要求；利用区域模型对单站点MODIS PWV进行校正，其模型精度与单站点模型的精度相差不大，区域模型用于单站点PWV校正时大部分站点的精度可达到毫米级，某些情况下可以用区域模型代替单站点模型对MODIS PWV进行校正；利用Delaunay三角网插值分析，可以将校正后的点状数据扩展为面状数据，插值结果可以对气象预报起到一定的参考作用。

鉴于热带季风气候地区MODIS PWV与GNSS PWV的相关性低，为提高模型精度，可选用其他更合适的函数关系来构建基于GNSS 的MODIS PWV校正模型。利用模型校正的MODIS PWV用于获取全国MODIS PWV变化图时，为保证精度，插值需要考虑温度、气压、地形等因素对该地区PWV变化的影响。

致谢：感谢中国大陆构造环境监测网络中心为本文研究提供GMONOC站点GNSS观测数据和气象数据！