技术领域
[0001] 本
发明属于多源遥感数据融合领域，尤其是一种面向田块尺度作物生长监测的遥感影像时空融合方法。
背景技术
[0002] 基于卫星遥感技术的田块尺度作物生长监测对于有效实施精确作物管理至关重要。然而，目前仍缺乏既具备高空间
分辨率、又具备高
时间分辨率的遥感影像，以满足田块尺度作物生长动态监测的需求。一方面，由于复杂的地理条件以及作物生产规模普遍偏小等原因，我国(尤其是长江中下游)的田块
破碎化程度高，田块面积较小，导致用中低空间分辨率遥感影像难以识别不同田块，不同田块间的作物类型、种植方式、生育期等可能存在较大差别；另一方面，受重访周期、
云雨
气候、影像购买成本等影响，在作物生长期内能够获取到的有效高空间分辨率影像有限，难以实现作物长势的连续动态监测。因此，为满足大范围、高
精度、快速变化的田间作物长势监测对高时空分辨率遥感数据的需求，解决目前遥感数据获取能
力不足的问题，迫切需要发展一种面向田块尺度作物生长监测的遥感影像时空融合
算法。
[0003] 面向田间作物监测，必须考虑作物种植区的空间异质性。虽然已有相关研究提出了适用于异质地表环境的时空融合算法，但这些方法均假设所有混合像元中同一组分在多源影像
覆盖时段的
光谱变化率是一致的，即仅考虑了不同类别间(如小麦与其他地表覆盖类别间)的异质性。然而，由于作物所处生育期、生产管理方式、品种等因素的差异，即使同一类作物的生长速率也必然存在差异。若假设属于同一作物所有像元的光谱变化率相同，则削弱了作物的光谱差异和空间变异。
[0004] 时空影像融合误差的主要来源是多尺度遥感影像的空间匹配误差。其中，几何配准误差对
像素级算法的融合精度影响甚大，而绝大多数时空融合模型都属于像素级算法。现有算法通常在遥感数据预处理过程中自动或手动进行几何配准，但并不能完全消除几何配准误差。另外，虽然
传感器的点扩散函数效应也会导致多源遥感影像之间的空间匹配误差，现有方法往往忽略了点扩散函数对时空融合精度的影响。已有研究表明，即使微小的空间匹配误差也会在融合过程中产生显著的偏差。
发明内容
[0005] 本发明所解决的技术问题在于提供一种面向田块尺度作物生长监测的遥感影像时空融合方法，解决构建过程中存在的技术问题，最终得到一组适用于小生产规模作物种植区的高时空分辨率卫星影像数据集，为实现区域性精细田块尺度下的作物(以小麦为例)动态长势监测提供数据支持。
[0006] 实现本发明目的的技术解决方案为：
[0007] 一种面向田块尺度作物生长监测的遥感影像时空融合方法，具体包括以下步骤：
[0008] 步骤1：获取一对准同步的中高分辨率影像及一景获取自作物生长季内且不同于中高分辨率影像对的获取时间的待融合中分辨率影像，并对上述影像进行预处理，其中，准同步指的是获取的时间差不超过5天；
[0009] 步骤2：基于预处理后的中高分辨率影像对，利用二维高斯函数和粒子群
优化算法，构建空间尺度匹配滤波函数，并应用于两景中分辨率影像，两景中分辨率影像分别是与高分辨率影像获取时间准同步的中分辨率影像以及待融合中分辨率影像；
[0010] 步骤3：基于滤波处理后的两景中分辨率影像，以像元为处理单元，计算中分辨率影像获取时间间隔内的光谱变化速率；
[0011] 步骤4：将高分辨率影像逐像元乘以对应的光谱变化速率，生成待融合中分辨率影像对应时间的高分辨率影像。
[0012] 进一步的，本发明的面向田块尺度作物生长监测的遥感影像时空融合方法，步骤1中预处理具体包括：
[0013] 步骤1-1：对非反射率产品进行
辐射定标、大气校正、几何校正，获得反射率数据；
[0014] 步骤1-2：将中、高分辨率影像重投影到相同的
坐标系；
[0015] 步骤1-3：利用最近邻插值法将统一坐标系后的中分辨率影像数据重
采样至高分辨率，得到与高分辨率影像一致的空间分辨率影像数据。
[0016] 进一步的，本发明的面向田块尺度作物生长监测的遥感影像时空融合方法，步骤2中构建空间尺度匹配滤波具体包括：
[0017] 步骤2-1：利用二维高斯函数模拟高空间分辨率传感器的系统点扩散函数；
[0018] 步骤2-2：通过旋转并移动点扩散函数，建立空间匹配滤波；
[0019] 步骤2-3：利用
粒子群优化算法，优化空间匹配滤波参数，得到最优空间匹配滤波。
[0020] 进一步的，本发明的面向田块尺度作物生长监测的遥感影像时空融合方法，步骤2-1中的点扩散函数的公式为：
[0021]
[0022]
[0023]
[0024]
[0025]
[0026] 其中，x、y分别表示东-西方向和南-北方向，FWHM为高斯函数的半峰全宽，X、Y分别为东-西方向和南-北方向上点扩散函数宽度的一半，σx、σy分别为东-西方向和南-北方向上的标准差，G为高斯函数，PSF为点扩散函数。
[0027] 进一步的，本发明的面向田块尺度作物生长监测的遥感影像时空融合方法，步骤2-2中的点扩散函数旋转的公式为：
[0028]
[0029] 其中，θ为点扩散函数以南-北为轴逆
时针旋转的
角度；
[0030] 位移过程是将SMF(FWHMx,FWHMy,θ)向东-西、南-北方向分别平移shiftx和shifty的距离。
[0031] 进一步的，本发明的面向田块尺度作物生长监测的遥感影像时空融合方法，步骤2-3中的粒子群优化算法的代价函数为：
[0032]
[0033]
[0034] 其中，RMSE为均方根误差， 为滤波处理后的中分辨率影像，bn代表第n个波段，t0为参考时间，F(bn,t0)为获取自参考时间t0、波段bn的高分辨率影像反射率数据，N为像元总数，C(bn,t0)为获取自参考时间t0、波段bn的经过步骤1预处理后的中分辨率影像反射率数据，SMF(bn)为波段bn的空间匹配滤波函数， 代表卷积运算。
[0035] 进一步的，本发明的面向田块尺度作物生长监测的遥感影像时空融合方法，步骤3中光谱变化速率的计算公式为：
[0036]
[0037] 其中，MMF(bn)为波段bn中各像元的光谱变化速率， 为滤波处理后的中分辨率影像，t1为待预测时间，t0为参考时间。
[0038] 进一步的，本发明的面向田块尺度作物生长监测的遥感影像时空融合方法，步骤4中生成高时空分辨率影像的公式为：
[0039] F(bn,t1)＝MMF(bn)×F(bn,t0)
[0040] 其中，F(bn,t1)为待预测时间t1、波段bn的经过预处理后的高分辨率影像。
[0041] 本发明采用以上技术方案与
现有技术相比，具有以下技术效果：
[0042] 1、本发明的面向田块尺度作物生长监测的遥感影像时空融合方法同时考虑了类间(如小麦与其他地表覆盖类别间)以及类内(如小麦种植区内)异质性，可重构更为准确的田块空间细节；
[0043] 2、本发明的面向田块尺度作物生长监测的遥感影像时空融合方法可降低多尺度遥感影像间的空间匹配误差，且具有更强的
稳定性；
[0044] 3、本发明的面向田块尺度作物生长监测的遥感影像时空融合方法操作步骤简单、运算速度快、融合精度高，适用于生产规模小的作物种植区，为实现区域性精细田块尺度下的小麦动态长势监测提供数据支持。
附图说明
[0045] 图1是本发明的面向田块尺度作物生长监测的中高分辨率遥感影像时空融合方法的示例影像，其中(a)、(b)是GF-1影像，分别获取自2015年3月12日、4月22日，(c)、(d)是Landsat-8影像，分别获取自2015年3月10日、4月27日，已通过最邻近插值法将Landsat-8影像重采样至与GF-1影像一致的空间分辨率。
[0046] 图2是本发明的面向田块尺度作物生长监测的中高分辨率遥感影像时空融合方法的示意图；
[0047] 图3是其他方法与本发明方法所生成的高分辨率影像与原始影像对比图，其中，(a)为空间和时间自适应反射融合模型STARFM的融合结果，(b)为灵活时空数据融合模型FSDAF的融合结果，(c)为Fit-FC算法的融合结果，(d)为本发明方法的融合结果，(e)为用于结果验证的原始影像，第二行(f-j)和第三行(k-o)分别表示在a-e
中子区域右方框和左方框的放大图，其中虚线椭圆表示融合图像中的预测异常区域；
[0048] 图4是其他方法与本发明方法所生成的高分辨率影像反射率与原始影像各波段反射率关系的
密度散点图，其中，(a-d)为蓝光波段，(e-h)为绿光波段，(i-l)为红光波段，(m-p)为
近红外光波段，从左到右的四列分别对应STARFM、FSDAF、Fit-FC、本发明方法。
具体实施方式
[0049] 下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。
[0050] 本发明
实施例以小麦种植区为例，研究区如图1所示，所采用的卫星数据为高分一号全色/多光谱(PMS/GF-1)影像与准同步的陆地卫星-8(OLI/Landsat-8)影像所对应的可见光(蓝、绿、红波段)和近红外波段，具体波段信息如表1所示。
[0051] 表1 PMS/GF-1波段和相应的OLI/Landsat-8波段信息
[0052]
[0053] 如图2所示，一种面向田块尺度作物生长监测的遥感影像时空融合方法，具体包括以下步骤：
[0054] 步骤1：获取一对准同步(获取时间差不超过5天)的中高分辨率影像及一景获取自作物生长季内(不同于中高分辨率影像对的获取时间)的待融合中分辨率影像，并进行预处理，具体步骤包括：
[0055] 步骤1-1：对非反射率产品进行必要的辐射定标、大气校正、几何校正，获得反射率数据；
[0056] Landsat-8 OLI Level-2数据的表面反射率产品可通过美国地质调查局(USGS)免费获得，其空间分辨率为30米，在本发明实施例中为中分辨率影像数据。对PMS/GF-1数据进行完整的预处理时，包括辐射定标、大气校正、几何校正、波段融合处理，所有操作都在ENVI 5.2中进行。使用FLAASH模型对辐射定标后的数据进行大气校正，得到反射率数据；使用有理多项式系数(RPC)正射校正和相关的数字高程模型(DEM)文件对反射图像进行几何校正；
使用最近邻扩散算法(NNDiffuse)对全色和多光谱图像进行波段融合，得到2米分辨率的高空间分辨率反射率影像。
[0057] 步骤1-2：将中、高分辨率影像重新投影到相同的坐标系UTM WGS-84；
[0058] 步骤1-3：利用最近邻插值法将步骤1-2中的中分辨率影像数据重采样至高分辨率，得到与高分辨率影像一致的空间分辨率影像数据，如图1所示。
[0059] 步骤2：基于预处理后的中高分辨率影像对，利用二维高斯函数和粒子群优化算法，构建空间尺度匹配滤波，，并应用于两景中分辨率影像(与高分辨率影像获取时间准同步的中分辨率影像以及待融合中分辨率影像)，具体步骤包括：
[0060] 步骤2-1：利用二维高斯函数模拟系统点扩散函数，设置FWHM初始范围为20-180m，其中，点扩散函数的公式为：
[0061]
[0062]
[0063]
[0064]
[0065]
[0066] 其中，x、y分别表示东-西方向和南-北方向，FWHM为高斯函数的半峰全宽，X、Y分别为东-西和南-北方向上点扩散函数宽度的一半，σx、σy分别为东-西方向和南-北方向上的标准差，G为高斯函数，PSF为点扩散函数。
[0067] 步骤2-2：设置旋转角度初始范围为0-90°，其中，点扩散函数旋转的公式为：
[0068]
[0069] 其中，θ为点扩散函数以南-北为轴逆时针旋转的角度。设置平移距离的初始范围为0-50m，将SMF(FWHMx,FWHMy,θ)向东-西、南-北方向分别平移shiftx和shifty的距离，通过旋转并移动点扩散函数，建立空间匹配滤波。
[0070] 步骤2-3：设置粒子群优化算法中最大
迭代次数为100，若50次迭代后代价函数没有显著变化，优化过程则自动停止，输出最优空间匹配滤波参数(FWHMx,FWHMy,θ,shiftx,shifty)，其中代价函数为：
[0071]
[0072]
[0073] 其中，RMSE为均方根误差， 为滤波处理后的中分辨率影像，bn代表第n个波段，t0为参考时间，F(bn,t0)为获取自参考时间、波段bn的高分辨率影像的反射率数据，N为像元总数，C(bn,t0)为获取自参考时间、波段bn的经过步骤1预处理后的中分辨率影像反射率数据，SMF(bn)为波段bn的空间匹配滤波， 代表卷积运算。
[0074] 步骤3：基于滤波处理后的两景中分辨率影像，以像元为处理单元，计算中分辨率影像获取时间间隔内的光谱变化速率，其公式为：
[0075]
[0076] 其中，MMF(bn)为波段bn中各像元的光谱变化速率， 为滤波处理后的中分辨率影像，t1为待预测时间，t0为参考时间。
[0077] 步骤4：将高分辨率影像数据逐像元乘以对应的光谱变化速率，生成待融合中分辨率影像对应时间的高时空分辨率影像，其公式为：
[0078] F(bn,t1)＝MMF(bn)×F(bn,t0)
[0079] 其中，F(bn,t1)为待预测时间、波段bn的经过预处理后的高分辨率影像。
[0080] 下面比较本发明提出的方法与STARFM、FSDAF、Fit-FC算法在研究区内对GF-1和Landsat-8时空融合的表现。
[0081] 现有技术和本发明方法所生成的高分辨率影像与原始影像对比图如图3所示，各方法得到的波段反射率与原始影像各波段反射率关系的密度散点图如图4所示，表2汇总了其他方法与本发明方法对研究区内的各波段反射率预测的统计评估情况。从图3、图4和表2可以看出，本发明方法所生成的时空融合结果从目视角度和统计角度来看，与原始影像最为接近，各波段的融合效果均是最好的，本发明的方法重构田块内部空间细节的能力最强、预测反射率的精度最高、总体性能最稳定。
[0082] 表2其他方法与本发明方法对研究区内的各波段反射率预测的统计评估[0083]
[0084] 其中，RMSE是均方根误差，CC是相关系数，粗体值表示对应波段的最低RMSE值。
[0085] 以上所述仅是本发明的部分实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这些改进应视为本发明的保护范围。
