基于GPS-IR的美国中西部地区NDVI时间序列反演

吴继忠，吴 玮



基于GPS-IR的美国中西部地区NDVI时间序列反演

吴继忠，吴 玮

（南京工业大学测绘科学与技术学院，南京211816）

基于AVHRR（advanced very high resolution radiometer）、MODIS（moderate-resolution imaging spectroradiometer）等卫星遥感影像获取的归一化植被指数（normalized difference vegetation index，NDVI）存在大气噪声、土壤背景、饱和度等固有问题。GPS（global positioning system）卫星播发的L波段信号对土壤和植被水分含量变化较为敏感，GPS-IR（GPS-interferometric reflectometry）利用测地型接收机和天线记录GPS反射信号的变化，进而反演测站环境参数。该文研究了利用GPS-IR反演区域NDVI时间序列的方法。采用4个GPS参考站2007－2015年近9 a的连续观测数据，由伪距和相位观测值计算了归一化微波反射指数（normalized microwave reflection index，NMRI），傅立叶变换显示NMRI具有明显的周期特性，其中年周期和半年周期分量普遍较为突出。利用三角多项式拟合剔除NMRI中由积雪和降雨引起的粗差点后，其波动与同时间段内MODIS NDVI的趋势一致。一元线性回归结果显示NMRI与NDVI之间存在显著线性相关，相关系数在0.697～0.818（<0.001），NDVI反演误差的均方根误差在0.059～0.079，表明GPS-IR反演区域NDVI时间序列是可行的，该研究为获取准实时、低成本和高时间分辨率的NDVI提供了新的思路。

模型；植被；遥感；GPS-IR；归一化植被指数；归一化微波反射指数；反演；相关分析

0 引 言

归一化植被指数（normalized difference vegetation index，NDVI）可以反映植被生长情况、覆盖情况、生物量情况和植被种类情况，已广泛用于全球植被状态的定量化研究和应用[1-7]。迄今为止，NDVI都是基于卫星遥感影像来获取，常用的NDVI来源于AVHRR（advanced very high resolution radiometer）、MODIS（moderate-resolutionimaging spectroradiometer）等遥感数据产品。然而，NDVI也存在固有缺陷，特别是：1）受传感器自身因素和大气条件的影响，NDVI数据包含较大的噪声影响[8-12]；2）NDVI数据易受土壤质地、土壤水分含量等土壤背景噪声的影响[13-14]；3）在植被覆盖较高的条件下，NDVI数值易迅速饱和[15-17]。

全球定位系统（global positioning system，GPS）是一种无线电导航定位系统，GPS卫星发射的载波信号位于微波的波段，能穿透植被，且对于植被和土壤含水量的变化较为敏感，经地表植被和土壤反射后的载波信号可被GPS接收机接收并记录。GPS-IR（GPS -interferometric reflectometry）是近些年来发展的一种新方法，其主要思想是分析反射信号物理参数的变化，进而反演出环境变化信息，尤其是与地表土壤和植被相关的信息。Larson等研究利用GPS反射信号用于测量地表积雪厚度，提出了积雪厚度的正演模模型[18-19]。Chew等研究了GPS反射信号与土壤水分含量的回归模型[20-22]。Small等最早研究了GPS反射信号（或多路径信号）与植被状态之间的关系，大量实验结果表明，GPS伪距多路径均方根误差（root mean square error，RMS）与植被高度、植被含水量间都存在良好的负相关性[23]。Wan等研究建立了GPS信噪比数据振幅与草类植物含水量之间的线性回归模型，利用这一模型反演植物含水量误差小于1 kg/m2[24]。Larson等以GPS伪距和相位观测值为基础，提出了归一化微波反射指数（normalized microwave reflection index，NMRI），发现NMRI与植被水分含量之间存在较好的相关性[25]。上述研究已发现GPS反射信号与植被水分含量之间存在相关性，能否利用GPS反射信号来反演NDVI数据则未见相关研究成果。本文研究目的在于分析NMRI的时频分布特性，评价NMRI与NDVI的相关性，探索利用GPS-IR反演NDVI的可行性。论文利用4个不同区域的GPS参考站上长时间观测数据进行了验证和分析，其结论为区域NDVI数据获取提供了新的思路。

1 原理与方法

1.1 反射区与归一化微波反射指数计算

GPS-IR使用的观测设备是常规的测量型接收机及天线，不需要两幅天线，也无需更改天线的朝向。若GPS接收机天线高为，反射信号覆盖区域可用第一菲涅耳区进行描述，其形状是由短半轴和长半轴定义的椭圆

图1 第一菲涅耳区

Fig.1 First Fresnel zones

NMRI是评价反射信号振幅变化的一个综合性指标，其核心是计算1载波上伪距多路径指标MP1的RMS（root mean square）值，MP1的表示为[26]

式中1是1载波上伪距观测值，m；11575.42 MHz，21227.60 MHz；10.19 m，20.24 m；1和2是1、2载波相位观测值，周。与导航定位计算不同，计算MP1的RMS值不需要使用对应整周模糊度的真实值，在无周跳的情况下整周模糊度保持常数不变，其真值的绝对大小不影响RMS的计算，因此周跳的探测处理非常关键。采用经过改进的TurboEdit方法[27]将每颗单日观测数据周跳探测处理，得到重新划分的若干个“干净”弧段，将每个弧段的MP1做去均值化处理后进行弧段合并，计算每颗卫星单日MP1的RMS值并作加权平均，最终得到单日MP1的RMS值。NMRI的计算以MP1的RMS值为基础，其计算方法如下

（3）

式中RMSMP1是单日MP1的RMS值，max(RMSMP1)是RMSMP1序列数值由大到小排列前5%的RMSMP1的平均值，因此NMRI绝大部分值在0～1之间，有少数值为负数。

1.2 归一化植被指数计算

归一化植被指数由可见光波段和近红外波段二者反射率的反差来表征植被生物量的测度，其计算公式为

式中ch1、ch2分别为红光波段和近红外波段经过大气校正的地面反射率。红光波段（波长620～670 nm）处于入射辐射的光谱吸收区，近红外波段（波长841～876 nm）处于入射辐射的光谱反射区。NDVI的数值范围是[-1,1]，数值越大则绿色植被越密集。

2 试 验

2.1 数据来源与计算

GPS数据使用了美国板块边缘观测计划（Plate Boundary Observatory，PBO）中的4个参考站自2007－2015年近9 a的观测数据，数据采样间隔为30 s，卫星截止高度角为5°。由于PBO参考站用途是用于地壳运动监测，其观测环境较好，所用测站周边都是草地覆盖且没有高大障碍物的影响，部分参考站的观测环境视图如图2所示。

4个参考站的概略情况如表1所示。P042和P048在2007－2015年间使用的是Trimble NETRS接收机。P041在2007－2012年使用的是Trimble NETRS接收机，2013年2月更换成Trimble NETR9接收机，由于不同型号接收机对多路径信号的处理性能不同，为保持数据一致性，P041站仅使用2007－2012年的数据，类似的原因，P054使用2012－2015年的数据。将4个参考站在上述时间段内的观测数据下载后，按照式（3）计算单日RMSMP1。

表 1 GPS 参考站概况

目前MODIS NDVI被公认为数据质量较高的植被指数产品之一，本文使用的NDVI数据来自于美国地质调查局陆地过程分布式数据中心的MODIS植被指数产品MOD13Q1，其空间分辨率为250 m，时间分辨率为16 d。利用ENVI软件处理2007－2015年间的影像，各个参考站NDVI值取以测站为中心的3´3像素NDVI的均值。

2.2 数据分析

图3给出了4个参考站的RMSMP1时间序列，可以看出RMSMP1序列总体上具有明显的周期特性，这是植被生长随季节变化的体现；同时还可以看到图上也有明显偏离趋势线的粗差点，这些粗差点主要是受积雪、降雨等天气的影响而产生的[25]。由图3中P042和P054的变化还可以发现，2012年内有部分RMSMP1数值相对于其他年份明显偏大，且持续时间较长，这一现象极有可能是2012年美国中西部地区遭遇极端严重的干旱气候所造成的[28]，图4给出了各个参考站2007－2015年间各年份内累计降雨量，其中2012年的降雨量分别是前5年平均降雨量的68%、46%、103%、72%，可见P041、P042和P054 3个站在2012年的降雨量均低于往年的平均值，相应的RMSMP1数值也明显高于往年，这一特点在图上有清晰地体现，其直接原因是在干旱条件下，地表土壤和植被水分含量很低，伪距多路径效应增大，从而造成RMSMP1值增加[23]。

为揭示RMSMP1序列在频率域的分布规律，利用快速傅立叶变换将离散时间序列转换为离散频谱，图5是快速傅立叶变换后得到的频谱图，由频谱图可以看出，RMSMP1序列中普遍存在明显的周期信号，P041、P048和P054三个站上有年周期信号和半年周期信号，P041和P054两个站的年周期信号比半年周期信号突出，而P048站的半年周期信号比年周期信号突出，振幅高出20%；P042站上仅有明显的年周期信号，其最大振幅对应的信号频率为1.02 a-1，周期约为358 d。

由于RMSMP1容易受到积雪、降雨等影响，在计算NMRI之前需要预处理。最直接的方法是参照气象观测记录，将降雨降雪时间段的数据直接剔除，这种方法只能人工进行，工作量较大。考虑到RMSMP1序列显著的周期性，用三角多项式对RMSMP1序列进行拟合，将拟合残差大于2倍中误差的观测值剔除，再次进行拟合并进行反复迭代，直到拟合残差无超限为止，图6显示了经过误差点检测后的RMSMP1序列，对比图3可看出粗差点有明显的减少。图7显示了P042站上的NMRI，由式（3）可知NMRI与RMSMP1是线性相关，因此NMRI和RMSMP1具有相同的周期特性。

NMRI的时间分辨率为1 d，高于NDVI的时间分辨率。为便于对比，需将二者的采样率进行统一。为此将采样率高的NMRI进行插值处理，用三角多项式插值的方法获得与NDVI采样完全同步的NMRI数据。图8a－8d分别显示了4个参考站上NMRI和NDVI的散点分布图。

从图8可以看出，二者在时域内的总体变化趋势基本一致，NMRI和NDVI的峰值、谷值出现的时间吻合较好，可初步判定两者间存在相关性，这一特征说明由NMRI来反演NDVI是可行的。

2.3 反演模型的建立和验证

为便于模型验证，以时间顺序将试验数据的前60%用于建模，后40%用于模型检验。根据上述思路，分别将4个参考站的NMRI为自变量，NDVI为因变量，进行一元线性回归，表2给出了回归分析的结果。

表2 回归分析结果

由表2可以看出不同测站上的回归系数是不一样的，这与数据量有关，还与参考站硬件处理反射信号的方式有关，硬件的类型乃至固件的版本都会引起RMSMP1数值的尺度变化[29-30]，图3中不同测站上RMSMP1振幅的差异也能印证这一结论，但不同参考站上相关系数没有明显差异，介于0.697～0.818之间，在显著性水平选择为0.001的条件下，检验全部通过，表明NMRI和NDVI之间存在显著相关。考虑到一元线性回归的检验、检验和检验的等价性，检验和检验也必定通过。

将NMRI数据的后40%代入表2中建立的回归模型，得到计算的NDVI，将从遥感影像获取的NDVI值作为真实值，计算每个站上NDVI的反演误差的均方根误差，其大小分别为0.059、0.061、0.079、0.069，其数值与回归分析残差的均方根误差基本接近。为分析误差的区间分布规律，计算NDVI反演相对误差（即反演误差与真实值的比值），以20%的间隔统计每个区间内相对误差所占的比例，其统计结果见图9。

从图9可以看出，NDVI相对误差较小的所占比例较大，以±20%以内的相对误差为例，4个站上所占比例分别为60%、75%、70%、81%，随着相对误差的增大，其比例越来越小，误差分布总体上接近于正态分布，说明反演模型是有效的。

3 结果与讨论

NDVI作为一种重要的遥感参数，迄今为止均来源于AVHRR、MODIS等遥感数据产品，难以克服大气噪声、土壤背景和饱和度等固有缺陷。本文研究了利用GPS-IR（GPS- interferometric reflectometry）反演区域NDVI数据，并通过长时间的实测数据进行了检验和分析，主要结论为：

1）基于GPS-IR生成的NMRI具有明显的周期性，其中年周期和半年周期较为突出；

2）在干旱气候条件下土壤和植被水分含量低，造成反射信号振幅增大，NMRI值相应变小；

3）NMRI和NDVI在时域内波动趋势趋于一致，峰值、谷值出现的时间吻合较好，依此建立了NDVI反演的一元线性回归模型，相关系数在0.728～0.776（<0.05），NDVI拟合残差的均方根误差在0.056～0.091；

由于NMRI的获取更为便捷，利用GPS-IR反演区域NDVI时间序列具备可行性。与遥感卫星大尺度大范围的影像获取不同，GPS-IR记录的是测站周围地表反射信号的特征参数，其作用范围是以测站为中心的圆形区域，有效面积为数千平方米。目前，以GPS为代表的GNSS连续运行参考站已经成为一种空间信息基础设施，参考站的数量和密度在不断增加，GPS-IR可获取的观测范围也在不断增大，全球范围内大量分布的参考站可望成为潜在的NDVI传感器，且具有低成本、准实时、高时间分辨率的优点，不受传感器自身因素和大气条件的影响，但GPS-IR受雨雪、地形条件的影响，GPS-IR与土壤、植被间相互作用的物理机制尚不完全明确，上述问题还需要在后续工作中进一步深入探索。

[1] 边金虎，李爱农，宋孟强，等. MODIS 植被指数时间序列Savitzky-Golay滤波算法重构[J]. 遥感学报，2010，14(4)：725－741.

Bian Jinhu, Li Ainong, Song Mengqiang, et al. Reconstruction of NDVI time-series datasets of MODIS based on Savitzky-Golay filter[J]. Journal of Remote Sensing, 2010, 14(4): 725－741. (in Chinese with English abstract)

[2] Steven M, Malthus T, Baret F, et al. Intercalibration of vegetation indices from different sensor systems[J]. Remote Sensing of Environment, 2003, 88(4): 412－422.

[3] Lunetta R, Knight J, Edirwickrema J, et al. Land-cover change detection using multi-temporal MODIS NDVI data[J]. Remote Sensing of Environment, 2006, 105(2): 142－154.

[4] Montandon L, Small, E. The impact of soil reflectance on the quantification of the green vegetation fraction from NDVI[J]. Remote Sensing of Environment, 2008, 112(4): 1835－1845.

[5] Leon J, van L, Casady G. Using MODIS-NDVI for the modeling of post-wildfire vegetation response as a function of environmental conditions and pre-fire restoration treatments[J]. Remote Sensing, 2012, 4(3): 598－621.

[6] 杜加强，赵晨曦，贾尔恒·阿哈提，等. 近30 a新疆月NDVI动态变化及其驱动因子分析[J]. 农业工程学报，2016，32(5)：172－181.

Du Jiaqiang, Zhao Chenxi, Jiaerheng Ahati, et al. Analysis on spatio-temporal trends and drivers in monthly NDVI during recent decades in Xinjiang, China based two datasets [J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(5): 172－181. (in Chinese with English abstract)

[7] Otto M, Höpfner C, Curio J, et al. Assessing vegetation response to precipitation in northwest Morocco during the last decade: an application of MODIS NDVI and high resolution reanalysis data[J]. Theoretical and Applied Climatology, 2016, 123(1): 23－41.

[8] Roerink G J, Menenti M,Verhoef W. Reconstructing cloud free NDVI composites using Fourier analysis of time series[J]. International Journal of Remote Sensing, 2000, 21(9): 1911－1917.

[9] Atzberger C, Eilers P. Evaluating the effectiveness of smoothing algorithms in the absence of ground reference measurements [J]. International Journal of Remote Sensing, 2011, 32(13): 3689－3709.

[10] Chen J, Jönssonc P, Tamurab M, et al. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golayfilter[J]. Remote Sensing of Environment, 2004, 91: 332－344.

[11] Jennifer N, Gregory J. Noise reduction of NDVI time series: An empirical comparison of selected techniques [J]. Remote Sensing of Environment, 2009, 113(1): 248－258.

[12] 吴文斌，杨鹏，唐华俊，等. 两种NDVI 时间序列数据拟合方法比较[J]. 农业工程学报，2009，25(11)：183－188.

Wu Wenbin, Yang Peng, Tang Huajun, et al.Comparison of two fitting methods of NDVI time series datasets [J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2009, 25(11): 183－188. (in Chinese with English abstract)

[13] Cheng Tao, Rivard B, Sanchez-Azofeifa A. Spectroscopic determination of leaf water content using continuous wavelet analysis [J]. Remote Sensing of Environment, 2011, 115(2): 659－670.

[14] Libonati R, Dacamara C, Pereira C, et al. On a new coordinate system for improved discrimination of vegetation and burned areas using MIR/NIR information [J]. Remote Sensing of Environment, 2011, 115(6): 1464－1477.

[15] 齐述华，王长耀，牛铮. 利用温度植被旱情指数（TVDI）进行全国旱情监测研究[J]. 遥感学报，2003，7(5)：420－427.

Qi Shuhua, Wang Changyao, Niu Zheng. Evaluating soil moisture status in china using the temperature/vegetation dryness index (TVDI)[J]. Journal of Remote Sensing, 2003, 7(5): 420－427. (in Chinese with English abstract)

[16] Sandholt I, Rasmussen K, Andersen J. A simple interpretation of the surface temperaturevegetation index space for assessment of surface moisture status[J]. Remote Sensing of Environment, 2002, 79(2): 213－224.

[17] Liu Fang, Qin Qiming, Zhan Zhiming. A novel dynamic stretching solution to eliminate saturation effect in NDVI and its application in drought monitoring[J]. Chinese Geographical Science, 2012, 22(6): 683－694.

[18] Larson K, Nievinski F. GPS snow sensing: Results from the earthscope plate boundary observatory[J]. GPS Solutions, 2013, 17(3): 41－52.

[19] Nievinski F, Larson K. Forward modeling for near-surface reflectometry and positioning applications[J]. GPS Solution, 2014, 18(6): 309－322.

[20] Chew C, Small E, Larson K, et al. Effects of near-surface soil moisture on GPS snr data: Development of a retrieval algorithm for volumetric soil moisture[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(1): 537－543.

[21] Larson K, Small E, Gutmann E, et al. Using GPS multipath to measure soil moisture fluctuations: initial results[J]. GPS Solutions, 2008, 12(3): 173－177.

[22] 敖敏思，朱建军，胡友健，等. 利用SNR观测值进行GPS 土壤湿度监测[J]. 武汉大学学报：信息科学版，2015，40(1)：117－120.

Ao Minsi, Zhu Jianjun, Hu Youjian, 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 with English abstract)

[23] Small E, Larson K, Braun J. Sensing vegetation growth with reflected GPS signals[J]. Geophysical Research Letters, 2010, 37(12): 1－5.

[24] Wan Wei, Larson K. Using geodetic GPS receivers to measure vegetation water content[J]. GPS Solution, 2015, 19(2): 237－248.

[25] Larson K, Small E. Normalized microwave reflection index, part 1: A vegetation measurement derived from GPS networks[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(5): 1501－1511.

[26] Estey L, Meertens C. TEQC: The multi-purpose toolkit for GPS/GLONASS data[J]. GPS Solutions, 1999, 3(1): 42－49.

[27] 吴继忠，施闯，方荣新. TurboEdit单站GPS数据周跳探测方法的改进[J]. 武汉大学学报：信息科学版，2011，36(1)：29－33.

Wu Jizhong, Shi Chuang, Fang Rongxin. Improving the single station data cycle slip detection approach turboedit[J]. Geomatics and Information Science of Wuhan University, 2011, 36(1): 29－33. (in Chinese with English abstract)

[28] Mallya G, Zhao L, Song X, et al. 2012 Midwest drought in the United States[J]. Journal of Hydrologic Engineering, 2013, 18(7): 737－745.

[29] Bilich A, Larson K. Mapping the GPS multipath environment using the signal-to-noise ratio (SNR)[J]. Radio Science, 2007, 42(6): 1－16.

[30] Vázquez G, Bennett R, Spinler J. Assessment of pseudorange multipath at continuous GPS stations in Mexico[J]. Positioning, 2013, 4(3): 253－265.

Retrieving NDVI in midwestern America using GPS-interferometric reflectometry

Wu Jizhong, Wu Wei

(,,211816)

the NDVI (normalized difference vegetation index) data, routinely derived from the AVHRR (advanced very high resolution radiometer) or MODIS (moderate resolution imaging spectroradiometer) imagery, is a key indicator of vegetation status and a useful parameter in studies of terrestrial vegetation cover, it has been widely used in remote sensing studies to reflect regional and global vegetation dynamics. However, the inherent defects of NDVI, including the atmospheric noise, soil effects and saturation problems are unavoidable, and thus impede further analysis and have a risk to generating erroneous results. Global Positioning System-Interferometric Reflectometry (GPS-IR) is a bistatic radar remote sensing technique that relates temporal changes in reflected GPS signals to changes in environmental parameters surrounding a ground-based GPS site. All GPS satellites transmit signals at L-band, which is similar to those used in active microwave radar applications. L-band signals have a higher correlation with vegetation water content, therefore GPS reflections will be sensitive to water within and on the surface of vegetation, as well as water in soil and snow. The sensing footprint of GPS-IR is on the order of a thousand square meters, which depends on the antenna height and satellite elevation angle. Other than specially-designed antenna or receiver in order to estimate environmental parameters, GPS-IR utilizes geodetic-quality GPS receivers and antennas, which are currently used at many of the already-existing GPS stations. This article presents a new method to retrieve regional NDVI data using NMRI (normalized microwave reflection index), which is an index derived from GPS observations. An experiment was conducted to evaluate the feasibility of the NDVI retrieval using NMRI. In the experiment, continuous GPS observations of four plate boundary observatory GPS reference stations in midwestern America during the interval 2008-2012 and MOD13Q1 product within the same time from MODIS were used. In the first step, the NMRI time series were calculated with the GPS pseudoranges and carrier phase observations preprocessed with an improved Turboedit method, and then NDVI time series were extracted from MOD13Q1 product. In the second step, NMRI and NDVI were compared and analyzed. The temporal fluctuations of NMRI showed a clear periodicity as well as sudden drops, which were not compatible with the gradual process of vegetation change. Fast Fourier transform revealed that the annual and semi-annual periodicities exhibited dominant amplitude. To obtain cleaned NMRI data, trigonometric polynomial fitting method was adopted to remove outliers. A relatively high correlation coefficient between NMRI and NDVI was found, the coefficients of determination varied from 0.697 to 0.818 (with a significance level of<0.001), showing a near linear relationship involving these variables. With regression analysis, a linear retrieve model for NDVI could be established on each reference station, the root mean square of NDVI retrieve errors varied from 0.059 to 0.079. The outcomes of this study suggested that GPS-IR would be almost equally capable of retrieving regional NDVI data, in contrast, GPS-IR had the potential to be in near real time, with low price and high temporal resolution, and what’s more, existing GPS networks around the world had the potential to be the NDVI sensors, which could be regarded as a new opportunity to obtain NDVI data.

models; vegetation; remote sensing; GPS-interferometric reflectometry; normalized microwave reflection index; normalized difference vegetation index; retrieve; correlation analysis

10.11975/j.issn.1002-6819.2016.24.024

P228.4; P237.9

A

1002-6819(2016)-24-0183-06

2016-05-01

2016-11-29

国家自然科学基金资助项目（41504024）；江苏省测绘地理信息科研项目（JSCHKY201413）

吴继忠，男，湖北红安人，博士，副教授，硕士生导师，主要从事卫星导航定位应用研究。南京 南京工业大学测绘科学与技术学院，211816。Email：jzwumail@163.com