基于卫星遥感与回归克里格的地面日降雨量预测方法

摘要

本发明公开了一种基于卫星遥感与回归克里格的地面日降雨量预测方法。首先利用卫星遥感快速获取数据，根据降雨量分级，建立各级降雨量的地面观测值、TRMM、DEM及地理位置间的回归关系，得到各级回归估计值和回归残差。其次分析各级回归残差的空间集聚度，对回归残差进行趋势去除，并进行回归残差的克里格插值，得到1km的各级回归残差空间分布特征。然后分别将各级回归估计值与回归残差相加，得到1km的各级降雨的地面预测值。最终将各级降雨的地面预测值合并，得到1km的日降雨量预测值。本发明能准确地预测地面降雨量的时空分布特征，提高了地面日降雨的预测精度，又提高了预测的空间分辨率，解决水利部门进行地面降雨预报的关键问题。

说明书

技术领域

本发明涉及一种地面降雨预测方法，尤其是涉及一种基于卫星遥感与回归克里格的地面日降雨量预测方法。 

技术背景

降水时空变化较为剧烈，准确预测降水对水文、气象及灾害预报等具有重要意义。目前常用的降水数据主要通过地面观测站和卫星测雨等方法获取。地面观测站作为一种常规的测量降水手段，具有应用广泛、精度较高且技术成熟的特点，但是地面观测站布网密度满足相关研究工作的质量需要，且地面观测站观测的降水量仅代表地面观测站及周边一定范围的降水状况，在表述大面积降水分布特征时可能引起偏差。卫星测雨遥感技术能够提供较高时空分辨率的降水数据，覆盖空间范围更广，克服了地面雨量站和测雨雷达的局限，为监测全球降水变化提供有力的数据支持。 

近年来，随着气象卫星技术的发展，全球尺度高时空分辨率的测雨卫星应运而生，其中利用热带降雨测量卫星（Tropical Rainfall Measuring Mission,TRMM）上搭载的微波成像仪（TRMM Microwave Imager,TMI）及测雨雷达（Precipitation Radar,PR）验证降雨已经取得了长足的进步。TRMM降雨卫星是第一颗专门用于定量观测热带、亚热带降雨的气象卫星，能够以较高的时空分辨率，提供覆盖全球50°S～50°N以内的区域的降雨数据。 

国内外学者在利用TRMM卫星数据进行降水分析和降水资料验证等方面的工作上做出了大量贡献，Islam等（Islam M N，Uyeda H.Use of TRMM in determining the climatic characteristics of rainfall over Bangladesh.Remote Sensing of Environment，2007，108：264—276.）利用TRMM卫星数据在孟加拉国的对比结果精度达到96％，表明使用TRMM卫星反演降水能够很好的反映区域的降水状况。Huffman等（Huffman G J,Robert F A,David T B.The TRMM Multisatellite Precipitation Analysis(TMPA):Quasi-global,multiyear，combined-sensor precipitation estimates at fine scales.Journal of Hydrometeorology,2006,8(1):38-55.）在澳大利亚进行了TRMM卫星 数据与气象观测数据、雷达观测数据的一致性检测，结果表明TRMM卫星数据与观测数据在日、月尺度上均具有较高的线性相关特征。杨传国等（杨传国,余钟波,林朝晖.基于TRMM卫星雷达降雨的流域陆面水文过程.水科学进展，2009,20(4)：461-466.）利用TRMM卫星数据研究流域尺度陆面水文过程，结果表明利用TRMM卫星数据模拟的结果与利用地面观测数据模拟的结果相似。白爱娟等（白爱娟,方建刚,张科翔.TRMM卫星资料对陕西及周边地区夏季降水的探测.灾害学,2008,23(2):41-45.）采用TRMM卫星数据在陕西及其周边地区与地面气象站观测资料进行了对比，发现TRMM卫星数据对区域夏季降水过程的具有良好的探测能力。 

以上研究表明TRMM卫星数据在预测降雨方面具有很强的适应性和可行性。但是，TRMM卫星的原始分辨率较低（空间分辨率为0.25°，约30km），在预测区域尺度降雨方面具有局限性和误差，因此需要针对TRMM数据进行空间尺度的转换，以得到较高空间分辨率的降雨量预测值。Jia等（Jia,S.,Zhu,W.,L,A.,&Yan,T.(2011).Astatistical spatial downscaling algorithm of TRMM precipitation based on NDVI and DEM in the Qaidam Basin of China.Remote Sensing of Environment,115,3069-3079.）利用NDVI和DEM数据，对柴达木盆地的TRMM降雨数据进行降尺度研究，发现降尺度后的r²从0.72上升到0.96，取得较好成果。马金辉等（马金辉，屈创，张海莜，夏燕秋.2001-2010年石羊河流域上游TRMM降水资料的降尺度研究.地理科学进展,2013,32(9):1423-1432.）利用1km DEM数据，采用回归方程结合残差的Thiessen插值方法，将TRMM年降水量的空间分辨率提高到1km。经检验发现，降尺度后的数据精度有一定程度的提高。 

作为一种间接估测降雨的工具，TRMM卫星不仅存在空间尺度上的局限性，而且在应用到不同的地形及降雨条件时，同样具有很大程度上的不确定性及误差。Tian等（Tian,Y.D.,and Peters-Lidard,C.D.,2007,Systematic anomalies over inland water bodies in satellite-based precipitation estimates.Geophysical Research Letters,34,L14403,doi:10.1029/2007GL030787.）发现TRMM卫星数据应用到海洋及南美洲等地时，不确定性为40–60%，而应用到高纬度地区、海岸带、水体等地区时，不确定性达到了100–140%。As-Syakur等（As-Syakur, A.R.,Tanaka,T.,Prasetia,R.,Swardika,I.K.,and Kasa,I.W.,2011,Comparison of TRMM multisatellite precipitation analysis(TMPA)products and daily-monthly gauge data over Bali.International Journal of Remote Sensing,32,pp.8969-8982.）通过分析巴厘岛地面观测站点数据与TRMM卫星数据之间的关系，发现在干季时，TRMM卫星数据与地面感测数据具有较好的相关性，而这种相关性在湿季较差。 

因此，需要考虑地形、经纬度和距海岸线位置等因素，针对不同的降雨量级别，分别建立不同的预测模型，从而得到更准确地降雨量预测值。这种新的思路需要引进DEM、经纬度信息等，针对不同级别降雨量，对TRMM数据的误差进行分别分析，这就必须引入地统计及空间分析相关方法。 

回归克里格（regression kring）方法是一种结合了回归模型和残差空间插值的混合空间建模方法，以常规的多元回归分析技术和相关的辅助预测数据（如DEM）来分离趋势项，然后对残差进行克里格插值，最后将多元回归预测的趋势项和残差的克里格估计值相加进行估测，回归克里格方法中的趋势项是辅助预测数据的函数。该方法既分析主要影响因素也考虑随机因素，既模拟其空间分布趋势也模拟不确定性，不仅能提高模型的空间精度，预测精度也要优于其他模型。Lloyd等（Lloyd,C.D.2005.Assessing the effect of integrating elevation data into the estimation of monthly precipitation in Great Britain.Journal of Hydrology.308.pp.128-150.）将回归克里格方法用于地面站点观测降雨量的制图；Yemefack等（Yemefack,M.,D.G.Rossiter,and R.Njomgang.2005.Multi-scale characterization of soil variability within an agricultural landscape mosaic system in southern Cameroon.Geoderma.125.pp.117-143.）利用回归克里格方法进行土壤属性预测。 

发明内容

本发明的目的在于提供一种基于卫星遥感与回归克里格的地面日降雨量预测方法，结合降雨分级，分析不同模型中残差的空间趋势，建立不同的预测模型来估算地面降雨量的空间分布特征。 

本发明采用的技术方案的步骤如下： 

步骤（1）数据获取：获取TRMM气象卫星遥感影像数据及ASTERGDEM卫星遥感影像数据，同时收集地面观测站点的日降雨量观测值； 

步骤（2）数据预处理：将步骤（1）获取的TRMM气象卫星遥感影像数据的时间分辨率处理为与地面观测站点的日降雨量观测值的时间分辨率相等，同时将ASTER GDEM卫星遥感影像数据的空间分辨率重采样为1km，并将数据以文件形式存放； 

步骤（3）卫星遥感影像数据提取：根据地面观测站点的空间位置信息，分别提取步骤（2）中相应的TRMM气象卫星遥感影像数据及ASTER GDEM卫星遥感影像数据的栅格值，从而得到每个地面观测站点对应的卫星遥感影像数据的栅格值； 

步骤（4）日降雨量数据的分位数分级：根据步骤（1）收集的地面观测站点的日降雨量观测值及步骤（3）中获取的卫星遥感影像数据的栅格值，进行日降雨量数据的分位数分级； 

步骤（5）各级线性回归关系建立：根据步骤（4）的分位数分级，建立地面观测站点的日降雨量观测值、对应TRMM气象卫星遥感影像数据的栅格值、ASTER GDEM卫星遥感影像数据的栅格值以及经纬度坐标之间的线性回归关系，得到各级降雨量的回归估计值与回归残差； 

步骤（6）各级回归残差分析：分析步骤（5）中得到的各级降雨量的回归残差的空间自相关性及空间集聚性； 

步骤（7）各级回归残差处理：对步骤（6）中具有空间集聚性的各级回归残差进行趋势去除，然后根据ASTER GDEM卫星遥感影像数据对各级回归残差进行Block克里格插值，得到1km空间分辨率的各级回归残差的空间分布特征图； 

步骤（8）各级降雨量预测：将根据步骤（5）得到的各级降雨量的回归估计值与步骤（7）得到的1km空间分辨率的各级回归残差相加，得到1km空间分辨率的各级降雨量预测值； 

步骤（9）地面日降雨量预测：将根据步骤（8）所得的各级降雨量预测值合并，得到地面1km空间分辨率的日降雨量预测值； 

步骤（10）地面日降雨量预测精度评价：利用交叉检验的方法对步骤（9）中的地面日降雨量预测值进行预测精度验证，交叉检验选用均方根误差、平均误差以及偏差作为评价参数。 

步骤（1）中，所述的TRMM气象卫星遥感影像数据空间分辨率为0.25°×0.25°，时间分辨率为3小时，所述的ASTER GDEM卫星遥感影像数据的空间分辨率为30m。 

步骤（2）中，所述的时间分辨率为1天。 

步骤（3）中，所述的两种卫星遥感影像数据的栅格值的提取采用的是点对点的方法，即分别提取最靠近地面站点的TRMM气象卫星遥感影像数据及ASTER GDEM卫星遥感影像数据的栅格值。 

步骤（4）中，所述的日降雨量数据的分位数分级方法为四分位法。 

步骤（5）中，所述的线性回归关系建立采用的是多元线性回归模型，其具体公式如下, 

Z(x₀)＝α+β₁T(x₀)+β₂D(x₀)+β₃Lon(x₀)+β₄Lat(x₀)+ε(x₀) 

其中，Z(x₀)为地面观测站点x₀的降雨回归估计值，α、β₁、β₂、β₃、β₄均是回归系数，T(x₀)是地面观测站点x₀对应的TRMM卫星估测值，D(x₀)是地面观测站点x₀对应的ASTER GDEM卫星遥感影像数据的栅格值，Lon(x₀)为x₀的经度坐标值，Lat(x₀)为x₀的纬度坐标值，ε(x₀)为回归残差。 

步骤（6）中，所述的回归残差的空间自相关性采用半方差模型来判断， 

$\hat{\gamma }\left(h\right)=\frac{1}{2N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{\Sigma }}{[ϵ\left(x_0\right)-ϵ(x_0+h)]}^2$

其中，h是地面点之间的距离，N(h)是样本对数，ε(x₀)和ε(x₀+h)是地面点x₀和x₀+h的回归残差估计值。 

步骤（6）中，所述的回归残差的空间聚集性利用Moran’I系数来判断， 

$I=\frac{n}{S_0}\frac{\underset{i}{\Sigma }\underset{j}{\Sigma }{w}_{\mathrm{ij}}(x_i-\bar{x})(x_j-\bar{x})}{\underset{i}{\Sigma }{(x_i-\bar{x})}^2}$

其中，I为Moran’I系数，x_i为变量在单元i处的值，x_j为变量在单元j处的值，为变量x的平均值，w_ij为空间权重函数，n是由i、j决定的空间变量数。 

步骤（7）中，所述的各级回归残差进行趋势去除，是采用二阶趋势去除法去除。 

步骤（8）中，所述的各级降雨量预测采用以下公式计算， 

$\hat{Z}\left(x_0\right)=\hat{\alpha }\text{+}{\hat{\beta }}_{1}T\left(x_0\right)+{\hat{\beta }}_{2}D\left(x_0\right)+{\hat{\beta }}_3\mathrm{Lon}\left(x_0\right)+{\hat{\beta }}_4\mathrm{Lat}\left(x_0\right)+\widehat{ϵ}\left(x_0\right)$

其中，是各级降雨量的预测值，均是回归估计系数，T(x₀)是地面观测站点x₀对应的TRMM卫星估测值，D(x₀)是地面观测站点x₀对应的ASTER GDEM卫星遥感影像数据的栅格值，Lon(x₀)为x₀的经度坐标值，Lat(x₀)为x₀的纬度坐标值，为回归残差的估计值。 

与背景技术相比，本发明具有的有益效果是： 

本发明采用卫星遥感技术结合回归克里格方法，并利用四分位法进行降雨量分级设计预测模型，可以快速准确地获得地面的日降雨预测空间分布图，不仅大大节省预测时间，而且提高了预测的时空精度，对地面降雨观测站点较少地区的降雨预测尤为重要。解决了水文水利部门进行地面日降雨预报的重要问题，具有重要的理论、实践意义和推广应用价值。 

附图说明

图1是本实施方式的地面观测站点的日均降雨量空间分布特征图。 

图2是本实施方式的TRMM卫星估测的日均降雨量空间分布特征图。 

图3是本实施方式的回归残差的半方差图。 

图4是本实施方式的卫星遥感结合回归克里格方法的某日降雨量的空间分布特征图。 

图5是某一天的地面预测降雨量空间分布图。 

具体实施方式

下面结合附图和实施例对本发明作进一步说明。 

选取浙江省作为研究区域，对2011-2013年雨季（每年5月-10月）的日降雨量进行预测研究，最终得到某日1km空间分辨率的降雨量预测值。 

（1）数据获取：获取浙江省范围内30m空间分辨率的ASTERGDEM卫星遥感影像数据及2011年-2013年雨季（每年5月-10月）的TRMM气象卫星遥感影像数据，同时收集同期内浙江省1379个地面观测站点的日降雨量观测值。选择最靠近TRMM气象卫星遥感影像数据的栅格中心点的地面观测站点为验证站点，共153个；其余1226个地面观测站点为建模站点。地面观测站点的日均降雨量如图1所示。 

（2）数据预处理：将时间分辨率为3小时，空间分辨率为0.25° ×0.25°的TRMM气象卫星遥感影像数据计算为日降雨估测值，同时将30m空间分辨率的ASTER GDEM卫星遥感影像数据重采样为1km，并将计算结果以文件的形式存放。 

$T\left(d\right)=\underset{i=1}{\overset{n}{\Sigma }}T\left(3h\right)$

其中，T(d)为TRMM卫星的日降雨量，T(3h)为TRMM卫星的3小时降雨量，n等于8。 

（3）卫星遥感影像数据提取：根据地面观测站点的空间位置信息，分别提取步骤（2）中相应TRMM气象卫星遥感影像数据及ASTERGDEM卫星遥感影像数据的栅格值，从而得到1379个地面观测站点分别对应的卫星遥感影像数据的栅格值。TRMM卫星估测的日均降雨量如图2所示。 

（4）日降雨量数据的分位数分级：根据步骤（1）收集的地面观测站点的日降雨量观测值及步骤（3）中获取的卫星遥感影像数据的栅格值，通过四分位法，进行日降雨量数据的分位数分级。 

（5）各级线性回归关系建立：根据步骤（4）中日降雨量的四分位数分级，建立1226个建模地面观测站点的日降雨量观测值、对应TRMM气象卫星遥感影像数据的栅格值、ASTER GDEM卫星遥感影像数据的栅格值以及各站点的经纬度坐标之间的线性回归关系，得到各级降雨量的回归估计值与回归残差。 

Z(x₀)＝α+β₁T(x₀)+β₂D(x₀)+β₃Lon(x₀)+β₄Lat(x₀)+ε(x₀) 

其中，Z(x₀)为地面观测站点x₀的降雨回归估计值，α、β₁、β₂、β₃、β₄均是回归系数，T(x₀)是地面观测站点x₀对应的TRMM卫星估测值，D(x₀)是地面观测站点x₀对应的ASTER GDEM卫星遥感影像数据的栅格值，Lon(x₀)为x₀的经度坐标值，Lat(x₀)为x₀的纬度坐标值，ε(x₀)为回归残差。 

（6）各级回归残差分析：利用半方差模型分析步骤（5）中得到的各级降雨量回归残差的空间自相关性，并利用Moran’I系数分析各级降雨量回归残差的空间集聚程度。回归残差的分析结果如图3所示。 

$\hat{\gamma }\left(h\right)=\frac{1}{2N\left(h\right)}\underset{i=1}{\overset{N\left(h\right)}{\Sigma }}{[ϵ\left(x_0\right)-ϵ(x_0+h)]}^2$

其中，h是地面点之间的距离，N(h)是样本对数，ε(x₀)和ε(x₀+h)是地面点x₀和x₀+h的回归残差估计值。 

$I=\frac{n}{S_0}\frac{\underset{i}{\Sigma }\underset{j}{\Sigma }{w}_{\mathrm{ij}}(x_i-\bar{x})(x_j-\bar{x})}{\underset{i}{\Sigma }{(x_i-\bar{x})}^2}$

其中，I为Moran’I系数，x_i为变量在单元i处的值，x_j为变量在单元j处的值，为变量x的平均值，w_ij为空间权重函数，n是由i、j决定的空间变量数。 

（7）各级回归残差处理：采用二阶趋势去除法对步骤（6）中具有空间集聚性的各级回归残差进行趋势去除，然后根据ASTER GDEM卫星遥感影像数据对各级回归残差进行Block克里格插值，得到1km空间分辨率的各级回归残差的空间分布特征图。残差趋势去除后的结果如图4所示。 

（8）各级降雨量预测：将根据步骤（5）得到的各级降雨量的回归估计值与步骤（7）得到的1km空间分辨率的各级回归残差进行空间叠加，得到1km空间分辨率的各级降雨量预测值。 

$\hat{Z}\left(x_0\right)=\hat{\alpha }\text{+}{\hat{\beta }}_{1}T\left(x_0\right)+{\hat{\beta }}_{2}D\left(x_0\right)+{\hat{\beta }}_3\mathrm{Lon}\left(x_0\right)+{\hat{\beta }}_4\mathrm{Lat}\left(x_0\right)+\widehat{ϵ}\left(x_0\right)$

其中，是各级降雨量的预测值，均是回归估计系数，T(x₀)是地面观测站点x₀对应的TRMM卫星估测值，D(x₀)是地面观测站点x₀对应的ASTER GDEM卫星遥感影像数据的栅格值，Lon(x₀)为x₀的经度坐标值，Lat(x₀)为x₀的纬度坐标值，为回归残差的估计值。 

（9）地面日降雨量预测：将根据步骤（8）所得的各级降雨量预测值进行空间合并，得到某一天的1km空间分辨率的日降雨量预测值。某一天的地面预测降雨量空间分布如图5所示。 

（10）地面日降雨量预测精度评价：利用交叉检验的方法结合153个地面验证站点对模型的预测精度进行验证。选用均方根误差（RMSE）、平均误差（ME）以及偏差（Bias）作为评价参数。 

$\mathrm{RMSE}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{(Z_G-Z_S)}^2}$

$\mathrm{MS}=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}(Z_G-Z_S)$

$\mathrm{Bias}=\frac{\underset{i=1}{\overset{n}{\Sigma }}{Z}_S}{\underset{i=1}{\overset{n}{\Sigma }}{Z}_G}-1$

式中，Z_G是地面观测站点的观测降雨量,Z_S是模型预测降雨量,n为样本数。 

精度评价结果如表1所示： 

表1回归克里格方法预测结果与单纯TRMM卫星估计精度比较 

由表1可知，卫星遥感技术和回归克里格方法，结合四分位法进行降雨量分级设计预测模型，具有较高的预测精度。三个检验指标都有不同程度的提高。