CN115063332A - 一种构建高空间分辨率时序遥感数据的方法 - Google Patents

一种构建高空间分辨率时序遥感数据的方法 Download PDF

Info

Publication number
CN115063332A
CN115063332A CN202210748388.8A CN202210748388A CN115063332A CN 115063332 A CN115063332 A CN 115063332A CN 202210748388 A CN202210748388 A CN 202210748388A CN 115063332 A CN115063332 A CN 115063332A
Authority
CN
China
Prior art keywords
image
spatial
reflectivity
pixel
resolution
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
CN202210748388.8A
Other languages
English (en)
Other versions
CN115063332B (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.)
Hebei Normal University of Science and Technology
Original Assignee
Hebei Normal University of Science and Technology
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 Hebei Normal University of Science and Technology filed Critical Hebei Normal University of Science and Technology
Priority to CN202210748388.8A priority Critical patent/CN115063332B/zh
Publication of CN115063332A publication Critical patent/CN115063332A/zh
Application granted granted Critical
Publication of CN115063332B publication Critical patent/CN115063332B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/50Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/20Image preprocessing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/761Proximity, similarity or dissimilarity measures
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00Arrangements for image or video recognition or understanding
    • G06V10/70Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/764Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06VIMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00Scenes; Scene-specific elements
    • G06V20/10Terrestrial scenes
    • G06V20/13Satellite images
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10032Satellite or aerial image; Remote sensing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20172Image enhancement details
    • G06T2207/20182Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20212Image combination
    • G06T2207/20221Image fusion; Image merging
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30181Earth observation
    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Databases & Information Systems (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Astronomy & Astrophysics (AREA)
  • Remote Sensing (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及一种构建高空间分辨率时序遥感数据的方法,具体包含以下步骤:S1、低空间分辨率影像的混合像元分解;S2、基准期与预测期的回归模型拟合;S3、残差的分配和补偿;S4、拟合结果的空间滤波;S5、正反向估测结果的加权;S6、结合高空间分辨率影像的连续校正。本发明取得了更高的融合精度,极大地减弱了空间模糊现象,提供了更丰富的纹理特征和更好的目视效果,有利于地表信息连续动态监测及其他高空间分辨率时序遥感数据的应用。

Description

一种构建高空间分辨率时序遥感数据的方法
技术领域
本发明涉及遥感影像处理技术领域,尤其涉及一种构建高空间分辨率时序遥感数据的方法。
背景技术
遥感技术以其空间覆盖范围广、时间周期固定等优势为区域地表信息动态监测提供了重要的数据来源。然而,受到传感器本身和各种不利天气条件的限制,目前单一传感器很难兼具空间和时间分辨率。高时间分辨率的影像其空间分辨率较低,混合像元的影响制约了地表信息监测的精度;而高空间分辨率的影像回访周期长、时效性差,难以满足连续的动态监测。遥感数据时空融合技术是满足高时空应用需求的有效方式。遥感影像融合是指利用一对或多对准同步的粗-细分辨率影像作为基准期影像,结合待预测日期的粗分辨率影像估测待预测日期的细分辨率影像。
针对特定的应用场景学者们提出了诸多的时空融合算法,其中Fit-F方法仅仅需要起始日期的一对粗-细分辨率的基准期影像,并且计算效率高,对于有物候变化的场景表现出较好的预测精度,因此适用于地表信息连续的动态监测。但是该方法存在两方面的不足:一方面,受粗分辨率影像上混合像元的影响,该方法在异质性强但变化强度低的区域表现不佳,不擅长捕捉影像的结构和纹理;另一方面,由于仅利用一对基准期影像,对于基准期和预
说明书
测期的细分辨率影像间存在变化,而粗分辨率影像上不可见的情况,不能被捕捉到。因此,有必要改进现有的Fit-FC方法,使其更加适应于连续的地表信息动态监测。
此外,Sentinel-2A/B卫星具有5天的重访周期,即使受云等天气条件的影响造成数据缺失,其仍然为地表信息监测提供了较为密集的高空间分辨率影像。如何充分利用可获得的多幅高空间分辨率影像提供的信息,构建出更高质量的高空间分辨率时序遥感数据,也是亟待解决的问题。
发明内容
为了克服上述现有技术的局限,本发明提供了一种构建高空间分辨率时序遥感数据的方法。
本发明的技术方案如下:
一种构建高空间分辨率时序遥感数据的方法,具体包含以下步骤:
S1、低空间分辨率影像的混合像元分解:基于高空间分辨率影像的分类图,获取低空间分辨率影像上逐像元的各类地物的丰度,并利用像元分解模型计算每类地物的光谱反射率,获得粗影像的降尺度影像;
S2、基准期与预测期的回归模型拟合:基于步骤S1中得到的降尺度影像,建立基准期和预测期(
Figure 73132DEST_PATH_IMAGE001
)间反射率的线性回归模型,并应用于细影像,得到
Figure 459114DEST_PATH_IMAGE001
细影像的反射率拟合结果;
说明书
S3、残差的分配和补偿:将上述步骤S2得到的
Figure 97905DEST_PATH_IMAGE001
细影像的反射率升尺度至低空间分辨率,计算其与粗影像反射率的差值,即残差,并采用空间插值方法将残差分配至高空间分辨率,并补偿于
Figure 961956DEST_PATH_IMAGE001
的细影像;
S4、拟合结果的空间滤波:在一定窗口内筛选中心像元的光谱相似像元,利用相似性像元的光谱差异性和与中心像元的相对距离来加权计算预测期的中心像元值,从而得到空间滤波后的
Figure 369804DEST_PATH_IMAGE001
的细影像;
S5、正反向估测结果的加权:设置起始日期(
Figure DEST_PATH_IMAGE002
)和结束日期(
Figure 457846DEST_PATH_IMAGE001
)两个基准期,将上述步骤S1-S4分正、反两个方向进行,采用线性加权方法综合正、反向估测结果作为加权时序融合结果;
S6、结合高空间分辨率影像的连续校正:基于上述步骤S5获得的时序影像融合结果,构建地表信息的动态变化模型;利用可获得的多幅已知细影像对动态模型进行修正,使动态模型模拟值与已知细影像值的误差最小,最终得到校正后的时空融合结果。
优选的,步骤S1中具体包括:
S11:对获得的遥感影像进行大气校正、几何校正、投影转换、重采样、质量控制等预处理;
S12:采用监督分类或非监督分类法对细影像进行分类,利用公式(1)计算粗影像上每个像元内各类地物的丰度;
Figure DEST_PATH_IMAGE003
说明书
式(1)中,
Figure DEST_PATH_IMAGE004
表示i位置处混合像元内类别q的丰度,Q为混合像元中类别q的细像元数,S为混合像元中所有类别的细像元总数;
S13:混合像元的反射率值可以表示为像元中各类别的反射率与其丰度的线性组合,混合像元分解模型具体为:
Figure 162364DEST_PATH_IMAGE005
约束条件:
Figure DEST_PATH_IMAGE006
式(2)中,
Figure 423581DEST_PATH_IMAGE007
为第i个混合像元的反射率,
Figure DEST_PATH_IMAGE008
表示第i个混合像元中类别q的平均反射率,m为类别数,
Figure 358039DEST_PATH_IMAGE009
为残余误差。
S14:将上述步骤计算得到的各类别的平均反射率按照类别对应到相应地类的像元上,从而得到粗影像的降尺度数据;
优选的,步骤S2中具体包括:
根据公式(3)利用基准期和预测期的降尺度数据建立线性回归模型,需要在一定大小的移动窗口M内进行,利用最小二乘法求得移动窗口内中心像元的回归系数a和b,并应用于细影像,根据公式(4)估测得到
Figure 711922DEST_PATH_IMAGE001
的细影像;
Figure DEST_PATH_IMAGE010
Figure 528568DEST_PATH_IMAGE011
说明书
其中,
Figure DEST_PATH_IMAGE012
Figure 531159DEST_PATH_IMAGE013
分别为移动窗口M内
Figure 648020DEST_PATH_IMAGE002
Figure 546706DEST_PATH_IMAGE001
时期降尺度数据的反射率,
Figure DEST_PATH_IMAGE014
为拟合残余误差,
Figure 614762DEST_PATH_IMAGE015
Figure DEST_PATH_IMAGE016
分别表示
Figure 952203DEST_PATH_IMAGE002
Figure 392411DEST_PATH_IMAGE001
时期细影像的反射率。
优选的,步骤S3中具体包括:
采用平均聚合的方法将
Figure 321053DEST_PATH_IMAGE001
的细影像升尺度至低空间分辨率,利用公式(5)计算残差RD;将残差分配至高空间分辨率采用的空间插值方法可以选用适当的克里金插值法,并将高空间分辨率的残差rd补偿于与上述步骤S2得到的
Figure 518816DEST_PATH_IMAGE001
细影像的反射率;
Figure 958150DEST_PATH_IMAGE017
式(5)中,
Figure DEST_PATH_IMAGE018
表示低空间分辨率下的残差,
Figure 518444DEST_PATH_IMAGE019
Figure DEST_PATH_IMAGE020
分别为
Figure 555671DEST_PATH_IMAGE002
Figure 99784DEST_PATH_IMAGE001
时期粗像元的反射率,
Figure DEST_PATH_IMAGE021
Figure DEST_PATH_IMAGE022
分别为和
Figure 283422DEST_PATH_IMAGE001
时期粗像元中细像元的观测和估测反射率。
优选的,步骤S4中具体包括:
根据公式(6),综合基准期细影像的各波段相似像元和中心像元的光谱差异D,在一定窗口内,差异最小的前s个像元被筛选为相似像元;根据公式(7)计算邻近像元和中心像元的空间距离
Figure 229381DEST_PATH_IMAGE023
,根据公式(8)确定各相似像元的权重
Figure DEST_PATH_IMAGE024
;根据公式(9)通过加权窗口内相似性像元逐波段计算预测时期
Figure 437509DEST_PATH_IMAGE001
的中心像元值,从而
说明书
得到
Figure 468919DEST_PATH_IMAGE001
时期空间滤波后的细影像;
Figure 623956DEST_PATH_IMAGE025
Figure DEST_PATH_IMAGE026
Figure 988204DEST_PATH_IMAGE027
Figure DEST_PATH_IMAGE028
其中,
Figure 163970DEST_PATH_IMAGE029
Figure DEST_PATH_IMAGE030
分别表示细影像第b波段的i位置和中心位置像元的反射率,c为用于融合的影像波段数;w为窗口半径,
Figure 417097DEST_PATH_IMAGE031
为高空间分辨率下的残差分配值,
Figure DEST_PATH_IMAGE032
为w×w窗口内中心像元的正向预测结果。
优选的,步骤S5中具体包括:
依据
Figure 405519DEST_PATH_IMAGE002
细影像的分类图,确定正反向估测结果中每一地类的权重;具体地,以正向估测过程为例,根据公式(10),由
Figure 60492DEST_PATH_IMAGE002
估测
Figure 79263DEST_PATH_IMAGE033
的细影像,确定估测影像与实际观测影像间的回归系数k和l,并应用于
Figure 85265DEST_PATH_IMAGE001
的正向估测影像,根据公式(11)得到
Figure 910002DEST_PATH_IMAGE001
的似“真实影像”;根据公式(12)和(13)计算
Figure 560426DEST_PATH_IMAGE001
的正向估测影像和似“真实”影像的差距,得到正向估测的相对精度;
Figure DEST_PATH_IMAGE034
Figure 313881DEST_PATH_IMAGE035
说明书
Figure DEST_PATH_IMAGE036
Figure 541600DEST_PATH_IMAGE037
其中,
Figure DEST_PATH_IMAGE038
Figure 497923DEST_PATH_IMAGE039
分别表示
Figure 268433DEST_PATH_IMAGE033
的正向估测和真实观测值,
Figure DEST_PATH_IMAGE040
Figure 449579DEST_PATH_IMAGE041
分别表示
Figure 430174DEST_PATH_IMAGE001
的正向估测和似“真实”观测值,
Figure DEST_PATH_IMAGE042
表示正向估测与真实观测值之间的差距,
Figure 862292DEST_PATH_IMAGE043
表示正向估测的相对精度;
相似地,依照上述流程,可以获得反向估测的相对精度
Figure DEST_PATH_IMAGE044
;根据公式(14)和(15)确定各地类正反向估测的权重;根据公式(16)获得正反向估测加权的时序影像融合结果;
Figure 18467DEST_PATH_IMAGE045
Figure DEST_PATH_IMAGE046
Figure 644883DEST_PATH_IMAGE047
优选的,步骤S6中具体包括:
采用S-G滤波方法对时序融合数据进行滤波用于构建地表信息动态模型;采用一种连续校正的CACAO方法(公式(17)),利用多幅已知的细影像修正动态模型,将修正后的各波段反射率作为最终的时序数据融合结果;
Figure DEST_PATH_IMAGE048
式(17)中,BPture表示已知的细影像的各波段反射率,n为可获得的细影像的数量,BPpm表示动态模型模拟的各波段反射率,
说明书
scale和shift为需要计算的修正参数。
本发明采用以上技术方案,与现有技术相比,具有以下技术效果:
1、本发明结合了混合像元分解模型,极大地减弱了粗影像上混合像元对融合结果造成的空间模糊现象,得到的融合结果具有更丰富的纹理特征和更好的目视效果。
2、本发明同时考虑了基准期和预测期间拟合回归模型的残差和不同传感器间的误差,并且采用空间插值方法代替传统的插值法来分配残差,取得了更高的融合精度。
3、本发明综合了起始日期和结束日期对预测期的双向估测结果,有利于捕捉起始期和预测期间发生的地表覆盖类型变化,使得融合结果与实际影像更接近。
4、本发明充分利用了可获得的多幅已知高分辨率影像,降低了融合结果与实际影像间的误差,取得更好的时序融合结果。
附图说明
图1为本发明实施例的实现流程图;
图2为可获得的Sentinel-2影像的时间分布图;
图3不同方法融合影像反射率与实际观测影像反射率的散点密度图(以2019年8月4日为例);
图4为不同方法融合影像的空间特征图(以2019年8月4日为例)。
说明书。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明的具体实施方式作进一步详细描述。显然,以下实施例用于说明本发明,但不用来限制本发明的范围。
一、实施案例
如图1所示,一种构建高空间分辨率时序遥感数据的方法,具体包括以下步骤:
S1、低空间分辨率影像的混合像元分解;采用监督分类或非监督分类法对细影像进行分类,基于细影像的分类图,获取粗影像上逐像元的各类地物的丰度;利用像元分解模型计算每类地物的光谱反射率,获得粗影像的降尺度影像;
S11、根据获取到的遥感数据级别情况,进行大气校正、几何校正、投影转换、重采样、质量控制等预处理;
S12、采用监督分类或非监督分类法对基准期的细影像进行分类,根据研究区地表异质性确定分类的类别数;利用公式(1)计算粗影像上每个像元内各类地物的丰度。
Figure 643931DEST_PATH_IMAGE003
式(1)中,
Figure 286265DEST_PATH_IMAGE004
表示i位置处混合像元内类别q的丰度,Q为混合像元中类别q的细像元数,S为混合像元中所有类别的细像元总数;
说明书
S13、忽略不同传感器之间的地理配准和大气校正误差,混合像元的反射率值可以表示为像元中各类别的反射率与其丰度的线性组合,并且假设类别和丰度在基准期和预测期间不发生变化,在一定大小的窗口内,采用带约束条件的最小二乘法求解各类别的平均反射率。混合像元分解模型具体为:
Figure 156001DEST_PATH_IMAGE005
约束条件:
Figure 592799DEST_PATH_IMAGE006
式(1)中,
Figure 780941DEST_PATH_IMAGE007
为第i个混合像元的反射率,
Figure 289283DEST_PATH_IMAGE008
表示第i个混合像元中类别q的平均反射率,m为类别数,
Figure 623312DEST_PATH_IMAGE009
为残余误差;
S14、将上述步骤计算得到的各类别的平均反射率按照类别对应到相应地类的像元上,从而得到粗影像的降尺度数据,依次对所有预测期的粗影像进行混合像元分解。
S2、基准期与预测期的回归模型拟合。基于步骤1中得到的降尺度影像,建立基准期和预测期
Figure 355645DEST_PATH_IMAGE001
间反射率的线性回归模型,并应用于细影像,得到
Figure 407915DEST_PATH_IMAGE001
细影像的反射率拟合结果。
具体地,一定大小的移动窗口M内,根据公式(3)利用建立基准期和预测期的降尺度数据间线性回归模型,利用最小二乘法求得移动窗口内中心像元的回归系数a和b,并应用于细影像,根据公式(4)估测得到
Figure 516685DEST_PATH_IMAGE001
的细影像;
Figure 970800DEST_PATH_IMAGE010
说明书
Figure 375499DEST_PATH_IMAGE011
其中,
Figure 242961DEST_PATH_IMAGE012
Figure 765209DEST_PATH_IMAGE013
分别为移动窗口M内
Figure 198464DEST_PATH_IMAGE002
Figure 413545DEST_PATH_IMAGE001
时期降尺度数据,
Figure 299461DEST_PATH_IMAGE014
为拟合残余误差,
Figure 625401DEST_PATH_IMAGE015
Figure 417557DEST_PATH_IMAGE016
分别为
Figure 865856DEST_PATH_IMAGE002
Figure 380014DEST_PATH_IMAGE001
时期细影像的反射率。
S3、残差的分配和补偿。采用平均聚合的方法将上述步骤2得到的
Figure 368698DEST_PATH_IMAGE001
的细影像的反射率升尺度至低空间分辨率,利用公式(5)计算其与粗影像反射率的差值,即粗分辨率下的残差RD;此项残差包含了上述步骤2中回归模型的拟合残差以及不同传感器间的误差;采用空间插值方法将残差分配至高空间分辨率,空间插值方法可以根据拟合情况选用适当的克里金插值法,将高空间分辨率下的残差rd补偿于
Figure 651912DEST_PATH_IMAGE001
的细影像,即
Figure 333429DEST_PATH_IMAGE001
细影像的反射率加上残差rd作为
Figure 69304DEST_PATH_IMAGE001
残差校正后的融合结果;
Figure 566407DEST_PATH_IMAGE049
式(5)中,
Figure 563182DEST_PATH_IMAGE018
表示低空间分辨率下的残差,
Figure 290966DEST_PATH_IMAGE019
Figure 904350DEST_PATH_IMAGE020
分别为
Figure 110204DEST_PATH_IMAGE002
Figure 227064DEST_PATH_IMAGE001
时期粗像元的反射率,
Figure 952181DEST_PATH_IMAGE021
Figure 397069DEST_PATH_IMAGE022
分别为
Figure 796826DEST_PATH_IMAGE002
Figure 909139DEST_PATH_IMAGE001
时期粗像元中细像元的观测和估测反射率。
S4、拟合结果的空间滤波。在一定窗口内筛选中心像元的光谱相似像元,利用相似性像元的光谱差异性和与中心像元的相对
说明书
距离来加权计算预测期的中心像元值;
具体地,根据公式(6),在w×w的窗口内,综合基准期细影像的各波段相似像元和中心像元的光谱差异D,差异最小的前s个像元被筛选为光谱相似像元;根据公式(7)计算窗口内邻近像元和中心像元的空间距离
Figure 41043DEST_PATH_IMAGE023
,并根据公式(8)确定各相似像元的权重
Figure 629019DEST_PATH_IMAGE024
;各光谱相似像元拟合结果的线性组合视为空间滤波结果。根据公式(9)通过加权窗口内相似性像元及其权重逐波段计算预测时期
Figure 6036DEST_PATH_IMAGE001
的中心像元值,实现对步骤3所得的估测反射率的空间滤波,从而得到
Figure 972855DEST_PATH_IMAGE001
时期空间滤波后的细影像;
Figure DEST_PATH_IMAGE050
Figure 806819DEST_PATH_IMAGE026
Figure 616512DEST_PATH_IMAGE027
Figure 295755DEST_PATH_IMAGE028
其中,
Figure 117080DEST_PATH_IMAGE029
Figure 434796DEST_PATH_IMAGE030
分别表示细影像第b波段的i位置和中心位置像元的反射率,c为用于融合的影像波段数;w为窗口半径,
Figure 341573DEST_PATH_IMAGE031
为高空间分辨率下的残差分配值,
Figure 558927DEST_PATH_IMAGE032
为w×w窗口内中心像元的正向预测结果。
S5、正反向估测结果的加权:如果在
Figure 624972DEST_PATH_IMAGE001
时期之后可以获得第二
说明书
对粗-细影像,那么上述步骤1-4可以反向进行,正、反向估测结果的加权将有利于捕捉基准期与
Figure 941684DEST_PATH_IMAGE001
间发生的地表覆盖类别变化。设置起始日期(
Figure 725969DEST_PATH_IMAGE002
)和结束日期(
Figure 419119DEST_PATH_IMAGE033
)两个基准期,由
Figure 841135DEST_PATH_IMAGE002
的粗-细影像估测
Figure 922224DEST_PATH_IMAGE001
的细影像称为正向估测,由
Figure 803592DEST_PATH_IMAGE033
的粗-细影像估测
Figure 159487DEST_PATH_IMAGE001
的细影像称为反向估测;采用线性加权方法综合正、反向估测结果作为加权时序融合结果;
依据步骤S12得到的
Figure 75490DEST_PATH_IMAGE001
时期细影像的分类图,确定正、反向估测结果中每一地类的权重。具体地,以正向估测过程为例,由
Figure 124218DEST_PATH_IMAGE002
估测
Figure 492882DEST_PATH_IMAGE033
的细影像,根据公式(10),确定估测影像与实际观测影像间的回归系数k和l,并应用于
Figure 354266DEST_PATH_IMAGE001
的正向估测影像,根据公式(11)得到
Figure 983830DEST_PATH_IMAGE001
的似“真实影像”。根据公式(12)和(13)计算
Figure 344404DEST_PATH_IMAGE001
的正向估测影像和似“真实”影像的差距,得到正向估测的相对精度;
Figure 324999DEST_PATH_IMAGE034
Figure 225959DEST_PATH_IMAGE035
Figure 850975DEST_PATH_IMAGE036
Figure 742970DEST_PATH_IMAGE037
其中,
Figure 820647DEST_PATH_IMAGE038
Figure 853194DEST_PATH_IMAGE039
分别表示
Figure 332717DEST_PATH_IMAGE033
的正向估测和实际观测值,
Figure 97411DEST_PATH_IMAGE040
Figure 787018DEST_PATH_IMAGE041
分别表示
Figure 701885DEST_PATH_IMAGE001
的正向估测和似“真实”观测值,
Figure 196101DEST_PATH_IMAGE042
表示正向估测与实际观测值之间的差距,
Figure 803800DEST_PATH_IMAGE043
表示正向估测的相对精度;
相似地,依照上述流程,由
Figure 980703DEST_PATH_IMAGE033
估测
Figure 27156DEST_PATH_IMAGE002
的细影像,确定估测影像
说明书
与实际观测影像间的回归系数,并应用于
Figure 481271DEST_PATH_IMAGE001
的反向估测影像,得到
Figure 384505DEST_PATH_IMAGE001
的似“真实影像”,计算
Figure 189650DEST_PATH_IMAGE001
的反向估测影像和似“真实”影像的差距,可以获得反向估测的相对精度
Figure 337997DEST_PATH_IMAGE044
;根据公式(14)和(15)确定各地类正、反向估测的权重;根据公式(16)获得正、反向估测加权的时序影像融合结果。
Figure 646619DEST_PATH_IMAGE045
Figure 924016DEST_PATH_IMAGE046
Figure 809933DEST_PATH_IMAGE051
S6、结合高空间分辨率影像的连续校正:由于具有较高空间分辨率的Sentinel-2、高分系列等卫星同时提供了较高的时间分辨率,即便受到天气条件等影响,通常也可以获得多幅晴空下的影像,充分利用可获得的多幅细影像将有利于提升融合精度;基于上述步骤S5获得的加权时序影像融合结果,构建地表信息的动态变化模型;利用可获得的多幅已知细影像对动态模型进行修正,使动态模型模拟值与已知细影像值的误差最小,最终得到校正后的时空融合结果;
具体地,采用S-G滤波方法对步骤S5得到的加权时序融合数据进行滤波去噪,并用于构建地表信息的动态模型;采用一种连续校正的CACAO方法(公式(17)),利用多幅已知的细影像修正动态模型,通过最小化模型模拟值与实际观测值之间的误差,获得修正参数,将修正后的各波段反射率作为最终的时序遥感数
说明书
据融合结果;
Figure DEST_PATH_IMAGE052
式(17)中,
Figure 994926DEST_PATH_IMAGE053
表示已知的细影像的各波段反射率,n为可获得的细影像的数量,
Figure DEST_PATH_IMAGE054
表示动态模型模拟的各波段反射率,scale和shift为需要计算的修正参数,可通过最小二乘法获得;对于每一像元,可依据BPpm*(ti+shift)获取任意时间的反射率值。
二、测试数据与处理
本发明实施例以黑河流域中游大满超级站为例,所采用的卫星数据为Sentinel-2多光谱影像(S2-MSI)与准同步的中分辨率成像光谱仪影像所对应的6个波段,具体波段信息如表1所示。
MODIS数据选用2019年6月30日至2019年9月23日间逐天的MCD43A4产品,无云的Sentinel-2影像时间分布如图2所示。
说明书表1 S2-MSIandMODIS数据对应的波段信息
Figure 250065DEST_PATH_IMAGE055
S2-MSILevel-1影像从哨兵科学数据中心获取,作为细影像;
说明书
采用Sen2Cor(version2.8)进行大气校正,投影坐标为UTMWGS-84,同时生成场景分类图,采用SNAP软件将各波段的空间分辨率重采样为20m;MODIS数据从美国航空局网站下载,作为粗影像;采用MODIS投影转换工具(MRT)重投影为UTMWGS-84,空间分辨率重采样为460m。
2019年6月30日作为
Figure DEST_PATH_IMAGE056
基础期,2019年9月23日作为
Figure DEST_PATH_IMAGE057
基础期,采用Kmeans算法将基础期的Sentinel-2影像分为10个类别;将此分类图用于MODIS影像的混合像元分解,窗口设置为31×31;利用基准期的MODIS降尺度影像和Sentinel-2影像及预测期的MODIS降尺度影像分别进行正向和反向的估测,回归拟合中移动窗口内包含115×115个20m分辨率的像元,残差分配与空间滤波中窗口半径设为15,相似像元数为20;依据
Figure 760681DEST_PATH_IMAGE056
时期Sentinel-2影像的分类图,获得每一地类在正、反向估测中的权重,并加权得到2019年6月30日-2019年9月23之间预测期的时间序列影像。此时间序列影像经过S-G滤波后构成地表信息动态模型,利用剩余日期的Sentinel-2影像修正动态模型模拟结果,最终获得预测期连续修正后的高空间分辨率的时间序列遥感数据。
三、结果验证
采用留一影像交叉验证的方法评估融合影像的质量,即在连续校正过程中,预留一期的Sentinel-2影像用于结果验证,其余日期的影像用于修正动态模型,依次对比各日期的融合影像与实
说明书
际观测影像;定性评价采用空间可视化图评估融合影像的空间特征;定量评价选用融合影像与实际观测影像间的相关系数CC、质量指数QI、均方根误差RMSE以及线性拟合系数Slope和截距Intercept,计算公式如下:
Figure DEST_PATH_IMAGE058
(18)
Figure DEST_PATH_IMAGE059
(19)
Figure DEST_PATH_IMAGE060
(20)
CC用于评价估测影像反射率与实际观测影像反射率之间的相关性,负数表示负相关,正数表示正相关。RMSE表示融合误差,值越小表示误差越小;QI用于评价影像的结构相似性,值越接近1表示越相似。
下面以2019年8月4日为例,对比原Fit-FC方法与本发明方法在研究区内对Sentinel-2和MODIS影像时空融合的表现。表2为原Fit-FC方法与本发明方法估测影像的各波段反射率与实际观测反射率的定量评价指标,图3为两种方法融合影像和实际观测影像各波段反射率的散点密度图,图4为两种方法融合影像与实际观测影像的空间可视化图;可以看出,本发明方法融合生成的影像与实际观测影像更接近,消除了原算法产生的模糊现象,呈现出更丰富的纹理特征,极大地提高了融合精度指标,可以提供更准确的反射率。
说明书表2各波段反射率与实际观测反射率的定量评价指标
Figure 900937DEST_PATH_IMAGE061

Claims (7)

1.一种构建高空间分辨率时序遥感数据的方法,其特征在于,具体包含以下步骤:
S1、低空间分辨率影像的混合像元分解:基于高空间分辨率影像的分类图,获取低空间分辨率影像上逐像元的各类地物的丰度,并利用像元分解模型计算每类地物的光谱反射率,获得粗影像的降尺度影像;
S2、基准期与预测期的回归模型拟合:基于步骤S1中得到的降尺度影像,建立基准期和预测期(
Figure RE-DEST_PATH_IMAGE001
)间反射率的线性回归模型,并应用于细影像,得到
Figure RE-823175DEST_PATH_IMAGE001
细影像的反射率拟合结果;
S3、残差的分配和补偿:将上述步骤S2得到的
Figure RE-995530DEST_PATH_IMAGE001
细影像的反射率升尺度至低空间分辨率,计算其与粗影像反射率的差值,即残差,并采用空间插值方法将残差分配至高空间分辨率,并补偿于
Figure RE-884989DEST_PATH_IMAGE001
的细影像;
S4、拟合结果的空间滤波:在一定窗口内筛选中心像元的光谱相似像元,利用相似性像元的光谱差异性和与中心像元的相对距离来加权计算预测期的中心像元值,从而得到空间滤波后的
Figure RE-449569DEST_PATH_IMAGE001
的细影像;
S5、正反向估测结果的加权:设置起始日期(
Figure RE-297439DEST_PATH_IMAGE002
)和结束日期(
Figure RE-816145DEST_PATH_IMAGE001
)两个基准期,将上述步骤S1-S4分正、反两个方向进行,采用线性加权方法综合正、反向估测结果作为加权时序融合结果;
S6、结合高空间分辨率影像的连续校正:基于上述步骤S5获
得的时序影像融合结果,构建地表信息的动态变化模型;利用可获得的多幅已知细影像对动态模型进行修正,使动态模型模拟值与已知细影像值的误差最小,最终得到校正后的时空融合结果。
2.根据权利要求1所述的构建高空间分辨率时序遥感数据的方法,其特征在于,步骤S1中具体包括:
S11:对获得的遥感影像进行大气校正、几何校正、投影转换、重采样、质量控制等预处理;
S12:采用监督分类或非监督分类法对细影像进行分类,利用公式(1)计算粗影像上每个像元内各类地物的丰度;
Figure RE-509295DEST_PATH_IMAGE004
式(1)中,
Figure RE-DEST_PATH_IMAGE005
表示i位置处混合像元内类别q的丰度,Q为混合像元中类别q的细像元数,S为混合像元中所有类别的细像元总数;
S13:混合像元的反射率值可以表示为像元中各类别的反射率与其丰度的线性组合,混合像元分解模型具体为:
Figure RE-DEST_PATH_IMAGE007
约束条件:
Figure RE-633108DEST_PATH_IMAGE008
式(2)中,
Figure RE-DEST_PATH_IMAGE009
为第i个混合像元的反射率,
Figure RE-979776DEST_PATH_IMAGE010
表示第i个混合像元中类别q的平均反射率,m为类别数,
Figure RE-DEST_PATH_IMAGE011
为残余误差。
S14:将上述步骤计算得到的各类别的平均反射率按照类别对
应到相应地类的像元上,从而得到粗影像的降尺度数据。
3.根据权利要求1所述的构建高空间分辨率时序遥感数据的方法,其特征在于,步骤S2中具体包括:
根据公式(3)利用基准期和预测期的降尺度数据建立线性回归模型,需要在一定大小的移动窗口M内进行,利用最小二乘法求得移动窗口内中心像元的回归系数a和b,并应用于细影像,根据公式(4)估测得到
Figure RE-221664DEST_PATH_IMAGE001
的细影像;
Figure RE-DEST_PATH_IMAGE013
Figure RE-DEST_PATH_IMAGE015
其中,
Figure RE-46400DEST_PATH_IMAGE016
Figure RE-DEST_PATH_IMAGE017
分别为移动窗口M内
Figure RE-555879DEST_PATH_IMAGE002
Figure RE-479973DEST_PATH_IMAGE001
时期降尺度数据的反射率,
Figure RE-973271DEST_PATH_IMAGE018
为拟合残余误差,
Figure RE-DEST_PATH_IMAGE019
Figure RE-742644DEST_PATH_IMAGE020
分别表示
Figure RE-401902DEST_PATH_IMAGE002
Figure RE-231318DEST_PATH_IMAGE001
时期细影像的反射率。
4.根据权利要求1所述的构建高空间分辨率时序遥感数据的方法,其特征在于,步骤S3中具体包括:
采用平均聚合的方法将
Figure RE-211912DEST_PATH_IMAGE001
的细影像升尺度至低空间分辨率,利用公式(5)计算残差RD;将残差分配至高空间分辨率采用的空间插值方法可以选用适当的克里金插值法,并将高空间分辨率的残差rd补偿于与上述步骤S2得到的
Figure RE-316134DEST_PATH_IMAGE001
细影像的反射率;
Figure RE-675571DEST_PATH_IMAGE022
权利要求书
式(5)中,
Figure RE-DEST_PATH_IMAGE023
表示低空间分辨率下的残差,
Figure RE-800522DEST_PATH_IMAGE024
Figure RE-DEST_PATH_IMAGE025
分别为
Figure RE-471675DEST_PATH_IMAGE002
Figure RE-379588DEST_PATH_IMAGE001
时期粗像元的反射率,
Figure RE-485210DEST_PATH_IMAGE026
Figure RE-DEST_PATH_IMAGE027
分别为
Figure RE-781062DEST_PATH_IMAGE002
Figure RE-346035DEST_PATH_IMAGE001
时期粗像元中细像元的观测和估测反射率。
5.根据权利要求1所述的构建高空间分辨率时序遥感数据的方法,其特征在于,步骤S4中具体包括:
根据公式(6),综合基准期细影像的各波段相似像元和中心像元的光谱差异D,在一定窗口内,差异最小的前s个像元被筛选为相似像元;根据公式(7)计算邻近像元和中心像元的空间距离
Figure RE-792060DEST_PATH_IMAGE028
,根据公式(8)确定各相似像元的权重
Figure RE-DEST_PATH_IMAGE029
;根据公式(9)通过加权窗口内相似性像元逐波段计算预测时期
Figure RE-985144DEST_PATH_IMAGE001
的中心像元值,从而得到
Figure RE-717477DEST_PATH_IMAGE001
时期空间滤波后的细影像;
Figure RE-DEST_PATH_IMAGE031
Figure RE-DEST_PATH_IMAGE033
Figure RE-DEST_PATH_IMAGE035
Figure RE-65019DEST_PATH_IMAGE036
其中,
Figure RE-DEST_PATH_IMAGE037
Figure RE-642631DEST_PATH_IMAGE038
分别表示细影像第b波段的i位置和中心位置像元的反射率,c为用于融合的影像波段数;w为窗口半
径,
Figure RE-DEST_PATH_IMAGE039
为高空间分辨率下的残差分配值,
Figure RE-955801DEST_PATH_IMAGE040
为w×w窗口内中心像元的正向预测结果。
6.根据权利要求1所述的构建高空间分辨率时序遥感数据的方法,其特征在于,步骤S5中具体包括:
依据
Figure RE-734401DEST_PATH_IMAGE002
细影像的分类图,确定正反向估测结果中每一地类的权重;具体地,以正向估测过程为例,根据公式(10),由
Figure RE-398600DEST_PATH_IMAGE002
估测
Figure RE-DEST_PATH_IMAGE041
的细影像,确定估测影像与实际观测影像间的回归系数k和l,并应用于
Figure RE-920849DEST_PATH_IMAGE001
的正向估测影像,根据公式(11)得到
Figure RE-855569DEST_PATH_IMAGE001
的似“真实影像”;根据公式(12)和(13)计算
Figure RE-805070DEST_PATH_IMAGE001
的正向估测影像和似“真实”影像的差距,得到正向估测的相对精度;
Figure RE-DEST_PATH_IMAGE043
Figure RE-DEST_PATH_IMAGE045
Figure RE-DEST_PATH_IMAGE047
Figure RE-DEST_PATH_IMAGE049
其中,
Figure RE-222145DEST_PATH_IMAGE050
Figure RE-DEST_PATH_IMAGE051
分别表示
Figure RE-407139DEST_PATH_IMAGE041
的正向估测和真实观测值,
Figure RE-835846DEST_PATH_IMAGE052
Figure RE-DEST_PATH_IMAGE053
分别表示
Figure RE-48259DEST_PATH_IMAGE001
的正向估测和似“真实”观测值,
Figure RE-952630DEST_PATH_IMAGE054
表示正向估测与真实观测值之间的差距,
Figure RE-DEST_PATH_IMAGE055
表示正向估测的相对精度;
相似地,依照上述流程,可以获得反向估测的相对精度
Figure RE-410156DEST_PATH_IMAGE056
;根据公式(14)和(15)确定各地类正反向估测的权重;根据公式(16)获得正反向估测加权的时序影像融合结果;
Figure RE-693370DEST_PATH_IMAGE058
Figure RE-250254DEST_PATH_IMAGE060
Figure RE-DEST_PATH_IMAGE061
7.根据权利要求1所述的构建高空间分辨率时序遥感数据的方法,其特征在于,步骤S6中具体包括:
采用S-G滤波方法对时序融合数据进行滤波用于构建地表信息动态模型;采用一种连续校正的CACAO方法(公式(17)),利用多幅已知的细影像修正动态模型,将修正后的各波段反射率作为最终的时序数据融合结果;
Figure RE-DEST_PATH_IMAGE063
式(17)中,BPture表示已知的细影像的各波段反射率,n为可获得的细影像的数量,BPpm表示动态模型模拟的各波段反射率,scale和shift为需要计算的修正参数。
CN202210748388.8A 2022-06-29 2022-06-29 一种构建高空间分辨率时序遥感数据的方法 Active CN115063332B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210748388.8A CN115063332B (zh) 2022-06-29 2022-06-29 一种构建高空间分辨率时序遥感数据的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210748388.8A CN115063332B (zh) 2022-06-29 2022-06-29 一种构建高空间分辨率时序遥感数据的方法

Publications (2)

Publication Number Publication Date
CN115063332A true CN115063332A (zh) 2022-09-16
CN115063332B CN115063332B (zh) 2024-04-30

Family

ID=83203656

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210748388.8A Active CN115063332B (zh) 2022-06-29 2022-06-29 一种构建高空间分辨率时序遥感数据的方法

Country Status (1)

Country Link
CN (1) CN115063332B (zh)

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103177431A (zh) * 2012-12-26 2013-06-26 中国科学院遥感与数字地球研究所 一种多源遥感数据时空融合方法
CN104715467A (zh) * 2015-03-06 2015-06-17 中国科学院遥感与数字地球研究所 一种改进型多源遥感数据时空融合方法
CN105046648A (zh) * 2015-06-25 2015-11-11 北京师范大学 一种构建高时空遥感数据的方法
US20200082151A1 (en) * 2018-09-06 2020-03-12 National Central University Method of Top-of-Atmosphere Reflectance-Based Spatiotemporal Image Fusion Using Aerosol Optical Depth
CN111832518A (zh) * 2020-07-22 2020-10-27 桂林电子科技大学 基于时空融合的tsa遥感影像土地利用方法
US20210118097A1 (en) * 2018-02-09 2021-04-22 The Board Of Trustees Of The University Of Illinois A system and method to fuse multiple sources of optical data to generate a high-resolution, frequent and cloud-/gap-free surface reflectance product
CN113408929A (zh) * 2021-07-01 2021-09-17 福州大学 基于空间几何原理的四维遥感生态指数构建方法
CN114202699A (zh) * 2021-12-15 2022-03-18 南京大学 一种地表水体月度变化动态监测方法
CN114565843A (zh) * 2022-02-22 2022-05-31 中国电子科技集团公司第五十四研究所 一种时间序列遥感图像融合方法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103177431A (zh) * 2012-12-26 2013-06-26 中国科学院遥感与数字地球研究所 一种多源遥感数据时空融合方法
CN104715467A (zh) * 2015-03-06 2015-06-17 中国科学院遥感与数字地球研究所 一种改进型多源遥感数据时空融合方法
CN105046648A (zh) * 2015-06-25 2015-11-11 北京师范大学 一种构建高时空遥感数据的方法
US20210118097A1 (en) * 2018-02-09 2021-04-22 The Board Of Trustees Of The University Of Illinois A system and method to fuse multiple sources of optical data to generate a high-resolution, frequent and cloud-/gap-free surface reflectance product
US20200082151A1 (en) * 2018-09-06 2020-03-12 National Central University Method of Top-of-Atmosphere Reflectance-Based Spatiotemporal Image Fusion Using Aerosol Optical Depth
CN111832518A (zh) * 2020-07-22 2020-10-27 桂林电子科技大学 基于时空融合的tsa遥感影像土地利用方法
CN113408929A (zh) * 2021-07-01 2021-09-17 福州大学 基于空间几何原理的四维遥感生态指数构建方法
CN114202699A (zh) * 2021-12-15 2022-03-18 南京大学 一种地表水体月度变化动态监测方法
CN114565843A (zh) * 2022-02-22 2022-05-31 中国电子科技集团公司第五十四研究所 一种时间序列遥感图像融合方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
袁周米琪;张锦水;: "面向地表特征变化区域的时空遥感数据融合方法研究", 北京师范大学学报(自然科学版), no. 06, 15 December 2017 (2017-12-15) *
谢登峰;张锦水;孙佩军;潘耀忠;云雅;袁周米琪;: "结合像元分解和STARFM模型的遥感数据融合", 遥感学报, no. 01, 25 January 2016 (2016-01-25) *

Also Published As

Publication number Publication date
CN115063332B (zh) 2024-04-30

Similar Documents

Publication Publication Date Title
Zhao et al. A robust adaptive spatial and temporal image fusion model for complex land surface changes
CN107103584B (zh) 一种基于时空加权的生产高时空分辨率ndvi的方法
Song et al. Spatiotemporal satellite image fusion through one-pair image learning
CN110070518B (zh) 一种基于双路径支持下的高光谱图像超分辨率制图方法
US8879871B2 (en) Method for processing images using automatic georeferencing of images derived from a pair of images captured in the same focal plane
CN109685108B (zh) 一种生成高时空分辨率ndvi长时间序列的方法
CN103617629B (zh) 基于modis遥感影像的高分辨率遥感影像植被指数时间序列校正方法
CN112991288A (zh) 基于丰度图像锐化重构的高光谱遥感图像融合方法
CN108427964B (zh) 一种遥感图像与地球化学的融合方法及系统
CN110503137B (zh) 基于交叉融合的遥感影像时空融合基础图像对的确定方法
Gómez-Chova et al. Gridding artifacts on medium-resolution satellite image time series: MERIS case study
Wu et al. Spatiotemporal fusion with only two remote sensing images as input
Cenci et al. Describing the quality assessment workflow designed for DEM products distributed via the Copernicus Programme. Case study: The absolute vertical accuracy of the Copernicus DEM dataset in Spain
CN104637027A (zh) 顾及非局部特性与时空变化的遥感数据时空定量融合方法
Sihvonen et al. Spectral profile partial least-squares (SP-PLS): Local multivariate pansharpening on spectral profiles
CN115630308A (zh) 一种降尺度和融合联合的地表温度时空分辨率增强方法
CN111666896A (zh) 一种基于线性融合模型的遥感影像时空融合方法
CN113744134B (zh) 基于光谱解混卷积神经网络的高光谱图像超分辨方法
CN109359264B (zh) 一种基于modis的叶绿素产品降尺度方法及装置
CN115272144A (zh) 面向高光谱图像和多光谱图像的时空谱融合方法
CN109840539B (zh) 一种基于地块图斑的遥感时空数据融合方法
CN115063332B (zh) 一种构建高空间分辨率时序遥感数据的方法
CN109472237B (zh) 一种可见光遥感卫星影像的大气订正方法和系统
CN114639014B (zh) 一种基于高分辨率遥感影像的ndvi归一化方法
CN110991567A (zh) 一种由点到面的遥感瞬时地表温度数据检测方法

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