留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

利用MODIS温度产品进行秩修正滤波FRF时空插值

段悦 舒红 胡泓达

段悦, 舒红, 胡泓达. 利用MODIS温度产品进行秩修正滤波FRF时空插值[J]. 武汉大学学报 ● 信息科学版, 2016, 41(8): 1027-1033. doi: 10.13203/j.whugis20140495
引用本文: 段悦, 舒红, 胡泓达. 利用MODIS温度产品进行秩修正滤波FRF时空插值[J]. 武汉大学学报 ● 信息科学版, 2016, 41(8): 1027-1033. doi: 10.13203/j.whugis20140495
DUAN Yue, SHU Hong, HU Hongda. Using Fixed Rank Filtering to Make Spatio-Temporal Interpolation of MODIS Temperature[J]. Geomatics and Information Science of Wuhan University, 2016, 41(8): 1027-1033. doi: 10.13203/j.whugis20140495
Citation: DUAN Yue, SHU Hong, HU Hongda. Using Fixed Rank Filtering to Make Spatio-Temporal Interpolation of MODIS Temperature[J]. Geomatics and Information Science of Wuhan University, 2016, 41(8): 1027-1033. doi: 10.13203/j.whugis20140495

利用MODIS温度产品进行秩修正滤波FRF时空插值

doi: 10.13203/j.whugis20140495
基金项目: 

国家自然科学基金 41171313

湖北省自然科学基金 2014CFB725

苏州市科技计划应用基础研究 SYG201319

详细信息

Using Fixed Rank Filtering to Make Spatio-Temporal Interpolation of MODIS Temperature

Funds: 

The National Natural Science Foundation of China 41171313

Hubei Provincial Natural Science Foundation of China 2014CFB725

Suzhou Science and Technology Program of Applied Basic Research SYG201319

More Information
    Author Bio:

    DUAN Yue, master, specializes in spatio-temporal statistical analysis. E-mail:yue_duan@126.com

    Corresponding author: SHU Hong, PhD, professor. E-mail:shu_hong@whu.edu.cn
图(5) / 表(1)
计量
  • 文章访问数:  1023
  • HTML全文浏览量:  44
  • PDF下载量:  309
  • 被引次数: 0
出版历程
  • 收稿日期:  2015-04-13
  • 刊出日期:  2016-08-05

利用MODIS温度产品进行秩修正滤波FRF时空插值

doi: 10.13203/j.whugis20140495
    基金项目:

    国家自然科学基金 41171313

    湖北省自然科学基金 2014CFB725

    苏州市科技计划应用基础研究 SYG201319

    作者简介:

    段悦, 硕士, 研究方向为时空统计分析。yue_duan@126.com

    通讯作者: 舒红, 博士, 教授。shu_hong@whu.edu.cn
  • 中图分类号: P237

摘要: 普通克里金方法具有空间结构探索及插值分析(给出插值结果及其精度)的功能。但是,绝大部分的克里金方法主要用于空间插值,在时空插值方面的研究还较少。综合考虑具体的实验区域和观测数据构建基函数,使用秩修正滤波(fixed rank filtering,FRF)方法对MODIS平流层温度数据进行时空插值预测并将其结果与秩修正克里金(Fixed rank Kriging,FRK)方法的插值结果进行对比分析。实验结果表明,在空间数据(空间点)整体分布均匀且有已知点的情况下,FRK方法预测的数据精度更高,略优于FRF;而对于较大空间范围内缺失数据的情况,考虑温度在时间维上具有一定的相关性,FRF方法在缺失空间信息时能够引入更多时空信息从而获得较其他方法更高质量的插值结果。

English Abstract

段悦, 舒红, 胡泓达. 利用MODIS温度产品进行秩修正滤波FRF时空插值[J]. 武汉大学学报 ● 信息科学版, 2016, 41(8): 1027-1033. doi: 10.13203/j.whugis20140495
引用本文: 段悦, 舒红, 胡泓达. 利用MODIS温度产品进行秩修正滤波FRF时空插值[J]. 武汉大学学报 ● 信息科学版, 2016, 41(8): 1027-1033. doi: 10.13203/j.whugis20140495
DUAN Yue, SHU Hong, HU Hongda. Using Fixed Rank Filtering to Make Spatio-Temporal Interpolation of MODIS Temperature[J]. Geomatics and Information Science of Wuhan University, 2016, 41(8): 1027-1033. doi: 10.13203/j.whugis20140495
Citation: DUAN Yue, SHU Hong, HU Hongda. Using Fixed Rank Filtering to Make Spatio-Temporal Interpolation of MODIS Temperature[J]. Geomatics and Information Science of Wuhan University, 2016, 41(8): 1027-1033. doi: 10.13203/j.whugis20140495
  • 运用大量不完整且掺杂有噪声的数据进行插值分析是一项很具挑战性的工作。当面对大量的空间数据集时,由于传统的克里金方法需要对n×n阶协方差矩阵进行求逆运算,而当数据量n增长到一定程度时,该运算在一般计算机环境下难以实现,因此直接应用传统的方法进行计算有一定困难[1]。为此,国内外学者进行了大量研究,并提出了若干解决方案,其中包括协方差矩阵逼近[2]、克里金方程逼近[3]、协方差矩阵裁剪[4]和近似迭代[5]等方法。特别地,Cressie和Johannesson[6, 7]提出了一种较有影响的方法,即通过建立一个空间随机效应(spatial random effects, SRE)模型来构建具有高度灵活性的空间协方差矩阵,并对协方差矩阵精确求逆。进一步,Cressie等[8]将SRE模型推广到时空模型中并发展成为时空随机效应(spatio-temporal random effects, STRE)模型。该模型通过降低其中需进行复杂运算矩阵的秩来达到高效计算的目的,并同时使用卡尔曼滤波实现时间维上的快速递归更新。Cressie等[8]将这种通过构建STRE模型进行插值计算并获取其对应的估计误差的方法称为秩修正滤波(FRF)方法。

    秩修正克里金(fixed rank Kriging, FRK)方法是基于构建SRE模型进行插值的方法,它与FRF方法均能够达到简化计算并提高插值精度的目的。二者间区别在于,FRK仅使用某t时刻的空间数据对这个时刻的待测点进行估计,而FRF集合当前t时刻和t时刻之前的数据信息来对这个时刻各待测点值进行估计,即综合利用了时空维信息。先前实验已表明,FRF[9]及FRK[7]方法均能在较短的时间内处理大量数据,其算法的高效性已有考证,而带有时空信息的FRF方法在插值精度上能否较FRK有所提高?文献[8]通过构造一维模拟数据对比二者的插值精度,并且仅使用两个时刻的气溶胶光学厚度数据进行FRF时空插值实验;而文献[9]首先同样使用构造的一维模拟数据对比两种方法的插值精度,随后选用了多个时刻的气溶胶光学厚度数据经由FRF方法进行时空插值实验,获得多个时刻的插值结果。文献[8, 9]仅通过构造一维模拟数据的方式对两种方法进行对比实验,没有使用实际的观测数据。作为全球气候系统变化的一个重要组成部分,本文基于平流层温度的MODIS时空数据,分别使用FRF和FRK方法进行插值计算得到多个时刻的温度插值结果,并对这两种方法的插值精度进行对比分析。

    • FRF和FRK两种插值方法能够做到在保证精度的同时提高计算效率的原因主要在于,选取合理的空间基函数向量组St(以下简称基函数),保证了精度;以及协方差矩阵Σt求逆运算时,公式变换使得计算复杂度降低,效率提高。经一系列推导过程[8],可得:

      (1)

      式中,sDRdt∈1, 2, …,D为数据所在的空间范围;Z(s; t)为带有测量误差的观测值;Y(s; t)是去除测量误差后的真值;ε是高斯白噪声且ε(s; t)~N(0, σε2vε(s; t)),其中Y(s; t)又可分为两部分,即Y(s; t)=μt(s)+v(s; t),其中μt(·)是(时空)确定性趋势项;v(s; t)是空间变异随机项。作为随机效应向量,ηt是STRE模型中空间与时间信息的关键:

      (2)

      式中,St(·)≡(S1, t(·), …, Sr, t(·))′代表一组空间基函数向量(r值已知);ηt(·)是均值为0的高斯随机向量,其r×r阶协方差矩阵用Kt表示;ξ(s; t)代表由降秩引入的误差,ξ(s; t)~N(0, σξ2vξ(s; t)),本文中所有ηt(·)、ξ(·; ·)和ε(·; ·)均相互独立;第二式是一个一阶向量自回归过程,Ht+1是一阶自回归(或传播)矩阵,ζt+1独立于ηt,均值为0且方差vart+1)≡Ut+1

      给定数据Z(·; 1), …, Z(·; t),对Y(s0; t)进行预测的过程被称作滤波过程。对Z(t)求方差,则有:

      (3)

      nt表示第t个时刻Z(t)的数据量,var(Z(t))是一个nt×nt的正定矩阵;IntVε均为nt×nt阶对角阵,其中Int对角线上系数由指示函数确定,Vε对角线上nt个系数值vt(si, t)由测量误差方差构成。

      FRK和FRF方法均根据Sherman-Morrison-Woodbury[10]公式对Σt=StKtSt+Dt求逆, 这也是二者实现高效计算的关键。求逆过程如式(4)所示:

      (4)

      式中,协方差矩阵Σt具有高度灵活性,并且Σt的求逆运算只需对固定阶数(r×r)的正定矩阵Ktnt×nt阶对角阵Dt分别求逆获得,如此则在很大程度上简化了计算,将计算的时间复杂度从O(nt3)降为O(r3)。

    • FRK方法仅使用第t时刻获取的数据Z(t)对Y(s0; t)作预测,由Cressie和Johannesson在文献[6, 7]中给出:

      (5)

      式中,kt(s0)≡StKtSt(s0)+ct(s0), ct(s0)=σξ2(I(s0=s1, t), …, I(s0=snt, t))′。同时给出与待插值点对应的FRK预测误差方差:

      (6)
    • FRF方法与FRK方法最大的不同在于前者通过卡尔曼滤波将时间与空间信息结合进行预测。在FRF方法中,通过卡尔曼滤波器不断进行自回归运算更新ηt|tPt|t(这里的ηt|tPt|t与FRK中的ηtKt略有区别,其中t|tE(ηt|Z(1), …, Z(t)),均方预测误差矩阵Pt|tE[(ηt)(ηt)′])取值,可获得前t个时刻的插值结果:

      (7)

      (8)

      式中, (s0)=ct(s0)′(StPt|t-1St′+Dt)-1(Z(t)-μ(t)-Stt|t-1)。利用卡尔曼滤波器算法不断更新ηt|tPt|t的过程可用式(9)完成:

      (9)

      式中,前两个运算是卡尔曼滤波器预测公式,t|t-1指的是利用上一状态预测的结果,t-1|t-1是上一状态最优结果; 而Pt|t-1Pt-1|t-1分别是这两种状态结果对应的均方预测误差矩阵; Ut是系统过程的协方差矩阵。利用式(9)中前两个运算完成系统预测之后还需结合当前状态的测量值得到当前状态的最优估值t|t,这一过程由式(9)中第三个运算完成,Gt为卡尔曼增益。至此,已得到t状态下最优的估值t|t,但为了令卡尔曼滤波不断运行直到系统过程结束,还需更新t状态下的误差矩阵Pt|t,由式(9)中最后一式完成。如此过程,卡尔曼滤波算法便可由式(9)中的5个步骤自回归的运算得到所需的ηt|tPt|tGt值。

    • 平流层温度的变化趋势已被认为是研究地球大气温度垂直廓线变化及研究臭氧和水汽变化的关键。平流层温度的变化会引起其下层的对流层温度、气压场和风场的变化,进而引发对流层的换流和气候的变化,因而对平流层温度的变化及趋势估计引起了广泛的重视[11]。本实验选取平流层中下层250 hPa处的MODIS遥感反演温度数据作为实验数据,8 d为一周期(即每8 d内获取的数据取平均),起止时间为2013年1月1日~2013年3月5日,共8组。该数据的空间分辨率为1°×1°(部分点位缺失数据)。本文仅选用经度60°N~30°E,纬度40°S~50°N范围内的数据进行实验(共8 100个规则格网点数据,记为原始观测数据Z(t)all),将该范围记为D,经度20°W~20°E、纬度10°S~20°N的范围记D0DD0范围记为DV。前7 h分别在D内随机选取4 000个已知点数据作为使用的观测数据Z(t)obs,而第8个时刻在DV范围内选取4 000个已知点数据。图 1为8个时刻分别选取的观测点数据。

      图  1  8个时刻的观测点数据

      Figure 1.  Observation Data at Eight Times

    • 本文选用MODIS反演的平流层温度产品Z1~Z8作为观测数据,分别使用FRF和FRK两种方法对与Z1~Z8对应的T=1, …, 8时刻进行时空和空间插值计算,得到T=2, …, 8时刻的插值结果。

      在此过程中,由§1可知, 合理的基函数选取是保证插值精度的关键,也是计算空间协方差的关键步骤,而协方差又是获得最终结果的核心部分。为了能够最大程度表达温度的空间变异性,前人建议在不同空间尺度上选取基函数。其选择类型有很多种[12-17],其中Shi和Cressie在文献[14]中选用小波函数构造基函数St,文献[8, 9]中也同样采用此种方法构造St。而文献[7, 15]在确定中心点后使用双平方函数构造基函数St,但由于其中部分不同尺度上所选取的中心点有重合,可能会使预测误差图中呈现微小的振荡[16]。因此,曾有学者针对此问题进行改进使得无重合中心点出现[16, 17]。本文采用此种改进后的方法获取中心点。式(10)为通过双平方函数获取中心点后构造基函数的过程:

      (10)

      式中, vj(l)代表三个不同分辨率下基函数的中心点,j=1, …, nl=1, 2, 3,u代表对应的已知点。令该分辨率上任意两个基函数中心之间的最短距离为d,则rl=1.5d

      考虑实验区域和观测数据量的大小,在实验区域内分别从三个不同尺度上均匀选取中心点来构造基函数,以保证该基函数中包含有不同尺度上的空间变异信息,并且控制中心点数量(基函数数量)与数据量大小相适应。因为在保证n > r的情况下,若数据量n小但基函数数量r过多,会导致通过降秩来达到降低时间复杂度的效果减弱,所以当数据量不足够大时,选择的中心点数量也不宜过多。选择结果如图 2所示,分别用三种不同的符号表示,共9+36+144=189个中心点。获得中心点后,通过式(10)的计算即可得到基函数。

      图  2  不同尺度上中心点分布情况

      Figure 2.  Distribution of Center Points in Three Resolutions

      基函数St已知,FRK方法采用期望最大化(expectation-maximization, EM)算法迭代计算出K和σξ2的值即可代入计算其协方差矩阵求逆式(4)[15],完成插值计算。FRF方法则是通过给定协方差矩阵初始值,并计算得到P0|0HU。令HU不随时间变化,给定初值η0|0,将η0|0P0|0输入卡尔曼滤波器,不断更新ηt|tPt|t直到系统过程结束,根据式(7)、式(8)计算得到最终结果。

    • T=1,…,8时刻,D范围内,使用FRF方法对温度插值得到T=2,…,8时刻的结果,图 3所示为T=2, 4, 6时刻FRF的温度插值结果及标准差,其标准差最大达到0.3℃左右。在T=2~8时刻分别利用FRK方法进行空间插值计算。这里仅选取T=8时刻FRK的空间插值结果与通过FRF方法获得的T=8时刻时空插值结果进行对比,并绘制预测值、预测标准差及绝对误差(原始观测数据Z(8)all与预测结果Ŷ(8)之差的绝对值)图(图 4)。

      图  3  T=2, 4, 6时刻FRF预测值及标准差

      Figure 3.  FRF Predictions and Standard Errors for T=2, 4, 6

      图  4  三种插值方法对比

      Figure 4.  Comparison of Three Interpolation Methods

      图 4显示,整体上通过FRF方法获得时空插值结果的预测标准差与FRK的预测标准差大致相同,前者略优于后者。再来观察这两种方法分别在有已知数据DV和无已知数据D0范围内的插值情况。DV范围内FRF和FRK插值结果相近,但后者相较前者更为平滑;D0范围内二者的预测结果有所不同,FRF方法的插值结果呈现出与DV范围内气温分布较为一致的情况,变化趋于平稳,而FRK方法D0范围中间偏右部分一块类似“心型”不规则区域的温度明显高出周边。对比绝对误差图可知,整体上FRF方法的预测效果最优;在DV处其绝对误差略小于FRK,而在D0范围内,FRK的绝对误差明显大得多,尤其是类似“心型”区域,绝对误差由“心型”中心向外扩张,逐渐变小。

      因此,在观测数据点间距过大时(如D0范围),空间相关性非常小,以致预测的准确性很低。如图 4中FRK在T=8时刻的预测值,距其小矩形边框越远,预测越不准确,在中心处达到最大。实验随机选取4 000个观测点来插值8 100个点,4 000个观测点整体上较均匀地分布在该实验范围内,而在通过式(10)构造基函数时rl的最小值约为2 047.6 km,最大值约为6 041.1 km,绝大部分的观测点与周围最近的点间距都要小于rl的最小值,因此在本实验中,随机选取观测点对插值结果的影响几乎可以忽略不计。

      利用均方根误差(root-mean-square error, RMSE)指标对比各方法的预测结果误差,计算结果如表 1。分别对D0DV空间范围内计算均方根误差,观察到在无已知数据D0的范围内,FRF预测结果的RMSE值相对较小,而FRK方法的误差是它的4倍之多。可见,对于较大空间范围内缺失数据的情况,使用FRF方法更为妥当,即将前面几个时刻该区域的空间信息添加进来辅助预测;在随机分布有已知点处的DV范围内,FRK预测结果的RMSE值仅为0.25,是FRF的47.5%,可见对于整体范围数据分布较为均匀的情况,使用FRK方法能获得更高的预测精度;而在整个D范围内,FRF预测结果的RMSE值为0.646 2,FRK为它的2.6倍,可见从整体效果来看FRF方法更好。

      表 1  均方根误差

      Table 1.  Root-Mean-Square Error

      RMSE/℃ D0 DV D
      FRF 1.111 8 0.524 6 0.646 2
      FRK 4.402 4 0.249 3 1.710 0

      图 5T=8时刻的原始观测数据Z(8)all,将其与以上两种方法的插值结果对比易看出在D0范围内FRF的插值结果与实际更为接近,DV范围内FRK插值效果更好。

      图  5  T=8时刻原始观测数据

      Figure 5.  Original Observation Data at T=8

      为评价插值精度,结合图 1前8个时刻的数据,可明显观察到FRF方法的插值结果融合了前几个时刻的空间信息,因此更加确定FRF方法在D0范围内的预测结果为三者中最优。

      本实验使用的工作站处理器为Intel Core i3 540 3.07 GHz,内存为8 GB。在T=8时刻使用4 000个观测点经由FRK方法对8 100个点进行预测,其中模型参数估计的计算耗时53.0 s,获得预测结果及其误差标准差耗时0.9 s。而利用FRF方法从T=1时刻开始将每个时刻的4 000个观测点逐步加入卡尔曼滤波器计算得到每一时刻8 100个点的插值结果(T=1时刻除外),其中模型参数估计的计算耗时103.2 s(σξ2HUP),获得预测结果及误差标准差耗时34.1 s。

    • 本文结合具体的实验区域和实验数据在三个不同空间尺度上均匀选取中心点并控制合理的中心点数量,通过双平方函数构造空间基函数。分别利用FRF和FRK方法对T=1,…,8时刻250 hPa处的MODIS温度数据进行插值预测的对比实验。在空间随机分布有已知点的情况下,FRK方法预测的数据精度更高,略优于FRF;而对于较大空间范围内缺失数据的情况,由于温度具有一定的时间维相关性,因此若在某时刻出现此情形,可将前几个时刻该区域的空间信息添加进来以便增加有用信息辅助预测,这种情况下,使用FRF方法可获得更高质量的结果。综合以上两种情况,实验表明FRF方法具有更高的实用性,整体上能够达到较好的效果。同时,文献[8, 9]的实验结果也都证明了FRF较FRK的插值精度有所提高,在区域数据缺失处尤为显著。

      对比FRF和FRK这两种方法,其相同点首先在于提高计算效率的原理相同,均使用Sherman-Morrison-Woodbury公式,FRF方法是作为FRK在时间维上的扩展;而它们的不同点在于FRF的扩展加入了卡尔曼滤波器不断自回归更新参数,将t时刻之前多个时刻的信息代入参与共同预测该t时刻信息,这也使得FRF较FRK能够获得更高的插值精度。在计算机负载能够承受的情况下,这两种方法在处理大数据量时均能明显提高计算效率,而当获取的数据量或者数据量与面积之比较小时,计算效率反而会降低甚至增加计算量,同时插值精度也会受到一定影响。因此,当需要进行插值计算时还需要根据数据量等因素合理选择模型。

      实验过程中,二者均使用同一基函数,所以基函数的选择对于该对比实验基本无影响。但作为FRK和FRF方法的核心,不仅基函数的选择对插值精度有影响,其数量也直接影响计算耗时,因此后续研究主要围绕基函数展开。

参考文献 (17)

目录

    /

    返回文章
    返回