CN115203624A - 一种基于时序遥感的任意时刻地表环境综合评估方法 - Google Patents

一种基于时序遥感的任意时刻地表环境综合评估方法 Download PDF

Info

Publication number
CN115203624A
CN115203624A CN202210810131.0A CN202210810131A CN115203624A CN 115203624 A CN115203624 A CN 115203624A CN 202210810131 A CN202210810131 A CN 202210810131A CN 115203624 A CN115203624 A CN 115203624A
Authority
CN
China
Prior art keywords
remote sensing
pixel
index
time sequence
wave band
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.)
Granted
Application number
CN202210810131.0A
Other languages
English (en)
Other versions
CN115203624B (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.)
Ningbo University
Original Assignee
Ningbo University
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 Ningbo University filed Critical Ningbo University
Priority to CN202210810131.0A priority Critical patent/CN115203624B/zh
Publication of CN115203624A publication Critical patent/CN115203624A/zh
Application granted granted Critical
Publication of CN115203624B publication Critical patent/CN115203624B/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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N2021/1793—Remote sensing
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Health & Medical Sciences (AREA)
  • Mathematical Optimization (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于时序遥感的任意时刻地表环境综合评估方法，具体步骤为：（1）、对研究时期内选定区域的实拍遥感影像进行预处理，并将预处理后的全部影像排列构成影像时间序列；（2）、对影像中的每一像元分别建立初始时间序列模型；（3）、对像元的初始时间序列模型进行变化检测；（4）、获得任意时刻的预测遥感影像；（5）、将五个遥感指数作为选定区域的地表环境综合表征；（6）、将遥感指数去噪和归一化；（7）、获得地表环境综合指数；优点是本方法使地表环境定期评估工作成为可能，有效规避了由于影像缺失导致定期评估难以开展的问题，同时解决了由于影像成像时间差异导致的评估结果难以对比的难题。

Description

一种基于时序遥感的任意时刻地表环境综合评估方法
技术领域
本发明涉及地表环境的遥感评估方法，尤其涉及一种基于时序遥感的任意时刻地表环境综合评估方法。
背景技术
地表环境与人类生活品质以及当今社会可持续发展高度息息相关。半个世纪以来，我国快速城市化进程导致了一系列环境问题，如森林砍伐、多样性丧失、空气和水污染、土地荒漠化以及城市热岛效应等。鉴于城市化是人类文明发展过程中不可逆转的进程，城市化对生态环境的胁迫作用难以在短期内缓解。因此，及时、客观地探索城市化过程中地表环境的变化情况，对于制定与区域环境容量相协调的社会经济发展规划尤为重要，这成为近年来城市生态学研究的核心议题。
遥感技术具有覆盖范围广泛、观测尺度多样、光谱信息丰富、数据获取便利等优势，正日益成为国内外地表环境表征的重要手段。早期对地表环境描述主要侧重于基于遥感的单一指数的计算，主要强调地表环境的某个特定方面。典型的例子是归一化差异植被指数(NDVI)和地表温度(LST)经常被用来描述植被覆盖和城市热岛状况。而在大多数情况下，地表环境变化十分复杂，很难用单一的遥感指数来表示。因此，部分学者提出了综合多种基于遥感的指数的模型方案，例如森林干扰指数(FDI)、全球干扰指数(MGDI)和生态位建模(ENM)等。尽管这些综合指数能够更全面地反映地表环境状况，但仍然存在一系列挑战，如权重选取。实际上，确定各个指标对地表环境变化地作用大小并依比例分配权重是十分困难的，权重分配难免因主观经验差异而不同，不同权重又将直接影响到地表环境评价结果。此外，由于涉及一系列指数(如绿度、热度)，上述综合指数模型方案对年内物候非常敏感，多数研究采取在年度影像中挑选质量较好且季节相近的数据，或进行小范围的研究，但是对于云量较多的沿海城市区域往往存在可用数据缺失，严重限制了上述方法应用的准确性、可比性以及研究的区域范围和时序长度。
发明内容
本发明所要解决的技术问题是提供一种基于时序遥感的任意时刻地表环境综合评估方法，其自动化程度高、鲁棒性强，使得任意时刻地表环境评估成为可能，有效规避了由于影像缺失导致定期评估难以开展的问题，同时解决了影像成像时间差异导致的评估结果难以对比的难题。
本发明解决上述技术问题所采用的技术方案为：一种基于时序遥感的任意时刻地表环境综合评估方法，包括以下具体步骤：
(1)、对研究时期内选定区域的实拍遥感影像进行预处理，并将预处理后的全部影像按照成像时间顺序依次排列构成影像时间序列；
(2)、对影像中的每一像元根据时间序列在所有波段上分别建立初始时间序列模型；
(3)、对像元的初始时间序列模型进行变化检测，若变化没有发生，则保持初始时间序列模型；若变化发生，则将变化发生时间作为时间序列模型的断点，对像元的初始时间序列模型进行分段，得到分段时间序列模型，以反映对应像元的地表反射率变化；
(4)、在像元的时间序列模型的基础上，获得任意时刻的预测遥感影像；
(5)、基于预测遥感影像，将植被覆盖度VC、植被健康度VHI、土壤湿度Wet、地表土建度NDBSI和地表温度LST五个遥感指数作为选定区域的地表环境综合表征，并计算五个遥感指数的数值；
(6)、将遥感指数去噪和归一化；
(7)、获得地表环境综合指数。
进一步地，所述的步骤(1)的实拍遥感影像预处理的具体方法为：下载研究时期内选定区域中全部经过几何校正和辐射校正的实拍光学遥感影像数据，或下载全部原始实拍光学遥感影像数据后人工进行几何校正和辐射校正，然后利用遥感影像的质量评估波段，剔除每景影像中被云、云阴影或雪覆盖的无效观测值。
进一步地，所述的步骤(2)的具体建立方法为：每一像元的初始时间序列模型由线性函数和谐波函数组成，依据每一像元在时间序列中的有效观测值数量，从粗到细构建三种关系式，若12≤有效观测值数量＜18，则模型选用关系式(1)；若18≤有效观测值数量＜24，则模型选用关系式(2)；若有效观测值数量≥24，则模型选用关系式(3)；
Figure BDA0003740459250000021
Figure BDA0003740459250000031
Figure BDA0003740459250000032
其中，
Figure BDA0003740459250000033
表示像元的预测地表反射率；n表示初始时间序列模型所选用的关系式的标号，n＝1,2,3；i为遥感影像第i个波段；t为儒略日；w为频率，
Figure BDA0003740459250000034
c0,i和c1,i为截距系数和斜率系数，an,i和bn,i是第n次谐波系数，且c0,i、c1,i、an,i和bn,i这四个参数均通过最小二乘法确定最优取值。
进一步地，所述的步骤(3)的具体检测方法为：为了准确检测像元的初始时间序列模型的变化，变化的发生时间通过比较像元的初始时间序列模型预测值和后续有效观测值的残差resid与初始时间序列模型的均方根误差RMSE确定，即
Figure BDA0003740459250000035
其中：D为遥感影像的波段集合，i为遥感影像的第i个波段，residi为第i个波段的初始时间序列模型预测值与后续有效观测值的残差，RMSEi为第i个波段的像元的初始时间序列模型的均方根误差，
Figure BDA0003740459250000036
为阈值，等于自由度为影像波段数的逆卡方分布数；
将同一波段中同一像元的6个连续的后续有效观测值作为一组进行判别，若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值均大于阈值，则将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值出现的时间作为变化发生时间，并将该变化发生时间作为该像元在每个波段中的时间序列模型中的断点，对该像元在每个波段中的初始时间序列模型进行分段，并根据步骤(2)中的关系式得到分段时间序列模型，以反映该像元的地表反射率变化；
若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值没有全部大于阈值，则认为变化没有发生，并将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值加入对应波段的初始时间序列模型中，然后根据步骤(2)中的关系式将像元的初始时间序列模型拟合更新；
上述判别过程重复进行，直到全部有效观测值都位于像元的初始时间序列模型或分段时间序列模型中为止。
进一步地，所述的步骤(4)的具体方法为：逐波段逐像元定位像元的初始时间序列模型或分段时间序列模型，根据需要确定地表环境综合评估的儒略日t，再通过步骤(2)的关系式得到像元的预测地表反射率，然后将全部像元组成波段，将全部波段组成影像，获得任意时刻的预测遥感影像。
进一步地，所述的步骤(5)中的五个遥感指数：植被覆盖度VC、植被健康度VHI、土壤湿度Wet、地表土建度NDBSI和地表温度LST的关系式为：
Figure BDA0003740459250000041
Figure BDA0003740459250000042
其中，ρNIR和ρRed分别为预测遥感影像的近红外波段和红光波段，NDVI为归一化差异植被指数，NDVIveg和NDVIsoil分别为植被和土壤的归一化差异植被指数，两者分别为研究范围内NDVI的最大和最小值；
Figure BDA0003740459250000043
Figure BDA0003740459250000044
VHI＝PCA1(NDVI,NRI,NDSVI) (9)
其中，ρSWIR1和ρGreen分别为预测遥感影像的短波红外第1波段和绿光波段，NRI为氮元素的反射率指数，NDSVI为归一化差异植被衰亡指数，植被健康度VHI为通过将NDVI、NRI、NDSVI主成分变换后的第一分量PCA1；
Figure BDA0003740459250000045
其中，D是遥感影像波段集合，i为预测遥感影像第i个波段，ai为波段线性组合系数，ρi为第i个波段的影像地表反射率；
Figure BDA0003740459250000046
Figure BDA0003740459250000047
Figure BDA0003740459250000048
其中，ρBlue为预测遥感影像的蓝光波段，IBI为建成区指数，SI为土壤指数；
Figure BDA0003740459250000051
Figure BDA0003740459250000052
其中，h为普朗克常数，h＝6.26×10-34J-sec；c为光速，c＝2.998×108m/sec；K为斯忒潘-玻尔兹曼常数，K＝1.38×10-23J/K；α为中间变量，λ为热红外波段中心波长，ε为地表出射率；Tb为亮度温度。
进一步地，所述的步骤(6)的具体方法为：
对遥感指数VC、VHI、Wet、NDBSI、LST采用百分位去噪方法，设定百分位为a，将全部研究时期所包含的像元数值按从小到大排序，取a％作为有效数值下限，(100-a)％为有效数值上限，若像元数值小于有效数值下限，则取有效数值下限；若像元数值大于有效数值上限，则取有效数值上限；然后通过关系式(16)对五个遥感指数进行归一化，以屏蔽不同遥感指数数值范围差异的影响；
Figure BDA0003740459250000053
其中，Index为原始遥感指数，Indexlb和Indexub分别为遥感指数的下限和上限，Indexnom为归一化后的遥感指数。
进一步地，所述的步骤(7)的具体方法为：
将去噪归一化后的各遥感指数通过最小噪声分离变换，并取其第一分量作为初步的地表环境综合指数LSECI0，观察第一分量中各遥感指数的载荷数值，若表示植被覆盖度VC的载荷数值大于等于零，则LSECI0为最终的地表环境综合指数LSECI；若表示植被覆盖度VC的载荷数值为负，则将LSECI0取反获得最终的地表环境综合指数LSECI,其关系式为：
LSECI0＝MNF1(VC,VHI,Wet,NDBSI,LST) (17)
Figure BDA0003740459250000054
其中，MNF1为最小噪声分离变换的第一分量，q(VC)为在第一分量中植被覆盖度VC的载荷数值。
与现有技术相比，本发明的优点是：
(1)、本方法中通过对像元的时间序列模型构建和变化检测，能够生成任意时刻的预测遥感影像，攻克了常规光学影像受云雨覆盖而可用性不佳的瓶颈，使得云雨天气频发的沿海区域对地表长时间大范围连续观测成为可能。
(2)、提出了地表环境综合指数LSECI，采用最小噪声分离变换自动客观综合了植被覆盖度、植被健康度、土壤湿度、地表土建度和地表温度五个指标，解决了指标人为赋权的评价结果不确定性，能够在空间上量化评估结果并使评估结果具备跨区域可比性。
(3)、本方法使地表环境定期评估工作成为可能，有效规避了由于影像缺失导致定期评估难以开展的问题，同时解决了由于影像成像时间差异导致的评估结果难以对比的难题。而且该方法自动化程度高、鲁棒性强，适用结合多种平台、不同空间和时间分辨率的遥感数据，有望将地表环境综合评估扩展至精细的时间尺度。
附图说明
图1为本发明实例一的像元时序模型构建与变化检测过程，其中：
(a1～a3)为利用分段线性谐波函数模型捕捉土地利用/覆被变化(LUCC)，以SWIR1波段的反射率为例，绿色和蓝色曲线对应LUCC前后不同的时序模型分段；
(b1～b3)为从(a1～a3)中的红色矩形标记的时期获取的谷歌地球历史快照，黄色大头针对应于时序模型和变化检测的像元；
(c1～c2)为2010年7月1日和2018年7月1日由模型系数生成的预测Landsat卫星影像；
(d1～d2)为Globeland30土地利用分类图，红色十字对应于时序模型和变化检测的像元。
图2为通过不同方法归一化NDVI的比较，其中：
(a1～c1)为1987年和2017年的原始NDVI及其直方图分布；
(a2～c2)为1987年和2017年的归一化NDVI，使用单个年份的最大值和最小值(a3中的虚线标记)及其直方图分布计算；
(a3～c3)为1987年至2017年的归一化NDVI，在百分比去噪(在a3中用虚线标记)归一化后生成，以及它们的直方图分布。
图3为杭州市1986～2019年环境综合质量指数LSECI的空间分异，展现了对应年份7月1日的每年LSECI。
图4为杭州市2019年半月综合质量指数LSECI的空间分异，展现了相应月份的第10天和第20天，春、夏、秋、冬四季分别对应于第一至第四列。
具体实施方式
以下结合附图实施例对本发明作进一步详细描述。
如图所示，一种基于时序遥感的任意时刻地表环境综合评估方法，包括以下具体步骤：
(1)、下载研究时期内选定区域中全部经过几何校正和辐射校正的实拍光学遥感影像数据，或下载全部原始实拍光学遥感影像数据(如：Landsat、Sentinel、MODIS等)后人工进行几何校正和辐射校正，然后利用遥感影像的质量评估波段，剔除每景影像中被云、云阴影或雪覆盖的无效观测值，再将预处理后的全部影像按照成像时间顺序依次排列构成影像时间序列；
(2)、对影像中的每一像元根据时间序列在所有波段上分别建立初始时间序列模型，具体为：每一像元的初始时间序列模型由线性函数和谐波函数组成，依据每一像元在时间序列中的有效观测值数量，从粗到细构建三种关系式，若12≤有效观测值数量＜18，则模型选用关系式(1)；若18≤有效观测值数量＜24，则模型选用关系式(2)；若有效观测值数量≥24，则模型选用关系式(3)；
Figure BDA0003740459250000071
Figure BDA0003740459250000072
Figure BDA0003740459250000073
其中，
Figure BDA0003740459250000074
表示像元的预测地表反射率；n表示初始时间序列模型所选用的关系式的标号，n＝1,2,3；i为遥感影像第i个波段；t为儒略日；w为频率，
Figure BDA0003740459250000075
c0,i和c1,i为截距系数和斜率系数，代表逐渐的年际变化；an,i和bn,i是第n次谐波系数，代表周期性的年内变化，且c0,i、c1,i、an,i和bn,i这四个参数均通过最小二乘法确定最优取值；
(3)、对像元的初始时间序列模型进行变化检测：
为了准确检测像元的初始时间序列模型的变化，变化的发生时间通过比较像元的初始时间序列模型预测值和后续有效观测值的残差resid与初始时间序列模型的均方根误差RMSE确定，即
Figure BDA0003740459250000081
其中：D为遥感影像的波段集合，i为遥感影像的第i个波段，residi为第i个波段的初始时间序列模型预测值与后续有效观测值的残差，RMSEi为第i个波段的像元的初始时间序列模型的均方根误差，
Figure BDA0003740459250000082
为阈值，等于自由度为影像波段数的逆卡方分布数；
将同一波段中同一像元的6个连续的后续有效观测值作为一组进行判别，若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值均大于阈值，则将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值出现的时间作为变化发生时间，并将该变化发生时间作为该像元在每个波段中的时间序列模型中的断点，对该像元在每个波段中的初始时间序列模型进行分段，并根据步骤(2)中的关系式得到分段时间序列模型，以反映该像元的地表反射率变化；
若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值没有全部大于阈值，则认为变化没有发生，并将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值加入对应波段的初始时间序列模型中，然后根据步骤(2)中的关系式将像元的初始时间序列模型拟合更新；
上述判别过程重复进行，直到全部有效观测值都位于像元的初始时间序列模型或分段时间序列模型中为止；
(4)、在像元的时间序列模型的基础上，获得任意时刻的预测遥感影像，具体为：
逐波段逐像元定位像元的初始时间序列模型或分段时间序列模型，根据研究的需要确定地表环境综合评估的儒略日t，t可为研究时期内的任意时间，再通过步骤(2)的关系式得到像元的预测地表反射率，然后将全部像元组成波段，将全部波段组成影像，获得任意时刻的预测遥感影像；
(5)、基于预测遥感影像，将植被覆盖度VC、植被健康度VHI、土壤湿度Wet、地表土建度NDBSI和地表温度LST五个遥感指数作为选定区域的地表环境综合表征，并计算五个遥感指数的数值，具体关系式为：
Figure BDA0003740459250000091
Figure BDA0003740459250000092
其中，ρNIR和ρRed分别为预测遥感影像的近红外波段和红光波段，NDVI为归一化差异植被指数，NDVIveg和NDVIsoil分别为植被和土壤的归一化差异植被指数，两者分别为研究范围内NDVI的最大和最小值；
Figure BDA0003740459250000093
Figure BDA0003740459250000094
VHI＝PCA1(NDVI,NRI,NDSVI) (9)
其中，ρSWIR1和ρGreen分别为预测遥感影像的短波红外第1波段和绿光波段，NRI为氮元素的反射率指数，NDSVI为归一化差异植被衰亡指数，植被健康度VHI为通过将NDVI、NRI、NDSVI主成分变换后的第一分量PCA1；
Figure BDA0003740459250000095
其中，D是遥感影像波段集合，i为预测遥感影像第i个波段，ai为波段线性组合系数，ρi为第i个波段的影像地表反射率；
Figure BDA0003740459250000096
Figure BDA0003740459250000097
Figure BDA0003740459250000098
其中，ρBlue为预测遥感影像的蓝光波段，IBI为建成区指数，SI为土壤指数；
Figure BDA0003740459250000101
Figure BDA0003740459250000102
其中，h为普朗克常数，h＝6.26×10-34J-sec；c为光速，c＝2.998×108m/sec；K为斯忒潘-玻尔兹曼常数，K＝1.38×10-23J/K；α为中间变量，λ为热红外波段中心波长，ε为地表出射率；Tb为亮度温度；
(6)、对遥感指数VC、VHI、Wet、NDBSI、LST采用百分位去噪方法去噪，具体为：设定百分位为a，将全部研究时期所包含的像元数值按从小到大排序，取a％作为有效数值下限，(100-a)％为有效数值上限，若像元数值小于有效数值下限，则取有效数值下限；若像元数值大于有效数值上限，则取有效数值上限；然后通过关系式(16)对五个遥感指数进行归一化，以屏蔽不同遥感指数数值范围差异的影响；
Figure BDA0003740459250000103
其中，Index为原始遥感指数，Indexlb和Indexub分别为遥感指数的下限和上限，Indexnom为归一化后的遥感指数；
(7)、将去噪归一化后的各遥感指数通过最小噪声分离变换(MNF)，并取其第一分量作为初步的地表环境综合指数LSECI0，观察第一分量中各遥感指数的载荷数值，若表示植被覆盖度VC的载荷数值大于等于零，则LSECI0为最终的地表环境综合指数LSECI；若表示植被覆盖度VC的载荷数值为负，则将LSECI0取反获得最终的地表环境综合指数LSECI,其关系式为：
LSECI0＝MNF1(VC,VHI,Wet,NDBSI,LST) (17)
Figure BDA0003740459250000104
其中，MNF1为最小噪声分离变换的第一分量，q(VC)为在第一分量中植被覆盖度VC的载荷数值。
以下对本方法进行结果验证。
实例一：
以杭州市为实验区，以陆地资源卫星(Landsat)影像为数据源，生成年际(每年7月1日)的地表环境综合评估结果，具体为：
(1)、获取美国地质调查局USGS公布的L2C2 Landsat影像产品，每景影像已经进行***地几何和辐射校正，收集1985～2020年Landsat TM/ETM+/OLI影像产品共计4348景。每景Landsat影像使用了8个波段，包括6个光学波段(蓝光、绿光、红光、近红外、短波红外1和短波红外2)、1个热红外波段和1个质量评估波段。在质量评估波段中，数值3,4,5分别表示地表被云、雪、云阴影覆盖情况。以此为据，剔除每景影像中的云、雪、云阴影覆盖的像元，再按照成像时间先后顺序依次排列影像，构成影像时间序列；
(2)、对于光学波段和热红外波段，逐像元建立像元的初始时间序列模型：
首先收集24个有效观测数据，利用关系式(3)建立初始时间序列模型，采用最小二乘法确定时序模型的四个参数的取值，见附图1；
(3)、对像元的初始时间序列模型进行变化检测，考虑到蓝光波段波长短、潜在噪声多，本实例中以其它5个光学波段(绿光、红光、近红外、短波红外1和短波红外2)作为模型变化检测波段，因此，变化检测的阈值Vinv_x2＝15.1。依次将后续6个有效观测值作为一组，利用关系式(4)进行变化检测，若检测出变化，将后续6个有效观测值中的第1个有效观测值出现的时间作为变化发生时间，并将该变化发生时间作为该像元的时间序列模型中的断点，对初始时间序列模型进行分段，得到分段时间序列模型；否则，将后续6个有效观测值中的第1个加入初始时间序列模型中，并对像元的初始时间序列模型拟合更新，见附图1，此过程重复进行，直至遍历完所有有效观测值为止；
(4)、获得任意时刻的预测遥感影像：本实例假定1986～2019年每年的7月1日为地表环境综合评估时间，将其转换成对应儒略日，逐波段逐像元定位以上儒略日对应的时间序列模型，通过关系式(1～3)获得像元地表反射率，在此基础上，将全部像元组成波段，将全部波段组成影像，见附图1，通过以上过程，获得1986～2019年每年的7月1日预测遥感影像共34景；
(5)、遥感指数计算：基于34景预测遥感影像逐景计算植被覆盖度VC、植被健康度VHI、土壤湿度Wet、地表土建度NDBSI和地表温度LST，其中，适用于Landsat影像的土壤湿度Wet的关系式为：
WET＝0.1484ρBlue+0.3068ρGreen+0.2437ρRed+0.1886ρNIR-0.7184ρSWIR1-0.5352ρSWIR2(19)
在计算地表温度LST时，利用Globeland30全球土地覆盖产品，见附图1，确定每个像元的土地利用类型，并设定对应的地表出射率ε：0.97-森林，0.91-草地，0.92-人造地表，0.95-裸地，0.99-水体，0.97-其它植被；
(6)、遥感指数去噪和归一化：利用关系式(16)对各遥感指数所有时期的影像进行百分位去噪归一化，百分位数a设置为0.5，即遥感指数像元数值的有效区间为[0.5％,99.5％]，与常规最大最小值归一化方法相比，百分位去噪归一化方法能够更好保证归一化前后数值分布的相对一致性，见附图2；
(7)、地表环境综合指数合成：对于每一时期的各遥感指数，利用关系式(17)进行最小噪声分离变换，获得第一分量作为初步的地表环境综合指数LSECI0，第一分量对于遥感指数整体信息的贡献度为95.3％，说明第一分量能够综合代表各遥感指数表现的地表状况，第一分量中，由于植被覆盖度VC的载荷为-0.632，故利用关系式(18)将初步的地表环境综合指数取反获得最终的地表环境综合指数LSECI。采用以上过程获得的1986～2019年每年的7月1日地表环境综合评估结果见附图3所示，数值越大，代表地表环境质量越好。
实例二：
同样以杭州市为实验区，以陆地资源卫星(Landsat)影像为数据源，生成2019年半月际(每月10日和20日)的地表环境综合评估结果。
其余步骤与实例一相同，不同之处在于：步骤(4)中，本实例假定2019年每月10日和20日为地表环境综合评估时间，将其转换成对应儒略日，逐波段逐像元定位以上儒略日对应的时间序列模型，通过关系式(1～3)获得像元地表反射率，而后将全部像元组成波段，将全部波段组成影像，获得2019年每月10日和20日预测遥感影像共24景；最终获得的2019年半月际地表环境综合评估结果见附图4所示，这进一步证实了本发明能够进行任意时刻的地表环境综合评估。
本发明的保护范围包括但不限于以上实施方式，其保护范围以权利要求书为准，任何对本技术做出的本领域的技术人员容易想到的替换、变形、改进均落入本发明的保护范围。

Claims (8)

1.一种基于时序遥感的任意时刻地表环境综合评估方法，其特征在于包括以下具体步骤：
(1)、对研究时期内选定区域的实拍遥感影像进行预处理，并将预处理后的全部影像按照成像时间顺序依次排列构成影像时间序列；
(2)、对影像中的每一像元根据时间序列在所有波段上分别建立初始时间序列模型；
(3)、对像元的初始时间序列模型进行变化检测，若变化没有发生，则保持初始时间序列模型；若变化发生，则将变化发生时间作为时间序列模型的断点，对像元的初始时间序列模型进行分段，得到分段时间序列模型，以反映对应像元的地表反射率变化；
(4)、在像元的时间序列模型的基础上，获得任意时刻的预测遥感影像；
(5)、基于预测遥感影像，将植被覆盖度VC、植被健康度VHI、土壤湿度Wet、地表土建度NDBSI和地表温度LST五个遥感指数作为选定区域的地表环境综合表征，并计算五个遥感指数的数值；
(6)、将遥感指数去噪和归一化；
(7)、获得地表环境综合指数。
2.如权利要求1所述的一种基于时序遥感的任意时刻地表环境综合评估方法，其特征在于所述的步骤(1)的实拍遥感影像预处理的具体方法为：下载研究时期内选定区域中全部经过几何校正和辐射校正的实拍光学遥感影像数据，或下载全部原始实拍光学遥感影像数据后人工进行几何校正和辐射校正，然后利用遥感影像的质量评估波段，剔除每景影像中被云、云阴影或雪覆盖的无效观测值。
3.如权利要求1所述的一种基于时序遥感的任意时刻地表环境综合评估方法，其特征在于所述的步骤(2)的具体建立方法为：每一像元的初始时间序列模型由线性函数和谐波函数组成，依据每一像元在时间序列中的有效观测值数量，从粗到细构建三种关系式，若12≤有效观测值数量＜18，则模型选用关系式(1)；若18≤有效观测值数量＜24，则模型选用关系式(2)；若有效观测值数量≥24，则模型选用关系式(3)；
Figure FDA0003740459240000011
Figure FDA0003740459240000012
Figure FDA0003740459240000013
其中，
Figure FDA0003740459240000021
表示像元的预测地表反射率；n表示初始时间序列模型所选用的关系式的标号，n＝1,2,3；i为遥感影像第i个波段；t为儒略日；w为频率，
Figure FDA0003740459240000022
c0,i和c1,i为截距系数和斜率系数，an,i和bn,i是第n次谐波系数，且c0,i、c1,i、an,i和bn,i这四个参数均通过最小二乘法确定最优取值。
4.如权利要求3所述的一种基于时序遥感的任意时刻地表环境综合评估方法，其特征在于所述的步骤(3)的具体检测方法为：为了准确检测像元的初始时间序列模型的变化，变化的发生时间通过比较像元的初始时间序列模型预测值和后续有效观测值的残差resid与初始时间序列模型的均方根误差RMSE确定，即
Figure FDA0003740459240000023
其中：D为遥感影像的波段集合，i为遥感影像的第i个波段，residi为第i个波段的初始时间序列模型预测值与后续有效观测值的残差，RMSEi为第i个波段的像元的初始时间序列模型的均方根误差，
Figure FDA0003740459240000024
为阈值，等于自由度为影像波段数的逆卡方分布数；
将同一波段中同一像元的6个连续的后续有效观测值作为一组进行判别，若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值均大于阈值，则将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值出现的时间作为变化发生时间，并将该变化发生时间作为该像元在每个波段中的时间序列模型中的断点，对该像元在每个波段中的初始时间序列模型进行分段，并根据步骤(2)中的关系式得到分段时间序列模型，以反映该像元的地表反射率变化；
若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值没有全部大于阈值，则认为变化没有发生，并将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值加入对应波段的初始时间序列模型中，然后根据步骤(2)中的关系式将像元的初始时间序列模型拟合更新；
上述判别过程重复进行，直到全部有效观测值都位于像元的初始时间序列模型或分段时间序列模型中为止。
5.如权利要求3所述的一种基于时序遥感的任意时刻地表环境综合评估方法，其特征在于所述的步骤(4)的具体方法为：逐波段逐像元定位像元的初始时间序列模型或分段时间序列模型，根据研究的需要确定地表环境综合评估的儒略日t，再通过步骤(2)的关系式得到像元的预测地表反射率，然后将全部像元组成波段，将全部波段组成影像，获得任意时刻的预测遥感影像。
6.如权利要求1所述的一种基于时序遥感的任意时刻地表环境综合评估方法，其特征在于所述的步骤(5)中的五个遥感指数：植被覆盖度VC、植被健康度VHI、土壤湿度Wet、地表土建度NDBSI和地表温度LST的关系式为：
Figure FDA0003740459240000031
Figure FDA0003740459240000032
其中，ρNIR和ρRed分别为预测遥感影像的近红外波段和红光波段，NDVI为归一化差异植被指数，NDVIveg和NDVIsoil分别为植被和土壤的归一化差异植被指数，两者分别为研究范围内NDVI的最大和最小值；
Figure FDA0003740459240000033
Figure FDA0003740459240000034
VHI＝PCA1(NDVI,NRI,NDSVI) (9)
其中，ρSWIR1和ρGreen分别为预测遥感影像的短波红外第1波段和绿光波段，NRI为氮元素的反射率指数，NDSVI为归一化差异植被衰亡指数，植被健康度VHI为通过将NDVI、NRI、NDSVI主成分变换后的第一分量PCA1；
Figure FDA0003740459240000035
其中，D是遥感影像波段集合，i为预测遥感影像第i个波段，ai为波段线性组合系数，ρi为第i个波段的影像地表反射率；
Figure FDA0003740459240000036
Figure FDA0003740459240000037
Figure FDA0003740459240000041
其中，ρBlue为预测遥感影像的蓝光波段，IBI为建成区指数，SI为土壤指数；
Figure FDA0003740459240000042
Figure FDA0003740459240000043
其中，h为普朗克常数，h＝6.26×10-34J-sec；c为光速，c＝2.998×108m/sec；K为斯忒潘-玻尔兹曼常数，K＝1.38×10-23J/K；α为中间变量，λ为热红外波段中心波长，ε为地表出射率；Tb为亮度温度。
7.如权利要求1所述的一种基于时序遥感的任意时刻地表环境综合评估方法，其特征在于所述的步骤(6)的具体方法为：
对遥感指数VC、VHI、Wet、NDBSI、LST采用百分位去噪方法，设定百分位为a，将全部研究时期所包含的像元数值按从小到大排序，取a％作为有效数值下限，(100-a)％为有效数值上限，若像元数值小于有效数值下限，则取有效数值下限；若像元数值大于有效数值上限，则取有效数值上限；然后通过关系式(16)对五个遥感指数进行归一化，以屏蔽不同遥感指数数值范围差异的影响；
Figure FDA0003740459240000044
其中，Index为原始遥感指数，Indexlb和Indexub分别为遥感指数的下限和上限，Indexnom为归一化后的遥感指数。
8.如权利要求1所述的一种基于时序遥感的任意时刻地表环境综合评估方法，其特征在于所述的步骤(7)的具体方法为：
将去噪归一化后的各遥感指数通过最小噪声分离变换，并取其第一分量作为初步的地表环境综合指数LSECI0，观察第一分量中各遥感指数的载荷数值，若表示植被覆盖度VC的载荷数值大于等于零，则LSECI0为最终的地表环境综合指数LSECI；若表示植被覆盖度VC的载荷数值为负，则将LSECI0取反获得最终的地表环境综合指数LSECI,其关系式为：
LSECI0＝MNF1(VC,VHI,Wet,NDBSI,LST) (17)
Figure FDA0003740459240000051
其中，MNF1为最小噪声分离变换的第一分量，q(VC)为在第一分量中植被覆盖度VC的载荷数值。
CN202210810131.0A 2022-07-11 2022-07-11 一种基于时序遥感的任意时刻地表环境综合评估方法 Active CN115203624B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210810131.0A CN115203624B (zh) 2022-07-11 2022-07-11 一种基于时序遥感的任意时刻地表环境综合评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210810131.0A CN115203624B (zh) 2022-07-11 2022-07-11 一种基于时序遥感的任意时刻地表环境综合评估方法

Publications (2)

Publication Number Publication Date
CN115203624A true CN115203624A (zh) 2022-10-18
CN115203624B CN115203624B (zh) 2023-06-16

Family

ID=83579312

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210810131.0A Active CN115203624B (zh) 2022-07-11 2022-07-11 一种基于时序遥感的任意时刻地表环境综合评估方法

Country Status (1)

Country Link
CN (1) CN115203624B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116466368A (zh) * 2023-06-16 2023-07-21 成都远望科技有限责任公司 基于激光雷达和卫星资料的沙尘消光系数廓线估算方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200311842A1 (en) * 2019-03-28 2020-10-01 China Waterborne Transport Research Institute Method for tracking, monitoring and evaluating ecological impact of channel project based on long-term time series satellite remote sensing data
US20210027429A1 (en) * 2019-07-26 2021-01-28 Zhejiang University Of Technology Noise detection method for time-series vegetation index derived from remote sensing images
CN113988626A (zh) * 2021-10-28 2022-01-28 中国矿业大学(北京) 一种矿区生态环境遥感综合评价指数实现方法
CN114241334A (zh) * 2021-12-20 2022-03-25 南京大学 一种多因子集成的干旱区生态环境遥感监测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200311842A1 (en) * 2019-03-28 2020-10-01 China Waterborne Transport Research Institute Method for tracking, monitoring and evaluating ecological impact of channel project based on long-term time series satellite remote sensing data
US20210027429A1 (en) * 2019-07-26 2021-01-28 Zhejiang University Of Technology Noise detection method for time-series vegetation index derived from remote sensing images
CN113988626A (zh) * 2021-10-28 2022-01-28 中国矿业大学(北京) 一种矿区生态环境遥感综合评价指数实现方法
CN114241334A (zh) * 2021-12-20 2022-03-25 南京大学 一种多因子集成的干旱区生态环境遥感监测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
李红星 等: "基于遥感生态指数的武汉市生态环境质量评估", 云南大学学报(自然科学版) *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116466368A (zh) * 2023-06-16 2023-07-21 成都远望科技有限责任公司 基于激光雷达和卫星资料的沙尘消光系数廓线估算方法
CN116466368B (zh) * 2023-06-16 2023-08-22 成都远望科技有限责任公司 基于激光雷达和卫星资料的沙尘消光系数廓线估算方法

Also Published As

Publication number Publication date
CN115203624B (zh) 2023-06-16

Similar Documents

Publication Publication Date Title
Loiseau et al. Satellite data integration for soil clay content modelling at a national scale
Song et al. Predicting temperate conifer forest successional stage distributions with multitemporal Landsat Thematic Mapper imagery
CN112395808A (zh) 一种结合随机森林和协同克里金的生物量遥感制图方法
CN110927120B (zh) 一种植被覆盖度预警方法
CN112733596A (zh) 一种基于中高空间分辨率遥感影像融合的森林资源变化监测方法及应用
CN114821349B (zh) 顾及谐波模型系数和物候参数的森林生物量估算方法
CN115372282B (zh) 一种基于无人机高光谱影像的农田土壤含水量监测方法
CN114241331B (zh) 以UAV为地面和Sentinel-2中介的湿地芦苇地上生物量遥感建模方法
CN114819737B (zh) 公路路域植被的碳储量估算方法、***及存储介质
Fitzgerald et al. Directed sampling using remote sensing with a response surface sampling design for site-specific agriculture
CN115203624B (zh) 一种基于时序遥感的任意时刻地表环境综合评估方法
CN112215135B (zh) 矿区开采与治理成效监测方法及装置
Khudhur et al. Comparison of the accuracies of different spectral indices for mapping the vegetation covers in Al-Hawija district, Iraq
Singh et al. Environmental degradation analysis using NOAA/AVHRR data
Das et al. Mapping vegetation and forest types using Landsat TM in the western ghat region of Maharashtra, India
Fernandes et al. A multi-scale approach to mapping effective leaf area index in boreal Picea mariana stands using high spatial resolution CASI imagery
Salman et al. Detection of Spectral Reflective Changes for Temporal Resolution of Land Cover (LC) for Two Different Seasons in central Iraq
Nielsen et al. The distribution in time and space of savanna fires in Burkina Faso as determined from NOAA AVHRR data
Colaninno et al. Defining Densities for Urban Residential Texture, through Land use Classification, from Landsat TM Imagery: Case Study of Spanish Mediterranean Coast
CN117196922A (zh) 一种连续评估自然保护区生态质量的遥感方法
Zhao et al. Towards an Annually Updated Oil Palm Age Database (OPAD) with Multimodal Remote Sensing Features: Design and Preliminary Result
Gupta et al. Fraction of vegetation cover and its application in vegetation characterization in the Hazaribagh wildlife sanctuary, Jharkhand, India
Singh et al. Estimation of changes in land surface temperature using multi-temporal Landsat data of Ghaziabad District, India.
Wang Monitoring time-series changes in vegetation coverage based on landsat images
Blasch Multitemporal soil pattern analysis for organic matter estimation at croplands using multispectral satellite data

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