CN107977944A - 一种估测npp遥感图像数据生成方法 - Google Patents
一种估测npp遥感图像数据生成方法 Download PDFInfo
- Publication number
- CN107977944A CN107977944A CN201711312194.9A CN201711312194A CN107977944A CN 107977944 A CN107977944 A CN 107977944A CN 201711312194 A CN201711312194 A CN 201711312194A CN 107977944 A CN107977944 A CN 107977944A
- Authority
- CN
- China
- Prior art keywords
- npp
- remote sensing
- sensing image
- estimation
- original
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 61
- 241001269238 Data Species 0.000 claims abstract description 104
- 230000009466 transformation Effects 0.000 claims abstract description 22
- 239000011159 matrix material Substances 0.000 claims description 21
- 238000012545 processing Methods 0.000 claims description 18
- 230000008569 process Effects 0.000 claims description 10
- 230000008859 change Effects 0.000 claims description 6
- 230000007704 transition Effects 0.000 claims description 3
- 230000000474 nursing effect Effects 0.000 claims description 2
- 230000009467 reduction Effects 0.000 claims description 2
- 238000000354 decomposition reaction Methods 0.000 description 17
- 238000003672 processing method Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 238000001914 filtration Methods 0.000 description 5
- 238000004364 calculation method Methods 0.000 description 4
- 238000009826 distribution Methods 0.000 description 4
- 239000002028 Biomass Substances 0.000 description 3
- 230000003595 spectral effect Effects 0.000 description 3
- 239000000443 aerosol Substances 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000029553 photosynthesis Effects 0.000 description 2
- 238000010672 photosynthesis Methods 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 238000000714 time series forecasting Methods 0.000 description 2
- 206010003694 Atrophy Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000037444 atrophy Effects 0.000 description 1
- 230000001651 autotrophic effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000002457 bidirectional effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000008602 contraction Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000002592 echocardiography Methods 0.000 description 1
- 239000004459 forage Substances 0.000 description 1
- GVVPGTZRZFNKDS-JXMROGBWSA-N geranyl diphosphate Chemical compound CC(C)=CCC\C(C)=C\CO[P@](O)(=O)OP(O)(O)=O GVVPGTZRZFNKDS-JXMROGBWSA-N 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000035800 maturation Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000004091 panning Methods 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000029058 respiratory gaseous exchange Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000004611 spectroscopical analysis Methods 0.000 description 1
- 230000006641 stabilisation Effects 0.000 description 1
- 238000011105 stabilization Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
本发明提供了一种估测NPP遥感图像数据生成方法，该方法结合小波变换与EKF方法生成高质量的估测NPP遥感图像数据。首先,使用离散小波变换方法减少主要由云和大气变化引起的MOD IS NPP噪声，并使用距离加权方法对干扰原始NPP遥感图像数据进行处理；基于离散小波变换的结果获得EKF的平均值、振幅、相位，生成估测NPP遥感图像数据。该估测NPP遥感图像数据较原始NPP遥感图像数据相比，与实际NPP数据更为接近，具有更为良好的准确性。
Description
技术领域
本发明涉及遥感数据处理领域，具体涉及一种估测NPP遥感图像数据生成方法。
背景技术
卡尔曼最初提出的滤波理论只适用于线性系统,Bucy，Sunahara等人提出并研究了扩展卡尔曼滤波(Extended Kalman Filter，简称EKF)，将卡尔曼滤波理论进一步应用到非线性领域。
NPP净初级生产力(net primary productivity)是指植物在单位时间单位面积上由光合作用产生的有机物质总量中扣除自养呼吸后的剩余部分，是生产者能用于生长、发育和繁殖的能量值，反映了植物固定和转化光合产物的效率，也是生态系统中其他生物成员生存和繁衍的物质基础。
MODIS中分辨率成像光谱仪(moderate-resolution imagingspectroradiometer)是当前世界上新一代“图谱合一”的光学遥感仪器，具有36个离散光谱波段，光谱范围宽，从0.4微米到14.4微米全光谱覆盖。
MOD17是净光合作用8天合成产品是MODIS的一种L4级的产品，尽管MOD17 的优化算法得到改进,但MODIS NPP卫星数据还包括各种噪声组件如云、大气变化和双向反射分布的因素。这些因素会影响作物产量,草料,森林生产等等的监控。因此,减少图像噪声和获得高质量的NPP遥感图像数据是非常重要的。
现有的NPP遥感图像数据处理方法很难处理厚重的云层及多云的影像等产生的噪声，并且预测的下一刻数据对上一时刻数据依赖性较强。
因此，需要一种估测NPP遥感图像数据处理方法，减少由云、大气等因素引起的MODIS NPP噪声，并生成较为准确的NPP遥感图像数据。
发明内容
为了克服现有NPP遥感图像数据处理方法的缺陷，本发明提供一种估测NPP 遥感图像数据处理方法，结合小波变换与EKF方法生成高质量的估测NPP遥感图像数据。
相应的，本发明提供的估测NPP遥感图像数据生成方法，包括以下步骤，
获取目标区域在目标时间段内的原始NPP遥感图像数据；
基于离散小波变换处理所述原始NPP遥感图像数据，获得去噪NPP遥感图像数据；
构建NPP时间序列三重调制余弦函数yk＝μk+αkcos(ωk+Φk)+νk，其中yk为去噪NPP遥感图像数据在k时间的值，νk为k时间的噪声样本，ωk为角频率，αk为振幅，Φk为相位，μk为平均值；
输出估测NPP遥感图像数据。
基于所述NPP时间序列三重调制余弦函数，输出估测NPP遥感图像数据。
优选的实施方式，所述原始NPP遥感图像数据为MOD17数据产品。
优选的实施方式，所述估测NPP遥感图像数据生成方法还包括：
基于时距加权函数NPPi＝α×NPPj1+β×NPPj2处理原始NPP遥感图像数据中第i天的干扰原始NPP遥感图像数据；
其中，NPPi为第i天的估算NPP数据，NPPj1为第j1天的原始NPP遥感图像数据，，NPPj2为第j2天的原始NPP遥感图像数据；j1＜i＜j2，为第一时间距离权重，为第二时间距离权重；j1和j2分别为距第i天最近的且具有清晰的原始NPP遥感图像数据的时间点。
优选的实施方式，所述离散小波变换包括以下步骤：
基于离散小波分解处理所述原始NPP遥感图像数据得出原始小波系数；
基于信号噪声比和均方根误差筛选出用于降噪的最优阈值函数；
基于所述最优阈值函数处理所述原始小波系数得出估计小波系数；
基于所述估计小波系数重构得到去噪NPP遥感图像数据。
优选的实施方式，所述EKF预测包括以下步骤：
构建观察模型yk+1＝H·xk+1+vk；其中，yk+1为k+1时间的去噪NPP遥感图像数据，xk+1设定为一个估计值，表示在k+1时间的先验估计，定义为 [μk+1,αk+1,Φk+1]T；H为观测算子；vk～N(O,R)为零均值高斯观测噪声；
预测阶段：
预测状态估计：
预测协方差估计：
其中，F为状态转换矩阵，Pk+1为先验估计误差协方差矩阵，为一种后验估计，为误差协方差矩阵，wk～N(O,Q)为零均值高斯过程噪声，为雅可比矩阵的函数F在关于x的评估
更新阶段：
更新状态估计
近似最佳卡尔曼增益
更新协方差估计
其中，Kg(k)是k时刻卡尔曼增益矩阵，是xk时刻函数H评估的雅克比矩阵；μk,αk,Φk的值为时间的函数，先验初始状态参数[μ1,α1,Φ1]分别为
N为遥感像元素，NPPmax和NPPmin分别为目标时间段内最大和最小NPP；
筛选出初始参数R,P1,Q1,xk+1。
优选的实施方式，所述输出估测NPP遥感图像数据包括所述输出估测NPP 遥感图像数据包括基于初始参数R、P1、Q1、xk+1迭代，输出不同时间点k的 NPP遥感图像数据。
本发明提供的估测NPP遥感图像数据生成方法，小波变换用于减少由云和大气变化引起的MODIS NPP噪声，基于小波变换获得EKF方法的平均值、振幅、相位，从而生成高质量的估测NPP遥感图像数据时间序列预测数据；与现有处理方法相比，本发明提供的NPP时间序列数据处理方法所生成的NPP时间序列预测数据更贴近实际数据，具有更为良好的准确性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。
图1为本发明实施例的估测NPP遥感图像数据生成方法流程图；
图2为本发明实施例的实地调查数据表；
图3为本发明实施例的目标区域的一幅原始NPP遥感图像数据图；
图4为离散小波变换的方法流程图；
图5为离散小波变换的原理图；
图6为本发明实施例基于离散小波变换处理原始NPP遥感图像数据过程示意图；
图7为本发明实施例阈值函数的SNR和RMSE的汇总表；
图8为本发明实施例的去噪前目标区域的局部放大图；
图9为本发明实施例的去噪后目标区域的局部放大图；
图10为本发明实施例的原始NPP遥感图像数据、去噪NPP遥感图像数据和估测NPP遥感图像数据不同天数时的NPP数据折线图；
图11为本发明实施例的估测NPP遥感图像数据与实际NPP数据的误差折线图；
图12为本发明实施例的一组估测NPP遥感图像数据。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
本发明实施例中，提供一种估测NPP遥感图像数据生成方法，该方法结合小波变换与EKF方法生成高质量的估测NPP遥感图像数据。首先,使用离散小波变换方法减少主要由云和大气变化引起的MODIS NPP噪声，并使用距离加权方法对干扰原始NPP遥感图像数据进行处理，获得去噪NPP遥感图像数据；基于去噪NPP遥感图像数据和EKF获NPP时间序列三重调制余弦函数的平均值、振幅、相位，生成估测NPP遥感图像数据。该估测NPP遥感图像数据与原始 NPP遥感图像数据相比，更为接近实际NPP数据，具有较优的准确性。
图2示出了本发明实施例的实地调查数据表。为了验证本发明实施例提供的估测NPP遥感图像数据生成方法准确性，本发明实施例以华南地区作为施行对象。具体实施中，为了提供一个参考标准，在某一时间点，随机采集10个样本地块数据，地块规模大小为1km*1km，测量了每个样本地块内的树种、树高和胸径、和龄级(年轻、中年、成熟)。利用树种的经验回归模型预测估算地上树的生物量，一般形式如下：Biomass＝a×Db，D代表树胸径、a和b是经验系数,因树种而异，具体可参照图2所示的经验列表。
通过每个样本地块的树木生物量值叠加得到估测的地上生物量。使用地块估测生物量除以每个地块平均树龄(年)方法可得每个样本地块的估测地上 NPP值。1km2的地上NPP是五个样本NPP的平均值。每个地块树的平均年龄是由龄级中间值和龄级里面的树木数量加权平均计算所得。按照该方法，对研究时间段内随机时间点施行，以获得实际NPP数据，用于本发明实施例的实际参照和评判标准。
本发明实施例的估测NPP遥感图像数据生成方法，包括以下步骤：
S101：获取目标区域在目标时间段内的原始NPP遥感图像数据；
本发明实施例的原始NPP遥感图像数据来源为8天合成的MOD17数据； MODIS是Terra和Aqua卫星上搭载的主要传感器，两颗卫星相互配合每1～2 天可重复观测整个地球表面，得到36个波段的观测数据；其中MOD17数据为陆地4级标准数据产品，内容为植被产品，GPP/NPP，空间分辨率1公里，频率为8天。
原始NPP遥感图像数据包括多组按时间排序的基于遥感获取的NPP数据，每一组NPP数据包括像元坐标、像元对应的NPP数据等数据信息；其中，每一组数据经计算机软件，如ArcGIS软件读取处理后，以遥感图像的形式进行展示。具体实施中，每一组NPP数据可表示为一个N×N二维矩阵[yij]＝xij+nij，其中，[i,j]表示NPP数据中的一个像元坐标，i＝1…N,j＝1…N,nij为噪声，由于噪声存在是不可估计和预测的，基于统计基础，噪声均值为零，具体实施中，可采用独立同分布高斯白噪声进行近似替代，分布方式和正太分布N(0,σ2)等同，xij为无噪声的光谱图像数据；即[yij]为无噪声NPP遥感图像数据和噪声遥感图像数据的叠加。
需要说明的是，本发明实施例所提供的NPP时间序列数据处理方法其针对对象为每一组NPP数据里的每一个像元，其预测对象为同一坐标的像元在时域上的变化，其整体表现为不同时间节点(8天合成)的NPP遥感图像数据的变化。
图3示出了本发明实施例的原始NPP时间序列数据的其中一组遥感图像数据的遥感图像，时间为2014年1月1日，分辨率为1km，目标区域为华南地区，占地大约65536平方公里(经纬度表示为：N22°36′30″——N24°44′ 0″和E112°53′21.25″——E113°22′13.28″)；目标时间段为2013年1 月9日至2015年12月18日；在该地区的主要植被类型是亚热带常绿阔叶林和常见的农作物；NPP范围从0至7。
S102：基于时距加权函数NPPi＝α×NPPj1+β×NPPj2处理原始NPP遥感图像数据中的干扰原始NPP遥感图像数据；其中，NPPi为第i天的估算NPP数据，NPPj1为具有清晰数据的且离第i天最近的第j1天的原始NPP遥感图像数据，NPPj2为具有清晰数据的且离第j2天最近的第j2天NPP观测数据，j1＜i＜j2，为第一时间距离权重，为第二时间距离权重；
原始NPP遥感图像数据中，其中一部分受多云、气溶胶等因素影响，NPP 数据图像大部分失真，已不具有参考价值，不能用于生成NPP时间序列数据，该类遥感数据称为干扰原始NPP遥感图像数据。具体可将云层覆盖区域连片面积达所述目标区域面积30％以上的遥感图像数据归类为干扰原始NPP遥感图像数据；云层覆盖区域连片是指在该连片区域内，任意一像元和其至少一个相邻像元均为云层覆盖区域，不具有有效的NPP数据；具体实施中，为了统计的便利性，通常以田字形的四个像元均为云层覆盖区域作为云层覆盖区域连片面积的计算单元标准，即在该连片区域内，任意呈田字型排布的四个像元均为云层覆盖区域纳入云层覆盖区域连片面积的计算，最终计算得出的云层覆盖区域连片面积如果占目标区域面积30％以上时，该时刻的原始NPP遥感图像数据为干扰原始NPP遥感图像数据。为了提高NPP数据的仿真精度，获得较为准确的原始NPP遥感图像数据，应舍弃该类图像数据，并基于时距加权函数生成替换数据。
假设第i天的NPP数据为干扰原始NPP遥感图像数据，通过检索查找出位于第i天之前的、且离第i天天数最少的、具有清晰NPP数据图像的NPP遥感图像数据，该NPP遥感图像数据为第j1天；查找出位于第i天之后的、且离第i天天数最少的、具有清晰NPP数据图像的NPP遥感图像数据，该NPP遥感图像数据为第j2天；基于时距加权函数NPPi＝α×NPPj1+β×NPPj2求出第i天的估算NPP数据，其中为第一时间距离权重，为第二时间距离权重。
经S102步骤处理后，原始NPP遥感图像数据被初步加工，经处理后的所有数据可用于参与后续步骤的计算；以下步骤所指的原始NPP遥感图像数据均指经步骤S102处理后的数据。
S103：基于离散小波分解处理所述原始NPP遥感图像数据得出原始小波系数；
小波变换的主要思想为把信号分解成原始小波经过移位和缩放之后的一系列小波，由这些小波来重构原始信号，因此小波是小波变换的基函数，小波系数反映的是不同缩放尺度和平移参数的小波基函数在重构原函数时所占的比重。
小波变换主要分为两大类，一类为离散小波变换，一类为连续小波变换。两者的主要区别在于，连续小波变换在所有可能的缩放和平移上操作，而离散小波变换采用所有缩放和平移值的特定子集。离散小波变换常用于信号编码而连续小波变换常用于信号分析。在许多领域中，离散小波变换可以替代傅里叶变换的作用。
图4示出了离散小波变换的方法流程示意图，图5示出了离散小波变换的具体步骤。离散小波变换主要是利用其特有的多分辨率性、去相关性和选基灵活性特点，使其在图像去噪领域中能发挥重大作用。离散小波图像去噪的实质是，将含噪图像进行多尺度小波分解，从时域变换到小波域，然后在各尺度下尽可能的提取或增强目标图像(频率)的小波系数，并除去噪声的小波系数，最后用小波逆变换重构信号，达到提取目标信号并去噪的目的。这种方法保留了大部分包含信号的小波系数,因此可以较好地保持图象细节。
离散小波变换进行图像去噪主要包括有以下步骤:
S201：输入含噪信号；含噪信号即为本发明实施例的原始NPP遥感图像数据；
S202：对含噪信号进行多尺度分解；多尺度主要包括行分解和列分解，分解后的行分解或列分解可继续进行二级分解、三级等多级分解；
S203：对分解后的含噪信号进行多尺度去噪；可以对每一次分解后的数据进行去噪处理，其中主要涉及到去噪阈值函数的筛选；
S204：基于小波逆变换得出去噪的重构信号；将分解后的数据进行重构，以获取原数据规模的数据；
S205：输出重构信号。
具体实施中，对于二维的原始NPP遥感图像数据，可以用分别在水平和垂直方向使用滤波的方法实现二维小波多分辨率分解，其中L表示低通滤波器， H表示高通滤波器；分解过程可描述为：首先对图像的每一行进行1D-DWT，获得原始图像在水平方向上的低频分量L和高频分量H，然后对变换所得数据的每一列进行1D-DWT，获得原始图像在水平和垂直方向上的低频分量LL、水平方向上的低频和垂直方向上的高频LH、水平方向上的高频和垂直方向上的低频HL以及水平和垂直反向上的高频分量HH。重构过程可描述为：首先对变换结构的每一列进行离散小波逆变换，在对变换所得数据的每一行进行以为离散小波逆变换，即可获得重构图像。具体实施中，图像的小波分解是一个将信号按照低频和有向高频进行分离的过程，分解过程中还可以根据需要对得到的 LL分量进行进一步的小波分解，直至达到要求。
图6示出了本发明实施例的三阶小波分解示意图。在本发明实施例中，根据小波变换进行三级分解去噪，由于高频成分分解更为详细，比一二等级去噪效果更好，另外分解太多等级计算量会随着增大，因此采用三级分解，3个小波等级和3个方向(HH、HL和LH)将被使用。根据分解示意图可知，HL部分波段在水平方向部分波段(右上角部分波段的一半),而那些LH部分波段中在垂直方向部分波段(左下部分波段的一半)，HH部分波段在对角线方向部分波段 (右下部分波段的一半)，LLk部分波段较低频率成分(左上角一半的部分波段)。
S104：筛选最优阈值函数得出估计小波系数；
小波变换具有很强的数据去相关性,能够使信号的能量在小波域集中在少量的小波系数中,而噪声却分布在整个小波域,对应大量的小波系数。经小波分解后,信号的小波系数的幅值要大于噪声,然后就可以用阈值的方法把信号小波系数保留,而使大部分噪声的小波系数减为0。小波域值收缩法去噪的具体处理过程是:将含噪信号在各尺度上进行小波分解,保留大尺度低分辨率下的全部小波系数；对于各尺度高分辨率下的小波系数,可以设定一个阈值,幅值低于该阈值的小波系数全部置0高于该阈值的小波系数或者完,整保留,或者做相应的收缩处理；最后将处理后获得的小波系数利用小波逆变换进行重构,恢复出有效的信号。
小波阈值去噪优点基于高阶统计量的小波阈值去噪方法由于高阶统计量对高斯噪声不敏感,能够排除协方差矩阵高斯白噪声和有色噪声的影响,因而在平滑噪声的同时能更准确地反映原图像的细节信息；利用高阶统计量描述图像的纹理信息对图像进行平滑滤波,可以更好地保留图像细节。缺点小波阈值去噪虽效果较好,但由于将幅值较大的小波系数萎缩会导致图像的边缘模糊, 因此结合小波变换和高阶统计量的特点,利用小波函数和信号相关函数的三重相关系数代替小波系数计算阈值,再通过小波阈值收缩方法对图像进行去噪处理效果会更好一些。
离散小波变换的关键问题归结为寻找最优阈值去除噪声,图像去噪即是从含噪图像去预测无噪图像，使得预测图像与无噪图像的均方误差(MSE)最小，值越小，去噪效果越好。为了获得一个光谱数据之间的估计最小均方误差，在本文中,我们选择六个阈值函数去除图像中的噪声。以2013年1月9日的原始 NPP遥感图像数据为例,通过六个阈值进行处理。
具体实施中，存在多种阈值函数，为了选取最优阈值函数，需要基于信号噪声比SNR和均方根RMSE实现阈值函数筛选；
计算公式如下所示，
其中，Xi为预测值，为真值，n为像元数；
S为信号值，N为噪声值；
本发明实施例列举六个阈值函数供筛选，分别为：
软阈值函数：
硬阈值函数：y＝x*(|x|＞t)；
半软阈值函数：σ＝2，t0＝σ*t，t1＝|x|-t，y0＝sign(x)*t2， y＝x*(|x|＞＝t0)+y0*(|x|＜t0)；
自定义阈值函数一：a＝0.1，σ＝0.4，t0＝σ*t，t1＝|x|-a*t，t2＝|x|-t0，y0＝sign(x)*t4；
自定义阈值函数二：k＝2，y＝sign(x)*t0*(|x|＞＝t)+t1*(|x|＜t)；
自定义阈值函数三：y＝sign(x)*t0*(|x|＞＝t)；
以上公式中，y为经阈值函数处理后的小波系数估计值，x为含噪信号小波变换后的小波系数，t为阈值。
图7示出了阈值函数的SNR和RMSE值表。数据表明，本发明实施例中自定义阈值函数二更适合用于图像中的噪声。本发明实施例中，自定义阈值函数 2使用近似对称的紧支集正交小波(dbN)，消失矩为8，该小波具有较好的正则性，即使用该小波所引入的光滑误差不容易被察觉，使得信号重构过程比较光滑，该自定义阈值函数2的特点是消失矩越高，光滑性就越好，频域的局部化能力就越强，频带的划分效果越好。
S105：基于所述小波系数估计值重构得到去噪NPP遥感图像数据。
本发明实施例经自定义阈值函数二处理后的小波分解图像，基于小波逆变换进行重构，得到去噪NPP遥感图像数据；
图8示出了本发明实施例的一组原始NPP遥感图像数据的遥感图像和局部放大图；图9示出了本发明实施例的一组去噪NPP遥感图像数据的遥感图像和局部放大图。随机抽取原始NPP遥感图像数据中的一组或多组，分别生成原始 NPP遥感图像和经步骤S103～S105处理后的去噪NPP遥感图像，通过视觉对比，去噪后的影像中白点噪声较少，白色斑纹被抑制。原始NPP遥感图像数据中云和气溶胶在大气中似乎更多的噪声，而小波变换的结果平滑且噪声量较少。
需要说明的是，步骤S103～S105针对对象为每一组NPP遥感图像数据，即一个时间点的NPP遥感图像数据。
S106：构建NPP时间序列三重调制余弦函数yk＝μk+αkcos(ωk+Φk)+vk；其中yk代表原始NPP遥感图像数据k时间的值，vk为k时间噪声样本，ωk＝2πf表示角频率，αk表示振幅，Φk表示相位。方程表示基于ωk、μk、αk、Φk的余弦函数。角频率可以计算ωk＝2πf，f是由年度植被生长周期决定的。在这项研究,因为Modis NPP数据从MOD17数据产品获得的，时间分辨率为8天，f设置为8/365。
此外，一个像元的μk、αk、Φk的值是与时间相关的函数，主要依赖于yk，因此，需要使用EKF算法来模拟估测的μk、αk、Φk值。
需要说明的是，该步骤中构建的NPP时间序列三重调制余弦函数中，μk、αk、Φk为未知的，需要通过步骤S107的EKF求出。
S107：使用了EKF算法来模拟估测的μk、αk、Φk值。
扩展卡尔曼滤波器EKF是解决许多估测问题最好的工具,EKF预测公式可以很容易通过线性化过程模型获得。
k+1时的模型观察yk+1＝H·xk+1+vk；
yk+1在k+1时估测的归一化植被指数，近红外波段的反射值与红光波段的反射值之差比上两者之和，xk+1表示在k+1时先验估计，依据经验事先给予一个估计值，定义为xk+1＝[μk+1αk+1Φk+1]T,H是观测算子，用于关联模型状态变量的观察，不需要线性的，vk是零均值高斯观测噪声vk～N(O,R)。离散时间EKF预测算法有两个主要阶段：
预测阶段预测方程
预测状态估计
预测协方差估计
F表示状态转换矩阵，Pk+1是先验估计误差协方差矩阵，代表了一种后验估计，是经验的误差协方差矩阵，Wk是零均值高斯过程噪声Wk～N(O,Q)，是雅可比矩阵的函数F在关于x的评估。
更新状态方程
近似最佳卡尔曼增益
更新协方差估计
Kg(k)是k时刻卡尔曼增益矩阵，是xk时刻函数H评估的雅可比矩阵。μk、αk、Φk的值是时间的函数,必须估计给定yk，k∈1…M。选择137张影像, 所以M＝100。先验初始状态参数(μk、αk、Φk)被设置为：
其中，N是遥感像元素，NPPmax和NPPmin分别为研究区内最大和最小NPP。
选择先验的初始估计Φ1＝0°假设的开始日期是2013年1月18日代表植被交替后的阶段。
本发明实施例中，根据NPP模型和NPP样本的布局，多次调整EKF设置参数，然后，初始参数设置为：R＝0.01,P1＝diag[0.015,0.015,0.02963], xk+1＝[4.5612,2.98425,0]，其中，xk+1包括了μk+1αk+1Φk+1的数据，通过回代至公式和迭代计算，可生成估测NPP遥感图像数据。
S108：生成估测NPP遥感图像数据。
图12示出了本发明实施例的一组估测NPP遥感图像数据，由于实际运算中数据角度，图12只列出了其中一部分。需要说明的是，该组数据代表时间点为1时的同一组估测NPP遥感图像数据中，不同像元的NPP数据，整体的估测NPP遥感图像数据量较为庞大，具体实施中以估测NPP遥感图像进行表示。
图10示出了本发明实施例的原始NPP遥感图像数据、去噪NPP遥感图像数据、估测NPP遥感图像数据对比折线图。为了验证估测NPP遥感图像数据的准确性，将原始NPP遥感图像数据、去噪NPP遥感图像数据、估测NPP遥感图像数据绘制在同一图中，具体比对的为整体的NPP；将估测NPP遥感图像数据与去噪NPP遥感图像数据比较，计算每组NPP数据的绝对误差百分比，绝对误差百分比＝|估测值-实测值|/实测值*100％，统计结果。
图11示出了估测NPP遥感图像数据与原始NPP遥感图像误差百分比折线图。初始设置参数导致误差的变化非常大，由于预先给定的初始值与实测值误差比较大造成的，利用卡尔曼滤波(EKF)算法来模拟估测的三重调制余弦函数中μk、αk、Φk值，随着观测数据增加，估算值变得稳定更加接近真值。20 天以后,误差开始变得小。直到大约120天,误差趋向与0％。这些较低的误差表明本文提出的NPP估测方法在估测精度高的NPP结果比只有EKF的方法有更大的潜力。
结合本发明实施例的估测NPP遥感图像数据生成方法生成估测NPP遥感图像数据，并与原始NPP遥感图像数据一起，分别计算与实际调查样本中每个样本地块绝对误差百分比。在1平方公里内五个样本的平均NPP分别为3.752 和3.716Mg.ha-1.year-1。其中，估测NPP遥感图像数据的绝对误差百分比为0.243％至4.79％；原始NPP遥感图像数据为29.18％至41.26％，本发明提供的估测NPP遥感图像数据与实际数据更为接近，NPP遥感图像数据更为有效。
本发明实施例中，提供一种估测NPP遥感图像数据生成方法，该方法结合小波变换与EKF方法生成高质量的估测NPP遥感图像数据。首先,使用离散小波变换方法减少主要由云和大气变化引起的MODIS NPP噪声，并使用距离加权方法对干扰原始NPP遥感图像数据进行处理；基于离散小波变换的结果获得 EKF的平均值、振幅、相位，生成估测NPP遥感图像数据。估测NPP遥感图像数据与原始NPP遥感图像相比，更为符合实际情况，具有更为良好的准确性。
以上对本发明实施例所提供的一种估测NPP遥感图像数据生成方法进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。
Claims (6)
1.一种估测NPP遥感图像数据生成方法，其特征在于，包括以下步骤，
获取目标区域在目标时间段内的原始NPP遥感图像数据；
基于离散小波变换处理所述原始NPP遥感图像数据，获得去噪NPP遥感图像数据；
构建NPP时间序列三重调制余弦函数yk＝μk+αkcos(ωk+Φk)+νk，其中,yk为去噪NPP遥感图像数据在k时间的值，νk为k时间的噪声样本，ωk为角频率，αk为振幅，Φk为相位,μk为平均值；
基于去噪NPP遥感图像数据，使用扩展卡尔曼滤波得出平均值μk、振幅αk和相位Φk；
输出估测NPP遥感图像数据。
2.如权利要求1所述的估测NPP遥感图像数据生成方法，其特征在于，所述原始NPP遥感图像数据为MOD17数据产品。
3.如权利要求1所述的估测NPP遥感图像数据生成方法，其特征在于，所述估测NPP遥感图像数据生成方法还包括：
基于时距加权函数NPPi＝α×NPPj1+β×NPPj2处理原始NPP遥感图像数据中第i天的干扰原始NPP遥感图像数据；
其中，NPPi为第i天的估算NPP数据，NPPj1为第j1天的原始NPP遥感图像数据，NPPj2为第j2天的原始NPP遥感图像数据；j1＜i＜j2，为第一时间距离权重，为第二时间距离权重；j1和j2分别为距第i天最近的且具有清晰的原始NPP遥感图像数据的时间点。
4.如权利要求1所述的估测NPP遥感图像数据生成方法，其特征在于，所述离散小波变换包括以下步骤：
基于离散小波分解处理所述原始NPP遥感图像数据得出原始小波系数；
基于信号噪声比和均方根误差筛选出用于降噪的最优阈值函数；
基于所述最优阈值函数处理所述原始小波系数得出估计小波系数；
基于所述估计小波系数重构得到去噪NPP遥感图像数据。
5.如权利要求1所述的估测NPP遥感图像数据生成方法，其特征在于，所述基于去噪NPP遥感图像数据，使用扩展卡尔曼滤波得出初始值μk、振幅αk和相位Φk包括以下步骤：
构建观察模型yk+1＝H·xk+1+vk；其中，yk+1为k+1时间的去噪NPP遥感图像数据，xk+1设定为一个估计值，表示在k+1时间的先验估计，定义为[μk+1,αk+1,Φk+1]T；H为观测算子；vk～N(O,R)为零均值高斯观测噪声；
预测阶段：
预测状态估计：
预测协方差估计：
其中，F为状态转换矩阵，Pk+1为先验估计误差协方差矩阵，为一种后验估计，为误差协方差矩阵，wk～N(O,Q)为零均值高斯过程噪声，▽Fx为雅可比矩阵的函数F在关于x的评估
更新阶段：
更新状态估计
近似最佳卡尔曼增益
更新协方差估计
其中，Kg(k)是k时刻卡尔曼增益矩阵，▽Hx是xk时刻函数H评估的雅克比矩阵；μk,αk,Φk的值为时间的函数，先验初始状态参数[μ1,α1,Φ1]分别为
Φ1＝0°；
N为遥感像元素，NPPmax和NPPmin分别为目标时间段内最大和最小NPP；
筛选出初始参数R,P1,Q1,xk+1。
6.如权利要求5所述的估测NPP遥感图像数据生成方法，其特征在于，所述输出估测NPP遥感图像数据包括基于初始参数R、P1、Q1、xk+1迭代，输出不同时间点k的NPP遥感图像数据。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711312194.9A CN107977944A (zh) | 2017-12-12 | 2017-12-12 | 一种估测npp遥感图像数据生成方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711312194.9A CN107977944A (zh) | 2017-12-12 | 2017-12-12 | 一种估测npp遥感图像数据生成方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN107977944A true CN107977944A (zh) | 2018-05-01 |
Family
ID=62009976
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711312194.9A Pending CN107977944A (zh) | 2017-12-12 | 2017-12-12 | 一种估测npp遥感图像数据生成方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107977944A (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109472834A (zh) * | 2018-10-23 | 2019-03-15 | 桂林电子科技大学 | 一种基于小波变换的卡尔曼滤波相位展开方法 |
| CN112037151A (zh) * | 2020-09-08 | 2020-12-04 | 天元大数据信用管理有限公司 | 基于小波分析的图像去噪方法 |
| CN112257225A (zh) * | 2020-09-16 | 2021-01-22 | 中国科学院地理科学与资源研究所 | 一种适用于高寒草地生态系统的npp计算方法 |
| CN113887324A (zh) * | 2021-09-10 | 2022-01-04 | 北京和德宇航技术有限公司 | 基于卫星遥感数据的火点检测方法 |
| CN114445467A (zh) * | 2021-12-21 | 2022-05-06 | 贵州大学 | 基于视觉的四旋翼无人机特定目标识别与追踪系统 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102663695A (zh) * | 2012-03-31 | 2012-09-12 | 重庆大学 | 基于小波变换的dr图像去噪方法及系统 |
| CN103530857A (zh) * | 2013-10-31 | 2014-01-22 | 清华大学深圳研究生院 | 基于多尺度的卡尔曼滤波图像去噪方法 |
| CN106021653A (zh) * | 2016-05-06 | 2016-10-12 | 华南农业大学 | 一种ndvi时间序列重构的方法及系统 |
| CN106026260A (zh) * | 2016-06-24 | 2016-10-12 | 南京航空航天大学 | 一种带有均衡电路的串连电池组soc估算方法 |
| CN107064932A (zh) * | 2017-02-28 | 2017-08-18 | 华南农业大学 | 一种基于时间序列sar遥感影像的建设用地变化检测方法 |
-
2017
- 2017-12-12 CN CN201711312194.9A patent/CN107977944A/zh active Pending
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102663695A (zh) * | 2012-03-31 | 2012-09-12 | 重庆大学 | 基于小波变换的dr图像去噪方法及系统 |
| CN103530857A (zh) * | 2013-10-31 | 2014-01-22 | 清华大学深圳研究生院 | 基于多尺度的卡尔曼滤波图像去噪方法 |
| CN106021653A (zh) * | 2016-05-06 | 2016-10-12 | 华南农业大学 | 一种ndvi时间序列重构的方法及系统 |
| CN106026260A (zh) * | 2016-06-24 | 2016-10-12 | 南京航空航天大学 | 一种带有均衡电路的串连电池组soc估算方法 |
| CN107064932A (zh) * | 2017-02-28 | 2017-08-18 | 华南农业大学 | 一种基于时间序列sar遥感影像的建设用地变化检测方法 |
Non-Patent Citations (5)
| Title |
|---|
| MARCOS A. M. LAIA 等: "A Novel Model for Combining Projection and Image Filtering using Kalman and Discrete Wavelet Transform in Computerized Tomography", 《2008 11TH IEEE INTEMATIONA1 CONFERENCE ON COMPUTATIONA1 SCIENCE AND ENGINEERING》 * |
| PERUMAL ARUMUGA 等: "image denoising using discrete wavelet transform", 《INTERNATIONAL JOURNAL OF COMPUTER SCIENCE AND NETWORK SECURITY》 * |
| SEDA POSTALCIOGLU 等: "comparison of kalman filter and wavelet filter for denoising", 《2005 INTERNATIONAL CONFERENCE ON NEURAL NETWORKS AND BRAIN》 * |
| XIAOLIANG LU 等: ""Removal of Noise by Wavelet Method to Generate High Quality Temporal Data of Terrestrial MODIS Products"", 《PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING》 * |
| 蒋雪冰: ""线性内插和扩展Kalman滤波组合法对NDVI时间序列的重构与预测"", 《中国优秀硕士学位论文全文数据库 基础科学辑》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109472834A (zh) * | 2018-10-23 | 2019-03-15 | 桂林电子科技大学 | 一种基于小波变换的卡尔曼滤波相位展开方法 |
| CN109472834B (zh) * | 2018-10-23 | 2023-04-14 | 桂林电子科技大学 | 一种基于小波变换的卡尔曼滤波相位展开方法 |
| CN112037151A (zh) * | 2020-09-08 | 2020-12-04 | 天元大数据信用管理有限公司 | 基于小波分析的图像去噪方法 |
| CN112257225A (zh) * | 2020-09-16 | 2021-01-22 | 中国科学院地理科学与资源研究所 | 一种适用于高寒草地生态系统的npp计算方法 |
| CN112257225B (zh) * | 2020-09-16 | 2023-07-14 | 中国科学院地理科学与资源研究所 | 一种适用于高寒草地生态系统的npp计算方法 |
| CN113887324A (zh) * | 2021-09-10 | 2022-01-04 | 北京和德宇航技术有限公司 | 基于卫星遥感数据的火点检测方法 |
| CN114445467A (zh) * | 2021-12-21 | 2022-05-06 | 贵州大学 | 基于视觉的四旋翼无人机特定目标识别与追踪系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107977944A (zh) | 一种估测npp遥感图像数据生成方法 | |
| Phung et al. | Monitoring rice growth status in the Mekong Delta, Vietnam using multitemporal Sentinel-1 data | |
| Feilhauer et al. | Mapping the local variability of Natura 2000 habitats with remote sensing | |
| CN103914678B (zh) | 基于纹理与植被指数的撂荒地遥感识别方法 | |
| CN101201937B (zh) | 基于小波重构与分解的数字图像增强方法及其装置 | |
| Pan et al. | Detailed and automated classification of land use/land cover using machine learning algorithms in Google Earth Engine | |
| Qiu et al. | Finer resolution estimation and mapping of mangrove biomass using UAV LiDAR and worldview-2 data | |
| CN109583378A (zh) | 一种植被覆盖度提取方法及系统 | |
| CN100557633C (zh) | 基于梳状波和拉普拉斯塔形分解的多源图像融合方法 | |
| CN112348812B (zh) | 林分年龄信息测量方法及装置 | |
| CN103854267A (zh) | 一种基于变分与分数阶微分的图像融合与超分辨率实现方法 | |
| CN101303764A (zh) | 基于非下采样轮廓波的多传感器图像自适应融合方法 | |
| US20140369568A1 (en) | Detecting Species Diversity by Image Texture Analysis | |
| Wang et al. | Biophysical estimation in tropical forests using JERS‐1 SAR and VNIR imagery. II. Aboveground woody biomass | |
| CN102063715A (zh) | 基于nsct和粒子群优化算法的台风云图融合方法 | |
| CN109063657A (zh) | 面向均值地域光谱单元的地上生物量估算和尺度转换方法 | |
| CN110032939B (zh) | 一种基于高斯混合模型的遥感时序数据拟合方法 | |
| CN114926748A (zh) | Sentinel-1/2微波与光学多光谱影像结合的大豆遥感识别方法 | |
| CN104143031B (zh) | 一种基于小波多尺度分解的植被指数时序数据重构方法 | |
| Foroumandi et al. | Climate change or regional human impacts? Remote sensing tools, artificial neural networks, and wavelet approaches aim to solve the problem | |
| Karna | Mapping above ground carbon using worldview satellite image and lidar data in relationship with tree diversity of forests | |
| Bejarano et al. | Combining optical and acoustic data to enhance the detection of Caribbean forereef habitats | |
| Yang et al. | Spatial-temporal dynamic monitoring of vegetation recovery after the Wenchuan earthquake | |
| Liang et al. | Maximum likelihood classification of soil remote sensing image based on deep learning | |
| CN116452023B (zh) | 基于低频微波雷达vod数据的公里级碳储量评估方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20180501 |
|
| RJ01 | Rejection of invention patent application after publication |