基于最优插值的土壤含水量遥感反演缺失数据插补

李 军， 董 恒， 王 祥， 游 林

(1.中国矿业大学(北京)地球科学与测绘工程学院，北京 100083；2.武汉理工大学资源与环境工程学院，武汉 430070)

0 引言

近几十年来，在全球气候变化的背景下干旱频繁发生，其影响范围之广，持续时间之长，受灾程度之重都十分罕见。据统计，每年因旱灾造成的全球经济损失高达60～80亿美元，远远超过了其他灾害[1]。我国是一个农业旱灾频发的国家，1950—2008年间，平均每年受旱面积为2 157万hm2，成灾面积为956万hm2，因旱灾损失粮食为158万t。而且干旱灾害发生次数逐年增加，特别是近几年接连不断发生在西南地区、山东省和长江中下游地区的特大干旱，给农业生产和经济发展带来巨大损失[2-3]。目前，利用遥感技术加强旱情实时监测已成为一个迫切的需求[4]。学者们针对各类卫星遥感数据，通过获取各类地表理化性质[5-9]监测旱情状况及变化，其中土壤含水量是反映干旱的一个最重要也是最直接的指标。通过遥感反演土壤含水量是干旱遥感监测业务中必不可少的环节。

遥感传感器在成像过程中，不可避免地会受到云、雪、气溶胶和传感器自身性能等因素的影响，使得遥感数据在空间上呈现不连续，即存在缺失数据的情况[10]。在反演土壤含水量时，若不能修补缺失数据将严重影响反演结果的实用价值与成果图件的美观效果，限制了产品的深入应用[11-12]。现有数据插补方法主要分为2大类。第一类方法是采用滤波法根据时间序列数据估算缺失数据。简易可行的一种滤波方法是Holben提出的最大值合成法[13]，即按照事先规定的时间间隔，将间隔内序列中各点按数值大小排序，选取各点序列中最大值作为新图像值。另一种效果较好的滤波方法是最佳指数斜率提取法[14-15]，采用滑动时间窗口识别时间序列值中的突变点，并替代噪声值。此后又相继出现了经验正交函数分解法[16]、时间窗口线性内插法[17]、基于非对称的高斯函数拟合法[18]、Savitzky-Golay滤波法[19]和局部最大值拟合法[20]等滤波方法。此类方法虽然具有较高的稳定性和准确性，但依赖于时间序列遥感数据，增加了产品生产的经济成本。并且很多遥感产品往往不具备完整的历史遥感资料，使得该方法的应用受到了限制。另一类方法是采用空间插值法根据邻近数据估算缺失数据。冯益明等和俞晓群等在解译遥感影像和分析海表叶绿素时利用Kriging插值法恢复了缺失数据[21-22]；杨金红等研究了利用线性插值法去除MODIS遥感影像中的条带噪声[23]。此外，还提出了邻行插值法[24]、GIS辅助数据下的影像缺失信息恢复方法[25]、基于纹理合成技术的数据修补方法[26]和矩匹配法[27]等典型的遥感缺失数据恢复方法。但这些方法只适用于特定数据源或图像条件，难以应用于业务部门的推广应用中。

针对上述各类方法存在的局限性，本文旨在围绕土壤含水量遥感反演中的数据缺失问题，提出利用最优插值法，结合气象站点的实测数据，对缺失像元进行插值与填充，为干旱监测业务流程的完整性及准确性提供技术支持。

1 最优插值法

最优插值法是Gandin在1963年提出的一种客观分析的方法[28]。最优插值理论被广泛应用于气象要素场的客观分析及数值天气预报和气象站网的设计中[28-29]。最优插值法里各已知点的内插权重不是预先确定的，而是根据它们对插值点所做贡献的大小，以一定数学方法求取的，用此方法计算的权重进行内插所产生的标准误差比任意选择的权重所造成的误差都小。从统计意义上讲是均方插值误差最小的线性插值法。一些学者已将最优插值法用于遥感作物指数计算[30]和海温预报资料同化[31-32]中，取得了较好效果。

本文将最优插值应用于土壤含水量Z的插补中，用下标a表示分析值，g表示初始值，ob表示气象站点处土壤含水量的观测值，下标i表示目标格点的序号，k表示气象台站的序号，则将离散的气象站的观测数据转化成均匀分布的网格分析值时的线性内插公式为

(1)

(2)

(3)

(4)

(5)

最优插值的本质就是要寻找出使分析误差方差达到最小时的权重因子Pk，根据∂E′/∂Pk=0，(k=1,2,…，n)，可导出

(6)

从而可得

(7)

于是Pk可以由式(6)求得，其中μik根据需求选定某种函数形式。

假设μkl随k和l这2点间的距离标量r呈指数衰减，相关系数可表示为

μkl=exp(-rkl/a) 。

(8)

式中a为常数系数，其取值越大，相关系数趋于0的距离(最大影响距离)也越大，此时对于密集的近处测站，它们的相关系数值相差很小，就不能很好反映出相关系数随距离变化的特性，最优权重系数也近似相等，将失去求相关的意义，故a的取值要根据各个测站的实际分布情况而定。本文研究实例中a取1 500 km。

2 研究方法

2.1 缺失像元识别

在土壤含水量的遥感反演过程中，遥感图像上云、雪覆盖等因素会导致反演结果出现缺失像素值。对缺失像元进行插补前需要识别缺失像元的空间位置。反演土壤含水量的干旱指数DI可表示为各类参数变量的函数,即

DI=f(X1,X2,…，Xn)，

(9)

式中Xi为计算干旱指数DI所使用的遥感或物理参数变量，如不同波段的地表反射率、植被指数和地表温度等。具体方法如下：

1)对于每个变量Xi，生成一幅标记有效值区域的掩模影像MaskXi。对于第k个像素，若像素值Xi,k∈[LXi,UXi]，则MaskXi,k=1，否则MaskXi,k=0，式中LXi和UXi分别为变量Xi的最小可能取值与最大可能取值。

2)合成最终掩模影像，对于任一像素，掩模影像合成公式为

Maskk=MaskX1,kMaskX2,k…MaskXn,k。

(10)

3)标记干旱指数DI中的缺失像元，令

DIfinalk=DIkMaskk+bgValue(1-Maskk) ，

(11)

式中：DIk为第k个像素处直接计算而得的干旱指数值；DIfinalk为最终的干旱指数值；bgValue为干旱指数的缺失像元填充值(如-999)。此时得到的干旱指数结果中，数值为bgValue的像元即为缺失像元。

2.2 缺失数据插补

最优插值法的核心是利用观测站点的观测值与初估值的偏差修订目标格点的值。针对干旱监测的特点，本文以气象站点为观测站点，以气象站的土壤体积含水量实测值为观测值。基于最优插值的缺失数据插补方法包括2个步骤：

1)初始背景场的构建。利用各气象站点Stationk(k=1,2,…,n)历年实测数据，计算气象站处土壤体积含水量的平均值ck，通过二阶反距离权重法插值构建土壤体积含水量的初始背景场，对于任一像元i，土壤体积含水量的初估值为

(12)

(13)

式中：Pk,i为二阶反距离权重；Dk,i为第k个站点离像元i的距离。

确定站点观测值选择策略后，便可对每个缺失像元，根据式(5)逐一计算最优插值的权重P，最后根据式(1)便可得到缺失像元的插值结果。

3 研究实例

3.1 研究区及数据源概况

本文以中国西北的宁夏回族自治区为实验区(如图1所示)，使用该区域内的MODIS地表反射率产品(MOD09A1)和反照率产品(MCD43B3)，由美国NASA的EOS数据中心(https: //wist.echo.nasa.gov/wist-bin/api/ims.cgi)下载，数据产品的时相为2001年3—9月间，是宁夏地区主要农作物的生长季节。使用前将原始数据转换为UTM投影，由红光波段与近红外波段计算归一化植被指数(normalized difference vegetation index,NDVI)，构建基于NDVI-Albedo特征空间的植被条件反射率干旱指数(vegetation condition albedo drought index,VCADI)[33]作为干旱遥感反演指数。此外，搜集实验区内16个国家级气象台站(图1)同时期的土壤水分自动监测数据。通过选取同时具有地面实测土壤体积含水量与遥感VCADI指数的样本，采用回归分析构建基于VCADI的土壤体积含水量反演模型，进而反演各时相土壤含水量。

图1 研究区及气象观测站点

3.2 缺失数据插补效果分析

本文选取了具有不同程度缺失数据的8个时相土壤体积含水量结果，即2001年的第97，105，113，153，161，185，249和257天的数据。应用本文方法，以16个气象台站的监测数据和遥感反演数据为基础，识别缺失像元并插补缺失数据。图2展示了应用插补方法恢复缺失像元前(左)后(右)的对比效果。

(a) 年积日97 (b) 年积日105

(c) 年积日113 (d) 年积日153

(e) 年积日161 (f) 年积日185

(g) 年积日249 (h) 年积日257

图22001年宁夏土壤体积含水量插补前(左)与插补后(右)对比

Fig.2SoilmoisturecontentofNingxiain2001before(left)andafter(right)interpolationprocess

从图2中可以看出：

1)对于像元缺失相对较少的时相，如97,113和153，插补整体效果较好，插补值与周围的旱情具有很好的空间连续性； 在时相113和153产品上，缺失数据既覆盖了高值区域，又覆盖了低值区域，插补后的数据仍然较好地反映了缺失区域的土壤体积含水量差异。

2)对于存在大面积缺失像元的时相，如105，161，185和249，本文方法仍能根据气象站点实测值估算缺失数据，有效地保证了产品的完整性，也在一定程度上体现了旱情的空间分布特征。

3)对于遥感数据几乎完全缺失的时相(如257)，采用了初始背景场对缺失值进行插补，其结果虽然没有分析价值，但仍然满足了业务化监测流程的数据完整性要求。

3.3 精度分析

为了定量分析最优插值方法用于土壤含水量缺失数据插补的精度及与其他方法的对比效果，本文选择了多期具有完整土壤含水量的数据作为参考(图3(a))，人工模拟矩形成块缺失数据(图3(a)中矩形框)，然后分别使用反向距离加权插值(图3(b))、Kriging空间插值(图3(c))和最优插值法(图3(d))插补缺失数据，并与原始土壤含水量进行对比分析。

(a) 土壤含水量原始数据 (b) 反向距离加权插值结果 (c)Kriging空间插值结果 (d)最优插值结果

图3最优插值与其他插值方法对比

Fig.3Comparisonbetweenresultsofoptimuminterpolationandotherinterpolationmethods

以土壤含水量原始数据为参照，计算3种插值方法用于各期数据的均方根误差的平均值。对于反向距离加权插值法，东、南、西、北矩形框缺失数据插补结果的均方根误差分别为0.036 1，0.039 2，0.088 1和0.079 8； 对于Kriging插值法，4块区域的均方根误差分别为0.036 9，0.034 8，0.074 4和0.108 8； 对于最优插值法，4块区域的均方根误差分别为0.031 4，0.032 7，0.054 9和0.057 8。可以看出，Kriging插值法在空间异质性较高的数据区域出现缺失时插补效果不如反向距离加权插值法，而在异质性低的区域效果较好，二者效果总体相当。与这2种方法相比，最优插值法得到的插补结果精度更高，虽然插补结果与实际数值仍然有一定偏差，但在成片数据缺失且不依赖长时序序列数据的情况下，能基本反映出数据缺失区域的整体状况。

此外，最优插值法还可以通过人工去掉数据点，模拟不同数据缺失率(10%，20%，30%和40%)的含水量(图4)。

(a) 缺失率为10%(b) 缺失率为20%(c) 缺失率为30%(d) 缺失率为40%

图4不同缺失率数据模拟

Fig.4Simulationofmissingdataatdifferentlevels

与原始数据进行对比，以均方根误差评价插补精度。当缺失率为10%时，插补误差为0.014 7； 当缺失率为20%时，插补误差为0.031 5； 当缺失率为30%时，插补误差为0.044 1； 当缺失率为40%时，插补误差为0.073 2。可以看出，在缺失比例很高时，插补结果仍然能达到较高精度,主要原因是模拟缺失数据是随机和均匀分布的，更利于依据周围数据插补真实结果。

4 结论

1)本文针对土壤含水量遥感反演过程中存在因云雪等导致的数据缺失问题，提出应用最优插值法，综合利用气象站点的实测和遥感反演数据为参考观测值，对缺失像元进行插值与填充。该方法不依赖于长时间序列的历史遥感数据，并且实测数据独立于遥感数据，降低了在成块缺失数据区域应用邻近数据估算所带来的误差。

2)选取宁夏回族自治区为研究区，以MODIS数据为数据源，根据VCADI指数与地面实测数据建立的土壤含水量反演模型，得到该区域具有不同程度数据缺失的多个时期土壤含水量。结合16个国家级气象台站观测数据，利用最优插值法对土壤体积含水量中的缺失数据进行插补。结果表明，该方法无论对少量数据缺失或大面积缺失都具有较好的效果。此外，本文分别利用反向距离加权插值法、Kriging空间插值法和最优插值法，对多期完整土壤含水量结果图进行矩形成块缺失数据含水量模拟，结果表明最优插值法具有更高插值精度。还通过不同数据缺失率模拟，测试了在高缺失率情况时最优插值法的插补效果。

3)本文方法适合于在分析区域内有足够数量实测站点且其空间分布相对均匀的情况，此外，面积较大的成块缺失数据会导致作为插补模型参考的实测数据与遥感反演数据均出现缺失，虽然本文方法给出了利用初始背景场进行填充的解决方案，但仅是为了确保产品的完整性，结果的使用需要慎重。本研究区主要包含的土壤类型为砂壤土、青壤土、黄壤土、粘土和黑垆土[34]，在进行数据插补时未区分土壤类型，虽然结果满足要求，但插补方法对不同土壤类型的敏感性是今后需要研究的重要问题。

参考文献(References)：

[1] Wilhite D A,Buchanan-Smith M.Drought as Hazard:Understanding the Natural and Social Context[M]//Wilhite D A.Drought and Water Crises Science,Technology,and Management Issues.Boca Raton:Taylor and Francis Group,2005.

[2] 翁白莎,严登华.变化环境下中国干旱综合应对措施探讨[J].资源科学,2010,32(2):309-316.

Weng B S,Yan D H.Integrated strategies for dealing with droughts in changing environment in China[J].Resources Science,2010,32(2):309-316.

[3] 周 磊,武建军,张 洁.以遥感为基础的干旱监测方法研究进展[J].地理科学,2015,35(5):630-636.

Zhou L,Wu J J,Zhang J.Remote sensing-based drought monitoring approach and research progress[J].Scientia Geographica Sinica,2015,35(5):630-636.

[4] 杨绍锷,闫娜娜,吴炳方.农业干旱遥感监测研究进展[J].遥感信息,2010(1):103-109.

Yang S E,Yan N N,Wu B F.Advances in agricultural drought monitoring by remote sensing[J].Remote Sensing Information,2010(1):103-109.

[5] 吴代晖,范闻捷,崔要奎,等.高光谱遥感监测土壤含水量研究进展[J].光谱学与光谱分析,2010,30(11):3067-3071.

Wu D H,Fan W J,Cui Y K,et al.Review of monitoring soil water content using hyperspectral remote sensing[J].Spectroscopy and Spectral Analysis,2010,30(11):3067-3071.

[6] 姚云军,秦其明,赵少华,等.基于ΔTs-Albedo光谱信息的土壤水分监测新指数研究[J].光谱学与光谱分析,2011,31(6):1557-1561.

Yao Y J,Qin Q M,Zhao S H,et al.New index for soil moisture monitoring based on ΔTs-Albedo spectral information[J].Spectroscopy and Spectral Analysis,2011,31(6):1557-1561.

[7] 孙 灏,陈云浩,孙洪泉.典型农业干旱遥感监测指数的比较及分类体系[J].农业工程学报,2012,28(14):147-154.

Sun H,Chen Y H,Sun H Q.Comparisons and classification system of typical remote sensing indexes for agricultural drought[J].Transactions of the Chinese Society of Agricultural Engineering,2012,28(14):147-154.

[8] 赵 娟,张耀鸿,黄文江,等.基于热点效应的不同株型小麦LAI反演[J].光谱学与光谱分析,2014,34(1):207-211.

Zhao J,Zhang Y H,Huang W J,et al.Inversion of LAI by considering the hotspot effect for different geometrical wheat[J].Spectroscopy and Spectral Analysis,2014,34(1):207-211.

[9] 王鹏新,吴高峰,白雪娇,等.基于Landsat数据的条件植被温度指数升尺度转换方法[J].农业机械学报,2015,46(7):264-271.

Wang P X,Wu G F,Bai X J,et al.Up-scaling transformation methods for vegetation temperature condition index retrieved from Landsat data[J].Transactions of the Chinese Society for Agricultural Machinery,2015,46(7):264-271.

[10] Bradley B A,Jacob R W,Hermance J F,et al.A curve fitting procedure to derive inter-annual phenologies from time series of noisy satellite NDVI data[J].Remote Sensing of Environment,2007,106(2):137-145.

[11] Cihlar J,Ly H,Li Z Q,et al.Multitemporal,multichannel AVHRR data sets for land biosphere studies-Artifacts and corrections[J].Remote Sensing of Environment,1997,60(1):35-57.

[12] 李 儒,张 霞,刘 波,等.遥感时间序列数据滤波重建算法发展综述[J].遥感学报,2009,13(2):335-341.

Li R,Zhang X,Liu B,et al.Review on methods of remote sensing time-series data reconstruction[J].Journal of Remote Sensing,2009,13(2):335-341.

[13] Holben B N.Characteristics of maximum-value composite images from temporal AVHRR data[J].International Journal of Remote Sensing,1986,7(11):1417-1434.

[14] Viovy N,Arino O,Belward A S.The best index slope extraction (BISE):A method for reducing noise in NDVI time-series[J].International Journal of Remote Sensing,1992,13(8):1585-1590.

[15] Lovell J L,Graetz R D.Filtering pathfinder AVHRR land NDVI data for Australia[J].International Journal of Remote Sensing,2001,22(13):2649-2654.

[16] Ganzedo U,Alvera-Azcárate A,Esnaola G,et al.Reconstruction of sea surface temperature by means of DINEOF:A case study during the fishing season in the Bay of Biscay[J].International Journal of Remote Sensing,2011,32(4):933-950.

[17] Park J,Tateishi R.Correction of time series NDVI by the method of temporal window operation[C]//Proceedings of 1998 Asian Conference on Remote Sensing,1998.

[18] Jonsson P,Eklundh L.Seasonality extraction by function fitting to time-series of satellite sensor data[J].IEEE Transactions on Geoscience and Remote Sensing, 2002,40(8):1824-1832.

[19] Savitzky A,Golay M J E.Smoothing and differentiation of data by simplified least squares procedures[J].Analytical Chemistry,1964,36(8):1627-1639.

[20] Akhter S,Sarkar I,Rabbany K G,et al.Adapting the LMF temporal splining procedure from serial to MPI/Linux clusters[J].Journal of Computer Science,2007,3(3):130-133.

[21] 冯益明,雷相东,陆元昌.应用空间统计学理论解译遥感影像信息“缺失”区[J].遥感学报,2004,8(4):317-322.

Feng Y M,Lei X D,Lu Y C.Interpretation of pixel-missing patch of remote sensing image with Kriging interpolation of spatial statistics[J].Journal of Remote Sensing,2004,8(4):317-322.

[22] 俞晓群,马翱慧.基于Kriging空间插补海表叶绿素遥感缺失数据的研究[J].测绘通报,2013(12):47-50.

Yu X Q,Ma A H.The spatial interpolation of missing remote sensing data in sea surface chlorophyll-a using Kriging[J].Bulletin of Surveying and Mapping,2013(12):47-50.

[23] 杨金红,顾松山,程明虎.插值法在去除MODIS遥感影像条带噪声中的应用[J].气象科学,2007,27(6):604-609.

Yang J H,Gu S S,Cheng M H.Application of interpolation method in destriping MODIS images[J].Scientia Meteorologica Sinica,2007,27(6):604-609.

[24] 熊贤成,杨春平,敖明武,等.MODIS影像条带噪声行的判断及去除研究[J].遥感技术与应用,2015,30(3):540-546.

Xiong X C,Yang C P,Ao M W,et al.A research on judging and removing stripe noises of MODIS image[J].Remote Sensing Technology and Application,2015,30(3):540-546.

[25] 陈仁喜,李鑫慧.GIS辅助数据下的影像缺失信息恢复[J].武汉大学学报(信息科学版),2008,33(5):461-464.

Chen R X,Li X H.Restoring lost information on remote sensing images based on accessorial GIS data[J].Geomatics and Information Science of Wuhan University,2008,33(5):461-464.

[26] 陈仁喜,李鑫慧,李盛阳.纹理合成技术在遥感影像缺失信息恢复中的应用[J].遥感信息,2009(5):15-18,86.

Chen R X,Li X H,Li S Y.Texture synthesis and it’s application in restoring missing information on remote sensing images[J].Remote Sensing Information,2009(5):15-18,86.

[27] 朱小祥,范天锡,黄 签.《神舟三号》成像光谱仪图像条带消除的一种方法[J].红外与毫米波学报,2004,23(6):451-454.

Zhu X X,Fan T X,Huang Q.Method to destripe imaging spectroradiometer data of SZ-3[J].Journal of Infrared and Millimeter Waves,2004,23(6):451-454.

[28] Gandin L S.Objective Analysis of Meteorological Fields[M].Gidromet:Almaty,Kazakhstan,1963.

[29] 李建通,张培昌.最优插值法用于天气雷达测定区域降水量[J].台湾海峡,1996,15(3):255-259.

Li J T,Zhang P C.Optimum interpolation method used for measuring regional precipition with weather Radar[J].Journal of Oceanography in Taiwan Strait,1996,15(3):255-259.

[30] 申广荣,田国良.基于GIS的黄淮海平原旱灾遥感监测研究——作物缺水指数模型的实现[J].生态学报,2000,20(2):224-228.

Shen G R,Tian G L.Remote sensing monitoring of drought in Huanghe,Huaihe and Haihe Plain based on GIS-the calculation of crop water stress index model[J].Acta Ecologica Sinica,2000,20(2):224-228.

[31] 朱 江,徐启春,王赐震,等.海温数值预报资料同化试验I.客观分析的最优插值法试验[J].海洋学报,1995,17(6):9-20.

Zhu J,Xu Q C,Wang C Z,et al. Assimilation experiment of prediction data of sea surface temperature I:Objective analysis of optimum interpolation[J].Acta Oceanologica Sinica,1995,17(6):9-20.

[32] 马寨璞,井爱芹.动态最优插值方法及其同化应用研究[J].河北大学学报(自然科学版),2004,24(6):574-580.

Ma Z P,Jing A Q.Dynamic interpolation and its application in data assimilation[J].Journal of Hebei University(Natural Science Edition),2004,24(6):574-580.

[33] Ghulam A,Li Z L,Qin Q M,et al.Exploration of the spectral space based on vegetation index and albedo for surface drought estimation[J].Journal of Applied Remote Sensing,2007,1(1):013529.

[34] 张秀珍,刘秉儒,詹硕仁.宁夏境内12种主要土壤类型分布区域与剖面特征[J].宁夏农林科技,2011,52(9):48-50,63.

Zhang X Z,Liu B R,Zhan S R.Distribution area and profile features of twelve soil types in Ningxia[J].Ningxia Journal of Agriculture and Forestry Science and Technology,2011,52(9):48-50,63.