CN115203624A - 一种基于时序遥感的任意时刻地表环境综合评估方法 - Google Patents
一种基于时序遥感的任意时刻地表环境综合评估方法 Download PDFInfo
- 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
Links
- 238000011156 evaluation Methods 0.000 title claims abstract description 36
- 238000000034 method Methods 0.000 claims abstract description 46
- 230000008859 change Effects 0.000 claims abstract description 45
- 238000011160 research Methods 0.000 claims abstract description 18
- 238000003384 imaging method Methods 0.000 claims abstract description 7
- 238000007781 pre-processing Methods 0.000 claims abstract description 5
- 238000002310 reflectometry Methods 0.000 claims description 20
- 239000002689 soil Substances 0.000 claims description 17
- 238000012937 correction Methods 0.000 claims description 12
- 230000036541 health Effects 0.000 claims description 10
- 230000003287 optical effect Effects 0.000 claims description 10
- 230000008569 process Effects 0.000 claims description 10
- 238000000926 separation method Methods 0.000 claims description 8
- 230000005855 radiation Effects 0.000 claims description 7
- IJGRMHOSHXDMSA-UHFFFAOYSA-N Atomic nitrogen Chemical compound N#N IJGRMHOSHXDMSA-UHFFFAOYSA-N 0.000 claims description 6
- 101000812677 Homo sapiens Nucleotide pyrophosphatase Proteins 0.000 claims description 6
- 101000841763 Homo sapiens Ubiquinol-cytochrome-c reductase complex assembly factor 2 Proteins 0.000 claims description 6
- 102100039306 Nucleotide pyrophosphatase Human genes 0.000 claims description 6
- 102100029513 Ubiquinol-cytochrome-c reductase complex assembly factor 2 Human genes 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 5
- 238000013441 quality evaluation Methods 0.000 claims description 4
- 238000012886 linear function Methods 0.000 claims description 3
- 229910052757 nitrogen Inorganic materials 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000011895 specific detection Methods 0.000 claims description 2
- 230000000737 periodic effect Effects 0.000 abstract description 6
- 238000001514 detection method Methods 0.000 description 7
- 238000010606 normalization Methods 0.000 description 5
- 230000007613 environmental effect Effects 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000003915 air pollution Methods 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000005184 irreversible process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000001303 quality assessment method Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
- 238000003911 water pollution Methods 0.000 description 1
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);
其中,表示像元的预测地表反射率;n表示初始时间序列模型所选用的关系式的标号,n=1,2,3;i为遥感影像第i个波段;t为儒略日;w为频率,c0,i和c1,i为截距系数和斜率系数,an,i和bn,i是第n次谐波系数,且c0,i、c1,i、an,i和bn,i这四个参数均通过最小二乘法确定最优取值。
进一步地,所述的步骤(3)的具体检测方法为:为了准确检测像元的初始时间序列模型的变化,变化的发生时间通过比较像元的初始时间序列模型预测值和后续有效观测值的残差resid与初始时间序列模型的均方根误差RMSE确定,即
其中:D为遥感影像的波段集合,i为遥感影像的第i个波段,residi为第i个波段的初始时间序列模型预测值与后续有效观测值的残差,RMSEi为第i个波段的像元的初始时间序列模型的均方根误差,为阈值,等于自由度为影像波段数的逆卡方分布数;
将同一波段中同一像元的6个连续的后续有效观测值作为一组进行判别,若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值均大于阈值,则将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值出现的时间作为变化发生时间,并将该变化发生时间作为该像元在每个波段中的时间序列模型中的断点,对该像元在每个波段中的初始时间序列模型进行分段,并根据步骤(2)中的关系式得到分段时间序列模型,以反映该像元的地表反射率变化;
若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值没有全部大于阈值,则认为变化没有发生,并将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值加入对应波段的初始时间序列模型中,然后根据步骤(2)中的关系式将像元的初始时间序列模型拟合更新;
上述判别过程重复进行,直到全部有效观测值都位于像元的初始时间序列模型或分段时间序列模型中为止。
进一步地,所述的步骤(4)的具体方法为:逐波段逐像元定位像元的初始时间序列模型或分段时间序列模型,根据需要确定地表环境综合评估的儒略日t,再通过步骤(2)的关系式得到像元的预测地表反射率,然后将全部像元组成波段,将全部波段组成影像,获得任意时刻的预测遥感影像。
进一步地,所述的步骤(5)中的五个遥感指数:植被覆盖度VC、植被健康度VHI、土壤湿度Wet、地表土建度NDBSI和地表温度LST的关系式为:
其中,ρNIR和ρRed分别为预测遥感影像的近红外波段和红光波段,NDVI为归一化差异植被指数,NDVIveg和NDVIsoil分别为植被和土壤的归一化差异植被指数,两者分别为研究范围内NDVI的最大和最小值;
VHI=PCA1(NDVI,NRI,NDSVI) (9)
其中,ρSWIR1和ρGreen分别为预测遥感影像的短波红外第1波段和绿光波段,NRI为氮元素的反射率指数,NDSVI为归一化差异植被衰亡指数,植被健康度VHI为通过将NDVI、NRI、NDSVI主成分变换后的第一分量PCA1;
其中,D是遥感影像波段集合,i为预测遥感影像第i个波段,ai为波段线性组合系数,ρi为第i个波段的影像地表反射率;
其中,ρBlue为预测遥感影像的蓝光波段,IBI为建成区指数,SI为土壤指数;
其中,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)对五个遥感指数进行归一化,以屏蔽不同遥感指数数值范围差异的影响;
其中,Index为原始遥感指数,Indexlb和Indexub分别为遥感指数的下限和上限,Indexnom为归一化后的遥感指数。
进一步地,所述的步骤(7)的具体方法为:
将去噪归一化后的各遥感指数通过最小噪声分离变换,并取其第一分量作为初步的地表环境综合指数LSECI0,观察第一分量中各遥感指数的载荷数值,若表示植被覆盖度VC的载荷数值大于等于零,则LSECI0为最终的地表环境综合指数LSECI;若表示植被覆盖度VC的载荷数值为负,则将LSECI0取反获得最终的地表环境综合指数LSECI,其关系式为:
LSECI0=MNF1(VC,VHI,Wet,NDBSI,LST) (17)
其中,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);
其中,表示像元的预测地表反射率;n表示初始时间序列模型所选用的关系式的标号,n=1,2,3;i为遥感影像第i个波段;t为儒略日;w为频率,c0,i和c1,i为截距系数和斜率系数,代表逐渐的年际变化;an,i和bn,i是第n次谐波系数,代表周期性的年内变化,且c0,i、c1,i、an,i和bn,i这四个参数均通过最小二乘法确定最优取值;
(3)、对像元的初始时间序列模型进行变化检测:
为了准确检测像元的初始时间序列模型的变化,变化的发生时间通过比较像元的初始时间序列模型预测值和后续有效观测值的残差resid与初始时间序列模型的均方根误差RMSE确定,即
其中:D为遥感影像的波段集合,i为遥感影像的第i个波段,residi为第i个波段的初始时间序列模型预测值与后续有效观测值的残差,RMSEi为第i个波段的像元的初始时间序列模型的均方根误差,为阈值,等于自由度为影像波段数的逆卡方分布数;
将同一波段中同一像元的6个连续的后续有效观测值作为一组进行判别,若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值均大于阈值,则将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值出现的时间作为变化发生时间,并将该变化发生时间作为该像元在每个波段中的时间序列模型中的断点,对该像元在每个波段中的初始时间序列模型进行分段,并根据步骤(2)中的关系式得到分段时间序列模型,以反映该像元的地表反射率变化;
若所有波段中同一像元的6个后续有效观测值的resid与RMSE的比值没有全部大于阈值,则认为变化没有发生,并将该像元在每个波段中的后续6个有效观测值中的第1个有效观测值加入对应波段的初始时间序列模型中,然后根据步骤(2)中的关系式将像元的初始时间序列模型拟合更新;
上述判别过程重复进行,直到全部有效观测值都位于像元的初始时间序列模型或分段时间序列模型中为止;
(4)、在像元的时间序列模型的基础上,获得任意时刻的预测遥感影像,具体为:
逐波段逐像元定位像元的初始时间序列模型或分段时间序列模型,根据研究的需要确定地表环境综合评估的儒略日t,t可为研究时期内的任意时间,再通过步骤(2)的关系式得到像元的预测地表反射率,然后将全部像元组成波段,将全部波段组成影像,获得任意时刻的预测遥感影像;
(5)、基于预测遥感影像,将植被覆盖度VC、植被健康度VHI、土壤湿度Wet、地表土建度NDBSI和地表温度LST五个遥感指数作为选定区域的地表环境综合表征,并计算五个遥感指数的数值,具体关系式为:
其中,ρNIR和ρRed分别为预测遥感影像的近红外波段和红光波段,NDVI为归一化差异植被指数,NDVIveg和NDVIsoil分别为植被和土壤的归一化差异植被指数,两者分别为研究范围内NDVI的最大和最小值;
VHI=PCA1(NDVI,NRI,NDSVI) (9)
其中,ρSWIR1和ρGreen分别为预测遥感影像的短波红外第1波段和绿光波段,NRI为氮元素的反射率指数,NDSVI为归一化差异植被衰亡指数,植被健康度VHI为通过将NDVI、NRI、NDSVI主成分变换后的第一分量PCA1;
其中,D是遥感影像波段集合,i为预测遥感影像第i个波段,ai为波段线性组合系数,ρi为第i个波段的影像地表反射率;
其中,ρBlue为预测遥感影像的蓝光波段,IBI为建成区指数,SI为土壤指数;
其中,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)对五个遥感指数进行归一化,以屏蔽不同遥感指数数值范围差异的影响;
其中,Index为原始遥感指数,Indexlb和Indexub分别为遥感指数的下限和上限,Indexnom为归一化后的遥感指数;
(7)、将去噪归一化后的各遥感指数通过最小噪声分离变换(MNF),并取其第一分量作为初步的地表环境综合指数LSECI0,观察第一分量中各遥感指数的载荷数值,若表示植被覆盖度VC的载荷数值大于等于零,则LSECI0为最终的地表环境综合指数LSECI;若表示植被覆盖度VC的载荷数值为负,则将LSECI0取反获得最终的地表环境综合指数LSECI,其关系式为:
LSECI0=MNF1(VC,VHI,Wet,NDBSI,LST) (17)
其中,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);
4.如权利要求3所述的一种基于时序遥感的任意时刻地表环境综合评估方法,其特征在于所述的步骤(3)的具体检测方法为:为了准确检测像元的初始时间序列模型的变化,变化的发生时间通过比较像元的初始时间序列模型预测值和后续有效观测值的残差resid与初始时间序列模型的均方根误差RMSE确定,即
其中:D为遥感影像的波段集合,i为遥感影像的第i个波段,residi为第i个波段的初始时间序列模型预测值与后续有效观测值的残差,RMSEi为第i个波段的像元的初始时间序列模型的均方根误差,为阈值,等于自由度为影像波段数的逆卡方分布数;
将同一波段中同一像元的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的关系式为:
其中,ρNIR和ρRed分别为预测遥感影像的近红外波段和红光波段,NDVI为归一化差异植被指数,NDVIveg和NDVIsoil分别为植被和土壤的归一化差异植被指数,两者分别为研究范围内NDVI的最大和最小值;
VHI=PCA1(NDVI,NRI,NDSVI) (9)
其中,ρSWIR1和ρGreen分别为预测遥感影像的短波红外第1波段和绿光波段,NRI为氮元素的反射率指数,NDSVI为归一化差异植被衰亡指数,植被健康度VHI为通过将NDVI、NRI、NDSVI主成分变换后的第一分量PCA1;
其中,D是遥感影像波段集合,i为预测遥感影像第i个波段,ai为波段线性组合系数,ρi为第i个波段的影像地表反射率;
其中,ρBlue为预测遥感影像的蓝光波段,IBI为建成区指数,SI为土壤指数;
其中,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)对五个遥感指数进行归一化,以屏蔽不同遥感指数数值范围差异的影响;
其中,Index为原始遥感指数,Indexlb和Indexub分别为遥感指数的下限和上限,Indexnom为归一化后的遥感指数。
8.如权利要求1所述的一种基于时序遥感的任意时刻地表环境综合评估方法,其特征在于所述的步骤(7)的具体方法为:
将去噪归一化后的各遥感指数通过最小噪声分离变换,并取其第一分量作为初步的地表环境综合指数LSECI0,观察第一分量中各遥感指数的载荷数值,若表示植被覆盖度VC的载荷数值大于等于零,则LSECI0为最终的地表环境综合指数LSECI;若表示植被覆盖度VC的载荷数值为负,则将LSECI0取反获得最终的地表环境综合指数LSECI,其关系式为:
LSECI0=MNF1(VC,VHI,Wet,NDBSI,LST) (17)
其中,MNF1为最小噪声分离变换的第一分量,q(VC)为在第一分量中植被覆盖度VC的载荷数值。
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)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116466368A (zh) * | 2023-06-16 | 2023-07-21 | 成都远望科技有限责任公司 | 基于激光雷达和卫星资料的沙尘消光系数廓线估算方法 |
Citations (4)
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 | 南京大学 | 一种多因子集成的干旱区生态环境遥感监测方法 |
-
2022
- 2022-07-11 CN CN202210810131.0A patent/CN115203624B/zh active Active
Patent Citations (4)
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)
Title |
---|
李红星 等: "基于遥感生态指数的武汉市生态环境质量评估", 云南大学学报(自然科学版) * |
Cited By (2)
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 |