CN114780904B - 一种端元自适应的山地植被覆盖度遥感反演方法 - Google Patents

一种端元自适应的山地植被覆盖度遥感反演方法 Download PDF

Info

Publication number
CN114780904B
CN114780904B CN202210686048.7A CN202210686048A CN114780904B CN 114780904 B CN114780904 B CN 114780904B CN 202210686048 A CN202210686048 A CN 202210686048A CN 114780904 B CN114780904 B CN 114780904B
Authority
CN
China
Prior art keywords
end member
spectral
pixel
adaptive
vegetation coverage
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.)
Active
Application number
CN202210686048.7A
Other languages
English (en)
Other versions
CN114780904A (zh
Inventor
边金虎
李爱农
张正健
雷光斌
南希
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Mountain Hazards and Environment IMHE of CAS
Original Assignee
Institute of Mountain Hazards and Environment IMHE of CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Institute of Mountain Hazards and Environment IMHE of CAS filed Critical Institute of Mountain Hazards and Environment IMHE of CAS
Priority to CN202210686048.7A priority Critical patent/CN114780904B/zh
Publication of CN114780904A publication Critical patent/CN114780904A/zh
Application granted granted Critical
Publication of CN114780904B publication Critical patent/CN114780904B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及山区多光谱遥感卫星影像的计算机处理技术领域，具体而言，涉及一种端元自适应的山地植被覆盖度遥感反演方法。步骤如下：获取端元光谱数据并构建端元光谱库；根据端元光谱库，对待分解像元执行逐像元基于最小光谱欧氏距离的自适应选择；进入迭代过程，从原始像元的光谱反射率中扣除已选择端元的影响，进行下一端元的自适应选择；通过预定义的拟合效果指数判定是否结束迭代，选择指数最小值对应的端元作为自适应选择端元；采用线性混合像元分解模型对自适应选择端元进行分解，得到植被覆盖度的反演结果。本方法能够实现复杂山地时间序列植被覆盖度高精度反演，并生成相应的误差报告，具有大区域、高效率、高精度、低成本的优点。

Description

一种端元自适应的山地植被覆盖度遥感反演方法
技术领域
本发明涉及山区多光谱遥感卫星影像的计算机处理技术领域，具体而言，涉及一种端元自适应的山地植被覆盖度遥感反演方法。
背景技术
植被覆盖度是指单位面积内植被垂直投影面积所占的比例。它不仅是表征地表植被覆盖状态的重要信息，同时还是生态系统退化、沙化过程的重要敏感因子。植被覆盖度是许多陆表模型、水文模型、土壤侵蚀模型和天气预报模型的关键输入变量。获得植被覆盖度的空间分布及其季相变化特征对于分析生态系统的功能及其变化，探索驱动因子以及评估在全球变化和人类活动的双重影响下生态系统的健康状况都至关重要。卫星遥感数据现已成为大区域植被覆盖度获取的主要数据源。由于遥感数据具有长时间序列、大区域和周期性覆盖，许多卫星遥感影像已经被广泛地应用于植被覆盖度的制图和动态模拟中。在过去的数十年里，国内外学者已经提出了多种基于遥感的植被覆盖度估算方法，分类方法、光谱分解法、回归树法、以及基于物理模型的反演方法等。在这些方法中，光谱分解法由于其具有很强的灵活性、较好的物理基础而被广泛采用。
一般来说，光谱分解法假设每一个像元是由几个部分（端元）组成的，因此一个像元可以被分解为几个端元及其相应面积比例的光谱组合。其中，绿色植被的面积比例被认为是植被的覆盖度。然而，在现实中，许多像元的光谱组合往往仅是端元库中的少数端元组成，尤其在复杂山地，地表更加破碎，如图1所示，图1展示了典型复杂山地无人机拍摄获取的高分辨率影像及对应的30米分辨率卫星影像像元格网，可以看出，30米分辨率卫星影像格网中包含了多种地物特征。同时，遥感像元端元的组成同样会随着植被的季相和地表条件变化而发生端元组合的变化。
因此，采用传统的算法原理开展地表覆被类型复杂、高时空异质性的山区多光谱遥感卫星影像拍摄的多光谱遥感卫星影像的植被覆盖度计算，无法准确的估算得到植被覆盖度。
因此，亟需提出一种新的关于山地植被覆盖度的遥感反演方法，以解决现有技术中存在的技术缺陷。
发明内容
本发明的目的在于解决背景技术中所指出的问题，提供一种能够实现高效率、高精度处理的端元自适应的山地植被覆盖度遥感反演方法。
本发明的实施例通过以下技术方案实现：一种端元自适应的山地植被覆盖度遥感反演方法，包括如下步骤：
S1.获取端元光谱数据并构建端元光谱库；
S2.根据端元光谱库，对原始像元中的待分解像元执行逐像元基于最小光谱欧氏距离的自适应选择，得到第一端元；
S3.进入迭代过程，从原始像元的光谱反射率中扣除已选择端元的影响，进行下一端元的自适应选择；
通过预定义的拟合效果指数判定是否结束迭代，选择指数最小值对应的端元作为自适应选择端元；
S4.采用线性混合像元分解模型对所述自适应选择端元进行分解，得到植被覆盖度的反演结果。
根据一种优选实施方式，所述方法在步骤S1之前还包括：定义端元的数量及每个端元的光谱特征。
根据一种优选实施方式，步骤S1还包括：基于传感器的光谱响应函数，对端元光谱数据进行地表反射率的转换，表达式如下：
Figure 966142DEST_PATH_IMAGE001
上式中，
Figure 335944DEST_PATH_IMAGE002
表示传感器在i波段的地表反射率，
Figure 757698DEST_PATH_IMAGE003
表示传感器在i波段的光谱响应函数，
Figure 85911DEST_PATH_IMAGE004
表示地表在波长
Figure 960326DEST_PATH_IMAGE005
处的反射率值，
Figure 399398DEST_PATH_IMAGE006
表示波长的积分，
Figure 675658DEST_PATH_IMAGE007
表示传感器在i波段的光谱上界，
Figure 909194DEST_PATH_IMAGE008
表示传感器在i波段的光谱下界。
根据一种优选实施方式，步骤S2具体包括：
为各待分解像元与端元光谱库中每一个端元计算光谱欧氏距离，表达式如下：
Figure 270905DEST_PATH_IMAGE009
上式中，
Figure 982509DEST_PATH_IMAGE010
表示待分解像元t与端元e的光谱欧氏距离，m表示光谱波段数，
Figure 378855DEST_PATH_IMAGE011
表示待分解像元t在波段b上的光谱反射率，
Figure 48871DEST_PATH_IMAGE012
表示端元e在波段b上的光谱反射率；
选择最小光谱欧氏距离的端元，作为第一端元。
根据一种优选实施方式，步骤S3中所述预定义的拟合效果指数的表达式如下：
Figure 897878DEST_PATH_IMAGE013
上式中，
Figure 413173DEST_PATH_IMAGE014
表示第k个端元的拟合效果指数，m表示光谱波段数，
Figure 664026DEST_PATH_IMAGE015
表示待分解像元t在波段i上的光谱反射率，
Figure 239364DEST_PATH_IMAGE016
表示端元e在波段i上的光谱反射率，
Figure 575667DEST_PATH_IMAGE017
表示端元e光谱反射率归一化后的贡献权重。
根据一种优选实施方式，所述贡献权重的计算表达式如下：
Figure 894653DEST_PATH_IMAGE018
上式中，
Figure 12DEST_PATH_IMAGE010
表示待分解像元t与端元e的光谱欧氏距离。
根据一种优选实施方式，步骤S3中所述通过预定义的拟合效果指数判定是否结束迭代具体包括：
若满足如下表达式，则结束迭代，表达式如下：
Figure 746251DEST_PATH_IMAGE019
根据一种优选实施方式，步骤S4中所述线性混合像元分解模型，表达式如下：
Figure 304272DEST_PATH_IMAGE020
上式中，
Figure 692528DEST_PATH_IMAGE021
表示像元I在波段
Figure 652394DEST_PATH_IMAGE023
的光谱反射率，m表示端元的数量，
Figure 569534DEST_PATH_IMAGE024
表示像元I的第k个端元的面积比例，
Figure 146009DEST_PATH_IMAGE025
表示第k个端元在波段
Figure 806797DEST_PATH_IMAGE023
的光谱反射率，
Figure 621170DEST_PATH_IMAGE026
表示线性混合像元分解模型的拟合残差。
根据一种优选实施方式，所述方法还包括，为所述线性混合像元分解模型设定约束条件，使得所有的混合分解结果之和为1，以及面积比例的非负性，表达式如下：
Figure 974790DEST_PATH_IMAGE027
根据一种优选实施方式，所述方法还包括，对反演结果进行验证，验证公式采用均方根误差RMSE、平均误差AD以及可决系数R2，具体表达式如下：
Figure 507403DEST_PATH_IMAGE028
Figure 706303DEST_PATH_IMAGE029
Figure 640761DEST_PATH_IMAGE030
上式中，N表示待分解像元数，
Figure 165283DEST_PATH_IMAGE031
表示像元I的模型估算值，
Figure 919613DEST_PATH_IMAGE032
表示像元I的地面实测参考值或无人机提取参考值，
Figure 187783DEST_PATH_IMAGE033
表示像元I的参考值均值。
本发明实施例的技术方案至少具有如下优点和有益效果：本发明所提供的山地植被覆盖度遥感反演方法，考虑到复杂山地是由许多复杂的地表覆被类型所构成，通过综合高分辨率卫星影像和新的自适应端元光谱分解模型，能够实现复杂山地时间序列植被覆盖度高精度反演，并生成相应的误差报告，具有大区域、高效率、高精度、低成本的优点；本方法对海量山地多光谱遥感影像的植被覆盖度反演非常有效，解决了传统线性光谱混合分解模型端元固定、反演精度低的问题，对高时空异质性的山地区域具有很高的有效性和可靠性，对于目前在轨的卫星数据反演，对山地植被监测、生态环境健康监测等具有十分重要的意义。
附图说明
图1为典型复杂山地无人机拍摄获取的高分辨率影像及对应的30米分辨率卫星影像像元格网示意图；
图2为线性光谱分解模型原理示意图；
图3为本发明实施例1提供的端元自适应选择植被覆盖度反演方法的流程图；
图4（a）为野外测量获得的草地反射率曲线及结合卫星光谱响应曲线转换后的地表反射率；
图4（b）为野外测量获得的不透水面反射率曲线及结合卫星光谱响应曲线转换后的地表反射率；
图4（c）为野外测量获得的水体反射率曲线及结合卫星光谱响应曲线转换后的地表反射率；
图4（d）为野外测量获得的裸土反射率曲线及结合卫星光谱响应曲线转换后的地表反射率；
图5（a）为生长季旺盛期HJ合成影像标准假彩色合成（RGB=4,3,2）；
图5（b）为无人机植被覆盖度参考数据；
图5（c）为本端元自适应选择方法模型估算结果；
图5（d）为传统线性光谱分解模型估算结果；
图6（a）为时间序列FVC在1月1日至1月16日的估算结果；
图6（b）为时间序列FVC在1月17日至2月1日的估算结果；
图6（c）为时间序列FVC在2月2日至2月17日的估算结果；
图6（d）为时间序列FVC在2月18日至3月5日的估算结果；
图6（e）为时间序列FVC在3月6日至3月21日的估算结果；
图6（f）为时间序列FVC在3月22日至4月6日的估算结果；
图6（g）为时间序列FVC在4月7日至4月22日的估算结果；
图6（h）为时间序列FVC在4月23日至5月8日的估算结果；
图6（i）为时间序列FVC在5月9日至5月24日的估算结果；
图6（j）为时间序列FVC在5月25日至6月9日的估算结果；
图6（k）为时间序列FVC在6月10日至6月25日的估算结果；
图6（l）为时间序列FVC在6月26日至7月11日的估算结果；
图6（m）为时间序列FVC在7月12日至7月27日的估算结果；
图6（n）为时间序列FVC在7月28日至8月12日的估算结果；
图6（o）为时间序列FVC在8月13日至8月28日的估算结果；
图6（p）为时间序列FVC在8月29日至9月13日的估算结果；
图6（q）为时间序列FVC在9月14日至9月29日的估算结果；
图6（r）为时间序列FVC在9月30日至10月15日的估算结果；
图6（s）为时间序列FVC在10月16日至10月31日的估算结果；
图6（t）为时间序列FVC在11月1日至11月16日的估算结果；
图6（u）为时间序列FVC在11月17日至12月2日的估算结果；
图6（v）为时间序列FVC在12月3日至12月18日的估算结果；
图6（w）为时间序列FVC在12月19日至12月31日的估算结果；
图7（a）为典型山地植被退化样带土地覆被图；
图7（b）为生长季初期HJ影像的端元选择结果；
图7（c）为生长季旺盛期HJ影像的端元选择结果；
图7（d）为生长季末期HJ影像的端元选择结果；
图7（e）为图7（b）局部放大示意；
图7（f）为图7（c）局部放大示意；
图7（g）为图7（d）局部放大示意；
图8（a）为无人机植被覆盖度参考数据与传统LSMM模型的密度散点图；
图8（b）为无人机植被覆盖度与本模型之间的密度散点图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
实施例1
经申请人研究发现，在众多的光谱分解模型中，线性光谱混合模型已经得到了广泛的应用；据研究结果表明，在一定程度上，线性光谱混合模型优于其他方法并能够获得精度较高的估算结果，传统线性混合模型中，对整体影像而言一般采用一套固定的端元进行像元光谱分解。
然而，在现实中，许多像元的光谱组合往往仅是端元库中的少数端元组成，尤其在复杂山地，地表更加破碎，如图1所示，图1展示了典型复杂山地无人机拍摄获取的高分辨率影像及对应的30米分辨率卫星影像像元格网，可以看出，30米分辨率卫星影像格网中包含了多种地物特征。同时，遥感像元端元的组成同样会随着植被的季相和地表条件变化而发生端元组合的变化。因此，理论上，采用真实的端元组合进行像元分解会相应的提高其分解精度。
此外，在线性混合模型中，端元的数量往往还会受到影像波段维度的限制。因此，在选择端元数量时还需要考虑到端元数量和模型最优精度之间的平衡。较多的端元能够更多的解释地表的异质性并因此提高模型的可应用性。然而，过多的端元也有可能导致模型对于端元的选择过于敏感而降低了其通用性。
因此，本发明实施例提供一种端元自适应的山地植被覆盖度遥感反演方法，参考图3所示，步骤具体如下：
首先，需要说明的是，对植被覆盖度反演，定义端元的数量及每个端元的光谱特征是求解线性混合像元分解模型的前提。本发明实施例考虑到复杂山地土地覆被特征和影像波段设置，选择植被、裸土、不透水面和水体四类端元作为光谱分解的基本端元。
进一步地，本发明实施例利用地物光谱仪野外测量获得的不同端元光谱数据进行端元光谱反射率的确认；野外测量获得的不同反射率曲线采用多次采样数据的平均值作为端元的光谱曲线。需要说明的是，由于卫星传感器有其各自的光谱响应函数和波段设置，因此测量获得的光谱数据需要根据传感器的波段响应函数进行地表反射率的转换，表达式如下：
Figure 242327DEST_PATH_IMAGE001
上式中，
Figure 672171DEST_PATH_IMAGE002
表示传感器在i波段的地表反射率，
Figure 179376DEST_PATH_IMAGE003
表示传感器在i波段的光谱响应函数，
Figure 985658DEST_PATH_IMAGE004
表示地表在波长
Figure 629129DEST_PATH_IMAGE005
处的反射率值，
Figure 495453DEST_PATH_IMAGE006
表示波长的积分，
Figure 224375DEST_PATH_IMAGE007
表示传感器在i波段的光谱上界，
Figure 834348DEST_PATH_IMAGE008
表示传感器在i波段的光谱下界。
进一步地，通过野外光谱测量获得全部端元光谱数据构建端元光谱库，基于端元光谱库中所有的端元光谱数据，对原始像元中的待分解像元执行逐像元基于最小光谱欧氏距离的自适应选择，得到第一端元；具体包括：为各待分解像元与端元光谱库中每一个端元计算光谱欧氏距离，表达式如下：
Figure 597905DEST_PATH_IMAGE009
上式中，
Figure 494185DEST_PATH_IMAGE010
表示待分解像元t与端元e的光谱欧氏距离，m表示光谱波段数，
Figure 975982DEST_PATH_IMAGE011
表示待分解像元t在波段b上的光谱反射率，
Figure 124067DEST_PATH_IMAGE012
表示端元e在波段b上的光谱反射率；
对比待分解像元与每个端元的光谱距离，选择最小光谱欧氏距离的端元，作为第一端元。
进一步地，进入迭代过程，定义拟合效果指数评价已选端元的光谱信息与待分解像元光谱信息的接近程度；从原始像元的光谱反射率中扣除已选择端元的影响，进行下一端元的自适应选择，确定待分解像元的所有端元组成；通过预定义的拟合效果指数判定是否结束迭代，选择指数最小值对应的端元作为自适应选择端元，预定义的拟合效果指数的表达式如下：
Figure 747989DEST_PATH_IMAGE013
上式中，
Figure 690538DEST_PATH_IMAGE014
表示第k个端元的拟合效果指数，m表示光谱波段数，
Figure 394051DEST_PATH_IMAGE015
表示待分解像元t在波段i上的光谱反射率，
Figure 611406DEST_PATH_IMAGE016
表示端元e在波段i上的光谱反射率，
Figure 83976DEST_PATH_IMAGE017
表示端元e光谱反射率归一化后的贡献权重，若满足如下表达式，则结束迭代，表达式如下：
Figure 463004DEST_PATH_IMAGE019
需要说明的是，随着端元的不断增加，拟合效果指数会逐渐降低并当像元最接近原始像元反射率值时达到最小；随着错误端元的不断引入，拟合效果指数会再次增加。若
Figure 653814DEST_PATH_IMAGE034
未能达到其最小值，则本算法将持续迭代，从未被选择的端元库中再次选择。通过上述步骤，可确定像元最终的端元构成。
上述贡献权重的计算表达式如下：
Figure 143702DEST_PATH_IMAGE018
考虑到端元的光谱矢量是非相交的，因此在本实施例中定义一个可调系数
Figure 736357DEST_PATH_IMAGE035
以增加算法的灵活度，该贡献权重的计算表达式如下：
Figure 286287DEST_PATH_IMAGE036
上式中，可调系数
Figure 698814DEST_PATH_IMAGE035
的值介于0和1之间，本实施例设置为0.5。需要注意的是，算法选择的第一端元的反射率会被首先用于计算拟合效果指数，并将其记为
Figure 992392DEST_PATH_IMAGE037
。
进一步地，根据线性混合像元分解模型假设，假设任一给定像元的地表反射率是其基本组成成分，即端元光谱以其面积比例为权重的线性组合，采用线性混合像元分解模型对所述自适应选择端元进行分解，得到植被覆盖度的反演结果。
所述线性混合像元分解模型，表达式如下：
Figure 439554DEST_PATH_IMAGE020
上式中，
Figure 425964DEST_PATH_IMAGE021
表示像元I在波段
Figure 591366DEST_PATH_IMAGE023
的光谱反射率，m表示端元的数量，
Figure 423056DEST_PATH_IMAGE024
表示像元I的第k个端元的面积比例，
Figure 724724DEST_PATH_IMAGE025
表示第k个端元在波段
Figure 616457DEST_PATH_IMAGE023
的光谱反射率，
Figure 3576DEST_PATH_IMAGE026
表示线性混合像元分解模型的拟合残差，其中，端元的面积比例即为该像元的植被覆盖度。
为了使模型的物理意义更加明显，本发明实施例为所述线性混合像元分解模型设定约束条件，使得所有的混合分解结果之和为1，以及面积比例的非负性，表达式如下：
Figure 638957DEST_PATH_IMAGE038
为了评价模型模拟结果的精度，本发明实施例对反演结果进行验证，验证公式采用均方根误差RMSE、平均误差AD以及可决系数R2来评价模型的全局误差、不同模型之前的偏差、反演数据与参考数据之间的相关关系。具体表达式如下：
Figure 326290DEST_PATH_IMAGE039
Figure 388924DEST_PATH_IMAGE040
Figure 263339DEST_PATH_IMAGE041
上式中，N表示待分解像元数，
Figure 436831DEST_PATH_IMAGE042
表示像元I的模型估算值，
Figure 713092DEST_PATH_IMAGE043
表示像元I的地面实测参考值或无人机提取参考值，
Figure 946627DEST_PATH_IMAGE033
表示像元I的参考值均值。
以下结合附图以及国产环境减灾卫星影像植被覆盖度反演为具体实施例对对本发明作进一步地具体描述，本发明所述方法包括以下步骤：
野外采用SVC-1024地物光谱仪实测的不同端元地物反射率曲线，野外测量获得的山区植被、裸土、不透水面和水体的反射率曲线如图4（a）至4（d）所示，图中示出了HJ-1A/B卫星CCD传感器对应波长400到1000nm的光谱数据。可以明显看出，不同端元的反射率曲线变化明显，通过地表反射率进行转换，转换后的对应波段地表反射率数据如图4中菱形点所示。误差棒代表多次测量结果的标准差。
通过依次计算待分解像元与各端元之间的欧式光谱距离、拟合效果指数、端元贡献权重、剩余光谱信号，通过算法迭代，计算多次迭代后的拟合效果指数并进行对比，确定最终待分解像元的端元组成。图5（a）至图5（d）示出了山区湿地典型退化样带三种代表植被不同生长阶段和湿地淹水阶段不同日期（生长季初期、生长季旺盛期和生长季末期）HJ-1A/B卫星合成影响的端元选择结果。图中同样示出了退化样带30米空间分辨率的土地覆盖图作为对比。图5（b）至5（d）中左上角为端元选择结果及其对应的放大效果。与土地覆盖图相比，不同土地覆被类型的端元信息随着时相有着显著的变化。
进一步地，基于线性光谱分解模型和端元自适应选择结果，将提出的自适应端元选择模型应用到16天融合的30米分辨率HJ-1A/B影像中，获得对应的时间序列FVC估算结果，线性光谱分解模型原理示意图参见图2
，线性光谱分解模型的原理如下：像元光谱=F（植被）+F（水体）+F（土壤），植被设置为0.5，土壤设置为0.3，水体设置为0.2，约束条件设置为0≤F≤1，∑F=1。如图6（a）至（w）所示，图6（a）至（w）示出了时间序列FVC在1月1日至12月31日的估算结果；所有图的颜色均拉伸到0至1。总体上，多时相遥感影像时空融合能够有效地增强遥感影像的时空连续性，通过S-G时相滤波，合成影像中的残云得到了进一步地剔除。本发明实施例模型得到的FVC的空间分布与地表的湿度总体一致性较强，表现出常见的FVC年度季相循环变化特征，变化趋势主要分布在草本植被区域。
进一步地评价模型估算精度，图7（a）至7（d）示出了典型山地植被退化样带的无人机获取的HJ合成影像图7（a），厘米级无人机植被覆盖度参考数据图7（b），反演获取的植被覆盖度结果7（c）以及传统线性光谱分解模型的分解结果7（d），图7（e）为图7（b）局部放大示意，图7（f）为图7（c）局部放大示意，图7（g）为图7（d）局部放大示意。从无人机植被覆盖度参考数据可以看出，植被覆盖度高值区域主要集中分布在草甸和湿草甸区域，而低值区域主要分布在草原和退化区域。总体上，模型估算结果与UAV参考影像的一致性更高。
图8（a）和8（b）示出了UAV参考FVC和不同模型间拟合效果的密度散点图。总体上，本方法模型估算得出的FVC和UAV参考FVC更接近1:1线，并显示出较小的估算偏差。
基于CPU：Intel(R) Core(TM) i9-10920X CPU @ 3.50GHz 3.50 GHz, RAM：16GB，软件平台Windows 10 64 bit的硬件平台，本方法实际测试结果为，单景影像的处理时间为7.8分钟，全年以16天为时间窗口的23期无云影像3小时内即可反演完成。本方案能够大大提高工作效率，可快速完成多幅卫星影像的处理工作，实现复杂山人力无法到达区域的植被覆盖度高效、高精度反演，减少人力、物力，在海量环境减灾卫星的数据处理中有十分重要的意义。
本发明所提供的山地植被覆盖度遥感反演方法，考虑到复杂山地是由许多复杂的地表覆被类型所构成，通过综合高分辨率卫星影像和新的自适应端元光谱分解模型，能够实现复杂山地时间序列植被覆盖度高精度反演，并生成相应的误差报告，具有大区域、高效率、高精度、低成本的优点；本方法对海量山地多光谱遥感影像的植被覆盖度反演非常有效，解决了传统线性光谱混合分解模型端元固定、反演精度低的问题，对高时空异质性的山地区域具有很高的有效性和可靠性，对于目前在轨的卫星数据反演，对山地植被监测、生态环境健康监测等具有十分重要的意义。
以上仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
1. 一种端元自适应的山地植被覆盖度遥感反演方法，其特征在于，包括如下步骤：
S1.获取端元光谱数据并构建端元光谱库；
S2.根据端元光谱库，对原始像元中的待分解像元执行逐像元基于最小光谱欧氏距离的自适应选择，得到第一端元；
S3.进入迭代过程，定义拟合效果指数评价已选则端元的光谱信息与待分解像元光谱信息的接近程度，从原始像元的光谱反射率中扣除已选择端元的影响，进行下一端元的自适应选择；
通过预定义的拟合效果指数判定是否结束迭代，当拟合效果指数达到最小时结束迭代，选择指数最小值对应的端元作为自适应选择端元；
S4.采用线性混合像元分解模型对所述自适应选择端元进行分解，得到植被覆盖度的反演结果。
2. 如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，所述方法在步骤S1之前还包括：定义端元的数量及每个端元的光谱特征。
3.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S1在构建端元光谱库之前还包括：基于传感器的光谱响应函数，对端元光谱数据进行地表反射率的转换，表达式如下：
Figure 308338DEST_PATH_IMAGE001
上式中，
Figure 285522DEST_PATH_IMAGE002
表示传感器在i波段的地表反射率，
Figure 681868DEST_PATH_IMAGE003
表示传感器在i波段的光谱响应函数，
Figure 86305DEST_PATH_IMAGE004
表示地表在波长
Figure 935312DEST_PATH_IMAGE005
处的反射率值，
Figure 716186DEST_PATH_IMAGE006
表示波长的积分，
Figure 967039DEST_PATH_IMAGE007
表示传感器在i波段的光谱上界，
Figure 542377DEST_PATH_IMAGE008
表示传感器在i波段的光谱下界。
4.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S2具体包括：
为各待分解像元与端元光谱库中每一个端元计算光谱欧氏距离，表达式如下：
Figure 613101DEST_PATH_IMAGE009
上式中，
Figure 932087DEST_PATH_IMAGE010
表示待分解像元t与端元e的光谱欧氏距离，m表示光谱波段数，
Figure 37446DEST_PATH_IMAGE011
表示待分解像元t在波段b上的光谱反射率，
Figure 783685DEST_PATH_IMAGE012
表示端元e在波段b上的光谱反射率；
选择最小光谱欧氏距离的端元，作为第一端元。
5.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S3中所述预定义的拟合效果指数的表达式如下：
Figure 607285DEST_PATH_IMAGE013
上式中，
Figure 729961DEST_PATH_IMAGE014
表示第k个端元的拟合效果指数，m表示光谱波段数，
Figure 955406DEST_PATH_IMAGE015
表示待分解像元t在波段i上的光谱反射率，
Figure 872547DEST_PATH_IMAGE016
表示端元e在波段i上的光谱反射率，
Figure 183442DEST_PATH_IMAGE017
表示端元e光谱反射率归一化后的贡献权重。
6.如权利要求5所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，所述贡献权重的计算表达式如下：
Figure 844231DEST_PATH_IMAGE018
上式中，
Figure 924182DEST_PATH_IMAGE010
表示待分解像元t与端元e的光谱欧氏距离。
7.如权利要求6所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S3中所述通过预定义的拟合效果指数判定是否结束迭代具体包括：
若满足如下表达式，则结束迭代，表达式如下：
Figure 12224DEST_PATH_IMAGE019
8.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S4中所述线性混合像元分解模型，表达式如下：
Figure 544837DEST_PATH_IMAGE020
上式中，
Figure 9316DEST_PATH_IMAGE021
表示像元I在波段
Figure 943774DEST_PATH_IMAGE023
的光谱反射率，m表示端元的数量，
Figure 202717DEST_PATH_IMAGE024
表示像元I的第k个端元的面积比例，
Figure 222626DEST_PATH_IMAGE025
表示第k个端元在波段
Figure 490796DEST_PATH_IMAGE023
的光谱反射率，
Figure 279760DEST_PATH_IMAGE026
表示线性混合像元分解模型的拟合残差。
9.如权利要求8所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，所述方法还包括，为所述线性混合像元分解模型设定约束条件，使得所有的混合分解结果之和为1，以及面积比例的非负性，表达式如下：
Figure 975184DEST_PATH_IMAGE038
10.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，所述方法还包括，对反演结果进行验证，验证公式采用均方根误差RMSE、平均误差AD以及可决系数R2，具体表达式如下：
Figure 482389DEST_PATH_IMAGE039
Figure 288671DEST_PATH_IMAGE040
Figure 932141DEST_PATH_IMAGE041
上式中，N表示待分解像元数，
Figure 816044DEST_PATH_IMAGE042
表示像元I的模型估算值，
Figure 544966DEST_PATH_IMAGE043
表示像元I的地面实测参考值或无人机提取参考值，
Figure 154939DEST_PATH_IMAGE033
表示像元I的参考值均值。

Claims (10)

1.一种端元自适应的山地植被覆盖度遥感反演方法，其特征在于，包括如下步骤：
S1.获取端元光谱数据并构建端元光谱库；
S2.根据端元光谱库，对原始像元中的待分解像元执行逐像元基于最小光谱欧氏距离的自适应选择，得到第一端元；
S3.进入迭代过程，定义拟合效果指数评价已选则端元的光谱信息与待分解像元光谱信息的接近程度，从原始像元的光谱反射率中扣除已选择端元的影响，进行下一端元的自适应选择；
通过预定义的拟合效果指数判定是否结束迭代，当拟合效果指数达到最小时结束迭代，选择指数最小值对应的端元作为自适应选择端元；
S4.采用线性混合像元分解模型对所述自适应选择端元进行分解，得到植被覆盖度的反演结果。
2.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，所述方法在步骤S1之前还包括：定义端元的数量及每个端元的光谱特征。
3.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S1在构建端元光谱库之前还包括：基于传感器的光谱响应函数，对端元光谱数据进行地表反射率的转换，表达式如下：
Figure FDA0003797564260000011
上式中，ρi表示传感器在i波段的地表反射率，fi(λ)表示传感器在i波段的光谱响应函数，ρ(λ)表示地表在波长λ处的反射率值，d(λ)表示波长的积分，λmax表示传感器在i波段的光谱上界，λmin表示传感器在i波段的光谱下界。
4.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S2具体包括：
为各待分解像元与端元光谱库中每一个端元计算光谱欧氏距离，表达式如下：
Figure FDA0003797564260000021
上式中，dt，e表示待分解像元t与端元e的光谱欧氏距离，m表示光谱波段数，ρt，b表示待分解像元t在波段b上的光谱反射率，ρe，b表示端元e在波段b上的光谱反射率；
选择最小光谱欧氏距离的端元，作为第一端元。
5.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S3中所述预定义的拟合效果指数的表达式如下：
Figure FDA0003797564260000022
上式中，Fk表示第k个端元的拟合效果指数，m表示光谱波段数，Ri，t表示待分解像元t在波段i上的光谱反射率，Ri，e表示端元e在波段i上的光谱反射率，Wt，e表示端元e光谱反射率归一化后的贡献权重。
6.如权利要求5所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，所述贡献权重的计算表达式如下：
Figure FDA0003797564260000031
上式中，dt，e表示待分解像元t与端元e的光谱欧氏距离。
7.如权利要求6所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S3中所述通过预定义的拟合效果指数判定是否结束迭代具体包括：
若满足如下表达式，则结束迭代，表达式如下：
Fk-1≥Fk≥Fk+1。
8.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，步骤S4中所述线性混合像元分解模型，表达式如下：
Figure FDA0003797564260000032
上式中，RI，λ表示像元I在波段λ的光谱反射率，m表示端元的数量，fI，k表示像元I的第k个端元的面积比例，ρk，λ表示第k个端元在波段λ的光谱反射率，εI，λ表示线性混合像元分解模型的拟合残差。
9.如权利要求8所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，所述方法还包括，为所述线性混合像元分解模型设定约束条件，使得所有的混合分解结果之和为1，以及面积比例的非负性，表达式如下：
Figure FDA0003797564260000041
10.如权利要求1所述的端元自适应的山地植被覆盖度遥感反演方法，其特征在于，所述方法还包括，对反演结果进行验证，验证公式采用均方根误差RMSE、平均误差AD以及可决系数R2，具体表达式如下：
Figure FDA0003797564260000042
Figure FDA0003797564260000043
Figure FDA0003797564260000044
上式中，N表示待分解像元数，
Figure FDA0003797564260000045
表示像元I的模型估算值，
Figure FDA0003797564260000046
表示像元I的地面实测参考值或无人机提取参考值，
Figure FDA0003797564260000047
表示像元I的参考值均值。
CN202210686048.7A 2022-06-17 2022-06-17 一种端元自适应的山地植被覆盖度遥感反演方法 Active CN114780904B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210686048.7A CN114780904B (zh) 2022-06-17 2022-06-17 一种端元自适应的山地植被覆盖度遥感反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210686048.7A CN114780904B (zh) 2022-06-17 2022-06-17 一种端元自适应的山地植被覆盖度遥感反演方法

Publications (2)

Publication Number Publication Date
CN114780904A CN114780904A (zh) 2022-07-22
CN114780904B true CN114780904B (zh) 2022-09-27

Family

ID=82420481

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210686048.7A Active CN114780904B (zh) 2022-06-17 2022-06-17 一种端元自适应的山地植被覆盖度遥感反演方法

Country Status (1)

Country Link
CN (1) CN114780904B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117218551B (zh) * 2023-08-04 2024-04-09 华南师范大学 一种基于误差分析的估测算法优化方法及装置
CN116778345B (zh) * 2023-08-18 2023-11-14 中国科学院、水利部成都山地灾害与环境研究所 山地森林火烧迹地植被盖度反演方法、装置及电子设备

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105224915A (zh) * 2015-09-07 2016-01-06 河海大学 一种高光谱图像混合像元分解方法
CN109035199A (zh) * 2018-06-21 2018-12-18 中国科学院西安光学精密机械研究所 基于空间特征的高光谱数据端元提取方法、计算机可读存储介质、电子设备
CN110288647A (zh) * 2019-06-25 2019-09-27 中国水利水电科学研究院 一种基于高分辨率卫星数据监测灌区灌溉面积方法
CN113310904A (zh) * 2021-06-15 2021-08-27 东南大学 一种植被覆盖下的土壤光谱还原的图像处理方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6608931B2 (en) * 2001-07-11 2003-08-19 Science Applications International Corporation Method for selecting representative endmember components from spectral data
CN102269576A (zh) * 2010-06-03 2011-12-07 曹春香 一种森林覆盖度及有效叶面积指数的主被动协同反演方法
CN103544477B (zh) * 2013-09-30 2016-11-02 北京师范大学 基于改进的线性光谱混合模型的植被覆盖度估算方法
CN104567754A (zh) * 2014-12-26 2015-04-29 河南省农业科学院农业经济与信息研究所 一种耦合星-地遥感的小麦叶面积指数估算方法
CN106124454B (zh) * 2016-06-30 2018-10-16 国交空间信息技术（北京）有限公司 一种基于遥感影像的沥青路面老化状况监测方法
CN107590800A (zh) * 2017-09-04 2018-01-16 东华理工大学 一种基于地物光谱库的高光谱遥感图像混合像元分解方法
CN112395989B (zh) * 2020-11-18 2023-07-14 北京师范大学 一种用于多卫星传感器的积雪覆盖度混合像元分解方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105224915A (zh) * 2015-09-07 2016-01-06 河海大学 一种高光谱图像混合像元分解方法
CN109035199A (zh) * 2018-06-21 2018-12-18 中国科学院西安光学精密机械研究所 基于空间特征的高光谱数据端元提取方法、计算机可读存储介质、电子设备
CN110288647A (zh) * 2019-06-25 2019-09-27 中国水利水电科学研究院 一种基于高分辨率卫星数据监测灌区灌溉面积方法
CN113310904A (zh) * 2021-06-15 2021-08-27 东南大学 一种植被覆盖下的土壤光谱还原的图像处理方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于面向对象多端元混解模型的植被覆盖度反演及其时空分布研究;李哲 等;《遥感技术与应用》;20181231;第33卷(第6期);第1149-1158页 *
端元可变非线性混合像元分解模型;李慧 等;《测绘学报》;20160131;第45卷(第1期);第80-86页 *

Also Published As

Publication number Publication date
CN114780904A (zh) 2022-07-22

Similar Documents

Publication Publication Date Title
CN114780904B (zh) 一种端元自适应的山地植被覆盖度遥感反演方法
CN107909607B (zh) 一种年度区域植被覆盖度计算方法
Weng et al. Generating daily land surface temperature at Landsat resolution by fusing Landsat and MODIS data
Mucher et al. Land cover characterization and change detection for environmental monitoring of pan-Europe
CN110751727B (zh) 一种基于Landsat长时间序列的合成图像构建方法
Kuusinen et al. Effects of forest age on albedo in boreal forests estimated from MODIS and Landsat albedo retrievals
CN105809148B (zh) 基于遥感时空谱融合的作物干旱识别及风险评估方法
CN110927120B (zh) 一种植被覆盖度预警方法
CN111337434A (zh) 一种矿区复垦植被生物量估算方法及系统
CN103544477A (zh) 基于改进的线性光谱混合模型的植被覆盖度估算方法
CN112348812A (zh) 林分年龄信息测量方法及装置
CN115077656B (zh) 水库水储量反演方法和装置
CN113466143B (zh) 土壤养分反演方法、装置、设备及介质
CN114821349A (zh) 顾及谐波模型系数和物候参数的森林生物量估算方法
CN112380980A (zh) 一种人工龙竹林lai遥感估测最佳尺度选择方法
Zhao et al. Retrieval of red solar-induced chlorophyll fluorescence with TROPOMI on the Sentinel-5 precursor mission
Cilek et al. The use of regression tree method for Sentinel-2 satellite data to mapping percent tree cover in different forest types
CN114062439A (zh) 一种利用时间序列遥感影像联合估算土壤剖面盐分的方法
Lovell et al. Analysis of POLDER-ADEOS data for the Australian continent: The relationship between BRDF and vegetation structure
Gross et al. Application of spatial resolution enhancement and spectral mixture analysis to hyperspectral images
Danoedoro et al. Combining pan-sharpening and forest cover density transformation methods for vegetation mapping using Landsat-8 Satellite Imagery
Mngadi et al. Quantifying carbon stock variability of species within a reforested urban landscape using texture measures derived from remotely sensed imagery
CN112052720B (zh) 一种基于直方图聚类的高时空归一化植被指数ndvi的融合模型
Scarth et al. Integrating high and moderate spatial resolution image data to estimate forest age structure
Nasirzadehdizaji et al. Crop mapping improvement by combination of optical and SAR datasets

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
GR01 Patent grant
GR01 Patent grant