CN113554011A - 遥感影像变化检测方法、装置、设备及存储介质 - Google Patents

遥感影像变化检测方法、装置、设备及存储介质 Download PDF

Info

Publication number
CN113554011A
CN113554011A CN202111104059.1A CN202111104059A CN113554011A CN 113554011 A CN113554011 A CN 113554011A CN 202111104059 A CN202111104059 A CN 202111104059A CN 113554011 A CN113554011 A CN 113554011A
Authority
CN
China
Prior art keywords
image
change
surface reflectivity
spectral
remote sensing
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
CN202111104059.1A
Other languages
English (en)
Other versions
CN113554011B (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.)
Henan Aerospace Hongtu Information Technology Co ltd
Aerospace Hongtu Information Technology Co Ltd
Original Assignee
Aerospace Hongtu Information Technology Co Ltd
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 Aerospace Hongtu Information Technology Co Ltd filed Critical Aerospace Hongtu Information Technology Co Ltd
Priority to CN202111104059.1A priority Critical patent/CN113554011B/zh
Publication of CN113554011A publication Critical patent/CN113554011A/zh
Application granted granted Critical
Publication of CN113554011B publication Critical patent/CN113554011B/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
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods

Abstract

本申请提供了一种遥感影像变化检测方法、装置、设备及存储介质，涉及遥感变化检测工程化技术领域，包括：对两个不同时序的原始影像分别进行影像处理，生成第一地表反射率影像和第二地表反射率影像；对第一地表反射率影像和第二地表反射率影像进行光谱特征分析，生成光谱差异性影像；对第一地表反射率影像和第二地表反射率影像进行影像线性变换，生成变化差异影像；基于光谱差异性影像和变换差异影像，生成基于光谱特征的变化影像；利用动态变化阈值从变化影像中划分出变化区域。本申请能够提高遥感影像变化的检测精度和检测速度。

Description

遥感影像变化检测方法、装置、设备及存储介质
技术领域
本申请涉及遥感变化检测工程化领域，尤其是涉及遥感影像变化检测方法、装置、设备及存储介质。
背景技术
实时精确地获取地表变化信息对于国土监测、环境保护、减灾防灾等领域有着重要意义。遥感对地观测技术具有大范围、长时间和周期性监测的优势，利用多时遥感影像来获取地表地物变化信息的变化检测是遥感技术应用最广泛的研究领域之一。
遥感变化检测是指基于覆盖同一地表区域的多时相遥感影像数据，采用图像处理和模式识别等方法提取变化信息，并分析和确定地表变化的特征与过程。
目前，遥感变化检测方法大体可分为影像代数法、影像变换法和分类检测法。影像代数法通过计算多时相遥感影像对应波段间的代数特征来衡量变化情况，主要是对两时影像进行差值或比值运算，以获取变化信息，但是该方法易受成像环境、遥感传感器参数差异等因素影响，在实际工程化应用中稳定性差、噪声大。影像变换法是利用统计模型，对多时相影像数据做数学变换，提取或突出变化地物特征信息，但是该方法主要基于未做大气校正的DN或辐亮度数据，无法有效减少大气对变化检测的影响，并且没有充分利用遥感影像特有的光谱、地理属性等特征信息，所以其变化检测的精度相对较低。分类检测法主要利用各种分类器通过监督分类技术进行变化检测，但是该方法需要大量的人工解译样本训练，所以自动化程度较低且算法相对复杂，不利于工程化应用和快速检测。
发明内容
有鉴于此，本申请提供一种遥感影像变化检测方法、装置、设备及存储介质，以解决现有遥感变化检测方法存在的自动化程度低且精度低的技术问题。
一方面，本申请实施例提供了一种遥感影像变化检测方法，包括：
对两个不同时序的原始影像分别进行影像处理，生成第一地表反射率影像和第二地表反射率影像；
对第一地表反射率影像和第二地表反射率影像进行光谱特征分析，生成光谱差异性影像；
对第一地表反射率影像和第二地表反射率影像进行影像线性变换，生成变化差异影像；
基于光谱差异性影像和变换差异影像，生成基于光谱特征的变化影像；
利用动态变化阈值从变化影像中划分出变化区域。
另一方面，本申请实施例提供了一种遥感影像变化检测装置，包括：
地表反射率影像生成单元，用于对两个不同时序的原始影像分别进行影像处理，生成第一地表反射率影像和第二地表反射率影像；
光谱差异性影像生成单元，用于对第一地表反射率影像和第二地表反射率影像进行光谱特征分析，生成光谱差异性影像；
变化差异影像生成单元，用于对第一地表反射率影像和第二地表反射率影像进行影像线性变换，生成变化差异影像；
变化影像生成单元，用于基于光谱差异性影像和变换差异影像，生成基于光谱特征的变化影像；
变化区域监测单元，用于利用动态变化阈值从变化影像中划分出变化区域。
另一方面，本申请实施例提供了一种电子设备，包括：存储器、处理器和存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现本申请实施例的遥感影像变化检测方法。
另一方面，本申请实施例提供了一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器运行时实现本申请实施例的遥感影像变化检测方法。
与现有技术相比，本申请的优势在于：
本申请实施例能够充分利用遥感影像的光谱特征，并且可以有效减少大气条件、光照、传感器差异、影像色差等因素对变化检测的影响，对水体、植被、建筑物等主要遥感地物变化敏感，具备适用性强、准确度高、自动化程度高等优势，能够为土地调查、生态变化、灾害检测和评估等领域方面的提供基础技术支撑。
附图说明
为了更清楚地说明本申请具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本申请的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本申请实施例提供的遥感影像变化检测方法的流程图 ；
图2为本申请实施例提供的生成两个不同时序的地表反射率影像的流程图；
图3（a）为本申请仿真实例的第一时序的原始影像；
图3（b）为本申请仿真实例的第二时序的原始影像；
图4为本申请仿真实例的光谱差异性影像；
图5为本申请仿真实例的变化差异影像；
图6为本申请仿真实例的变化影像；
图7为本申请仿真实例的动态阈值分割后的影像；
图8为本申请仿真实例的优化后的变化区域影像；
图9为本申请实施例提供的遥感影像变化检测装置的功能结构示意图；
图10为本申请实施例提供的电子设备的硬件结构示意图。
具体实施方式
为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本申请一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本申请实施例的组件可以以各种不同的配置来布置和设计。
首先对本申请实施例的设计思想进行简单介绍。
目前，遥感变化检测方法大体可分为影像代数法、影像变换法和分类检测法。其中，影像代数法易受成像环境、遥感传感器参数差异等因素影响，在实际工程化应用中稳定性差、噪声大。影像变换法的变化检测精度相对较低。分类检测法的自动化程度较低且算法相对复杂，不利于工程化应用和快速检测。
考虑到以上方法在实际应用中稳定性、准确性、简便性上有一定的限制，本申请实施例利用影像光谱特征提出了一种遥感影像变化检测方法，该方法遍历两期影像的每个像元，计算光谱差异性值，则可以构建用于表征两期影像差异性光谱特征的光谱差异性影像D1，影像值越大，变化的程度越大。利用典型关联分析算法对两期影像进行组合，并进行线性变换，有效减少光照、传感器参数、影像色差等带来的误差，同时可以最大程度的增强变化信息，进而构建差异影像D2，影像值越大，变化的程度越大。基于光谱差异性影像和影像变换的差异影像，构建顾及光谱特征的变化影像D。计算变化影像值域92.5%百分位所在的像元值，得到的分割结果。动态阈值分割结果有很多噪声，需要利用过滤技术对分割结果做进一步优化整合，得到最终的变化检测成果。
本申请实施例提出的遥感影像变化监测方法，能够充分利用遥感影像光谱特征，具有准确性高、自动化程度高以及便于工程化应用的优点。
在介绍了本申请实施例的设计思想之后，下面对本申请实施例提供的技术方案进行说明。
实施例一：
如图1所示，本申请实施例提出了一种遥感影像变化检测方法，包括：
步骤101：对两个不同时序的原始影像分别进行影像处理，生成地表反射率影像X和地表反射率影像Y；
在步骤101中，原始影像中包含影像的质量控制信息、卫星几何信息、传感器几何信息、获取时间等元数据。
如图2所示，该步骤包括：
步骤201：利用第一时序的原始影像包含的元数据对第一时序的原始影像进行大气校正，生成第一初始地表反射率影像；利用第二时序的原始影像包含的元数据对第二时序的原始影像进行大气校正，生成第二初始地表反射率影像；
在步骤201中，通过大气校正可以减少大气条件引入的误差；原始影像为原始多光谱或高光谱影像。
步骤202：利用第一时序的原始影像包含的质量控制信息对第一初始地表反射率影像进行有效区域掩膜处理，得到第一影像；利用第二时序的原始影像包含的质量控制信息对第二初始地表反射率影像进行有效区域掩膜处理，得到第二影像；
在步骤202中，通过有效区域掩膜处理可以去除影像的无效像元。
步骤203：对第一影像和第二影像进行影像配准；
在步骤203中，影像配准是为了保证两期影像在空间分布上的一致性。
步骤204：对配准后的第一影像和第二影像分别进行重采样和裁切，最终生成第一时序的地表反射率影像X和第二时序的地表反射率影像Y；
在步骤204中，影像X和影像Y值域为0到1，且影像X和影像Y都在同一区域、具备相同的大小、分辨率和波段数。
步骤102：对影像X和影像Y进行光谱特征分析，构建用于表征两期影像差异性光谱特征的光谱差异性影像D1；
在步骤102中，光谱特征分析主要为了充分利用反射率影像的光谱信息，分析两期影像的光谱差异性，用于增强变化信息。通过光谱角和光谱距离构建光谱差异性影像实现，具体的过程如下：
影像X和影像Y都有n个波段，大小都为H×W，H为影像的高度，W为影像的宽度；X和Y每个像元都有n个光谱值，即都有1维的n个元素的光谱向量。对每一个像元计算光谱差异性值d(i,j)，计算公式如下：
Figure P_210918102809673_673976001
其中，d(i,j)为(i,j)位置处像元的光谱差异性值，
Figure P_210918102809752_752684001
为影像X中位置(i,j)处像元光谱向量，
Figure P_210918102809799_799543002
为影像Y中位置(i,j)处像元光谱向量，
Figure P_210918102809830_830782003
为
Figure P_210918102809862_862008004
和
Figure P_210918102809894_894739005
的光谱角，
Figure P_210918102809941_941290006
为
Figure P_210918102809972_972874007
和
Figure P_210918102810004_004134008
的光谱距离，
Figure P_210918102810051_051225009
,
Figure P_210918102810084_084549010
，N为波段数，xn为
Figure P_210918102810131_131543011
向量第n波段光谱值，yn为
Figure P_210918102810162_162875012
向量第n波段光谱值。
计算各个像元光谱差异性值，则可以构建用于表征两期影像差异性光谱特征的光谱差异性影像D1，D1的值域为0到1，若像元光谱差异性值约接近1，则该像元发生变化的概率越大，越接近于0，则该像元发生变化的概率越小。
步骤103：对影像X和影像Y进行影像线性变换，构建差异影像D2；
在步骤103中，通过影像线性变换，计算修正系数，则可以有效减少光照、传感器参数、影像色差等带来的误差，同时可以最大程度的增强变化信息。
利用典型关联分析算法(Canonical Correlation Analysis，CCA)对两期影像进行组合，并进行线性变换，获取两期影像的变换系数，进而构建差异影像D2；具体实现如下：
将影像X转换为大小为N*M，均值为0的矩阵FX，其中M=H×W，矩阵FX的每一列为第一地表反射率影像的像元光谱向量；
将影像Y转换为大小为N*M，均值为0的矩阵FY，矩阵FY的每一列为第二地表反射率影像的像元光谱向量；
计算FX和FY的协方差
Figure P_210918102810194_194008001
、FY和FX的协方差
Figure P_210918102810225_225325002
、FX的方差
Figure P_210918102810272_272595003
以及FY的方差
Figure P_210918102810303_303395004
；
根据相关性系数公式，构建优化目标函数：
Figure P_210918102810334_334752001
其中，相关性系数向量
Figure P_210918102810381_381927001
，所对应的数值为两幅影像N个波段间对应的相关性系数，α为FX的线性变换系数矩阵，大小为N*N；b为FY的线性变换系数矩阵，大小为N*N；
由于目标函数的分子分母增大相同的倍数，优化目标结果不变，上述目标函数可转换为：
Figure P_210918102810412_412779001
构造Lagrange等式后求导，得到下列公式：
Figure P_210918102810459_459694001
对
Figure P_210918102810492_492886001
进行特征分解，解算出相关性系数向量ρ和矩阵α；
对ρ的各个元素从小到大进行排序，按照该顺序调整α的列向量的顺序；
然后将重排后的ρ和α，代入
Figure P_210918102810524_524121001
中，解算出矩阵b；
计算矩阵X'和矩阵Y'：
Figure P_210918102810570_570991001
利用矩阵X'和矩阵Y'，计算FX和FY之间的变化差异矩阵
Figure P_210918102810602_602362001
：
Figure P_210918102810649_649187001
其中，
Figure P_210918102810681_681291001
为X'第n波段的像元值，大小为1*M;
Figure P_210918102810713_713061002
为Y'第n波段的像元值，大小为1*M；
Figure P_210918102810743_743909003
的大小为1*M；
将
Figure P_210918102810790_790779001
的各个元素按照H*W的影像格式进行重新排列，生成差异影像D2。
步骤104：基于光谱差异性影像D1和变换差异影像D2，构建顾及光谱特征的变化影像；
对于光谱差异影像D1位置(i,j)处像元值
Figure P_210918102810822_822575001
，以及变化差异影像D2位置(i,j)处像元值
Figure P_210918102810853_853680002
，计算两者之间的乘积，作为变化影像D位置(i,j)处像元值
Figure P_210918102810886_886415003
：
Figure P_210918102810933_933337001
由此生成变化影像D。
步骤105：利用动态变化阈值从变化影像中划分出变化区域；
变化阈值是基于影像本身的动态变化阈值进行变化检测，即将变化影像D中的大于某一个阈值的像元，判定为变化像元。
通过大量统计，变化检测阈值取值为变化影像值域92.5%百分位所在的像元值。
步骤106：对变化区域进行优化。
利用过滤处理（Sieve）解决变化检测结果数据出现的孤岛或空洞问题。过滤处理使用斑点分组方法来消除被隔离或填补断裂的目标像元。
仿真实例
选取北京地区的两期Landsat8影像，时间分别为2014年9月4日、2019年9月2日，如图3（a）和图3（b）所示。每期影像都包含质量控制文件QA影像以及卫星、传感器的元数据信息。通过6s对Landsat8数据进行大气校正，并对大气校正后的数据进行掩膜、配准、重采样、裁切等处理，获取了两景分辨率为30米、大小为5000*5000、波段数为7的地表反射率影像X和Y。
遍历每个像元计算光谱差异性值，则可以构建用于表征两期影像差异性光谱特征的光谱差异性影像D1，如图4所示，影像值越大，变化的程度越大。利用典型关联分析算法对两期影像进行组合，并进行线性变换，有效减少光照、传感器参数、影像色差等带来的误差，同时可以最大程度的增强变化信息，进而构建差异影像D2，如图5所示，影像值越大，变化的程度越大。基于光谱差异性影像和影像变换的差异影像，构建顾及光谱特征的变化影像D，如图6所示。计算变化影像值域92.5%百分位所在的像元值，得到的分割结果，如图7所示。动态阈值分割结果有很多噪声，需要利用过滤技术对分割结果做进一步优化整合，得到最终的变化检测成果，如图8所示。
实施例二：
基于上述实施例，本申请实施例提供了一种遥感影像变化检测装置，参阅图9所示，本申请实施例提供的遥感影像变化检测装置300至少包括：
地表反射率影像生成单元301，用于对两个不同时序的原始影像分别进行影像处理，生成第一地表反射率影像和第二地表反射率影像；
光谱差异性影像生成单元302，用于对第一地表反射率影像和第二地表反射率影像进行光谱特征分析，生成光谱差异性影像；
变化差异影像生成单元303，用于对第一地表反射率影像和第二地表反射率影像进行影像线性变换，生成变化差异影像；
变化影像生成单元304，用于基于光谱差异性影像和变换差异影像，生成基于光谱特征的变化影像；
变化区域检测单元305，用于利用动态变化阈值从变化影像中划分出变化区域；
优化单元306，用于消除变化区域的孤岛或填补变化区域的空洞。
在一种可能的实施方式中，地表反射率影像生成单元301具体用于：
利用第一时序的原始影像包含的元数据对第一时序的原始影像进行大气校正，生成第一原始地表反射率影像；利用第二时序的原始数据包含的元数据对第二时序的原始影像进行大气校正，生成第二原始地表反射率影像；
利用第一时序的原始影像包含的质量控制信息对第一原始地表反射率影像进行有效区域掩膜处理，得到第一影像;根据第二时序的原始数据包含的质量控制信息对第二原始地表反射率影像进行有效区域掩膜处理，得到第二影像;
对第一影像和第二影像进行影像配准；
对配准后的第一影像进行重采样和裁切，生成第一地表反射率影像X；对配准后的第二影像进行重采样和裁切，生成第二地表反射率影像Y。
在另一种可能的实施方式中，光谱差异性影像生成单元302具体用于：
获取第一地表反射率影像X位置(i,j)处像元光谱向量
Figure P_210918102810964_964542001
，和第二地表反射率影像Y位置(i,j)处像元光谱向量
Figure P_210918102810996_996040002
，N为波段数，xn为
Figure P_210918102811042_042681003
的第n波段光谱值，yn为
Figure P_210918102811075_075336004
第n波段光谱值；第一地表反射率影像和第二地表反射率影像的大小均为H*W，1≤i≤H，1≤j≤W；
计算光谱向量
Figure P_210918102811107_107111001
和光谱向量
Figure P_210918102811138_138354002
之间的光谱角
Figure P_210918102811169_169606003
：
Figure P_210918102811200_200853001
计算光谱向量
Figure P_210918102811249_249214001
和光谱向量
Figure P_210918102811298_298124002
之间的距离
Figure P_210918102811344_344941003
：
Figure P_210918102811376_376179001
计算第一地表反射率影像和第二地表反射率影像的(i,j)位置处像元的光谱差异性值D1(i,j)：
Figure P_210918102811423_423157001
基于各个像元光谱差异性值，生成光谱差异性影像D1。
在另一种可能的实施方式中，变化差异影像生成单元303具体用于：
将第一地表反射率影像转换为大小为N*M，均值为0的矩阵FX，其中M=H×W，矩阵FX的每一列为第一地表反射率影像的像元光谱向量；
将第二地表反射率影像转换为大小为N*M，均值为0的矩阵FY，矩阵FY的每一列为第二地表反射率影像的像元光谱向量；
计算FX和FY的协方差
Figure P_210918102811470_470844001
、FY和FX的协方差
Figure P_210918102811518_518261002
、FX的方差
Figure P_210918102811549_549496003
以及FY的方差
Figure P_210918102811596_596521004
；SXY、SYX、SXX和SYY均为N*N的矩阵；
根据相关性系数公式，构建优化目标函数：
Figure P_210918102811642_642864001
其中，相关性系数向量
Figure P_210918102811691_691601001
，所对应的数值为两幅影像N个波段间对应的相关性系数，α为FX的线性变换系数矩阵，大小为N*N；b为FY的线性变换系数矩阵，大小为N*N；
上述目标函数可转换为：
Figure P_210918102811738_738511001
构造Lagrange等式后求导，得到下列公式：
Figure P_210918102811785_785353001
对
Figure P_210918102811816_816618001
进行特征分解，解算出相关性系数向量ρ和矩阵α；
对ρ的各个元素从小到大进行排序，按照该顺序调整α的列向量的顺序；
然后将重排后的ρ和α，代入
Figure P_210918102811863_863504001
中，解算出矩阵b；
计算矩阵X'和矩阵Y'：
Figure P_210918102811894_894786001
利用矩阵X'和矩阵Y'，计算FX和FY之间的变化差异矩阵
Figure P_210918102811941_941610001
：
Figure P_210918102811972_972956001
其中，
Figure P_210918102812019_019701001
为X'第n波段的像元值，大小为1*M;
Figure P_210918102812051_051020002
为Y'第n波段的像元值，大小为1*M；
Figure P_210918102812082_082614003
的大小为1*M；
将
Figure P_210918102812113_113485001
的各个元素按照H*W的影像格式进行重新排列，生成差异影像D2。
在另一种可能的实施方式中，变化影像生成单元304具体用于：
对于光谱差异影像D1位置(i,j)处像元值
Figure P_210918102812144_144791001
，以及变化差异影像D2位置(i,j)处像元值
Figure P_210918102812191_191653002
，计算两者之间的乘积，作为变化影像D位置(i,j)处像元值
Figure P_210918102812222_222894003
：
Figure P_210918102812254_254115001
由此生成变化影像D。
在另一种可能的实施方式中，变化区域检测单元305具体用于：
对于变化影像的各个像元，如果像元的像元值大于变化阈值，则判断该像元为变化像元，由此划分出变化区域。
需要说明的是，由于本申请实施例提供的遥感影像变化检测装置300解决技术问题的原理与本申请实施例提供的遥感影像变化检测方法相似，因此，本申请实施例提供的遥感影像变化检测装置300的实施可以参见本申请实施例提供的遥感影像变化检测方法的实施，重复之处不再赘述。
实施例三：
基于上述实施例，本申请实施例还提供了一种电子设备，参阅图10所示，本申请实施例提供的电子设备400至少包括：处理器401、存储器402和存储在存储器402上并可在处理器401上运行的计算机程序，处理器401执行计算机程序时实现本申请实施例提供的遥感影像变化检测方法。
本申请实施例提供的电子设备400还可以包括连接不同组件（包括处理器401和存储器402）的总线403。其中，总线403表示几类总线结构中的一种或多种，包括存储器总线、外围总线、局域总线等。
存储器402可以包括易失性存储器形式的可读介质，例如随机存储器（RandomAccess Memory，RAM）4021和/或高速缓存存储器4022，还可以进一步包括只读存储器（ReadOnly Memory，ROM）4023。
存储器402还可以包括具有一组（至少一个）程序模块4025的程序工具4024，程序模块4025包括但不限于：操作子系统、一个或者多个应用程序、其它程序模块以及程序数据，这些示例中的每一个或某种组合中可能包括网络环境的实现。
电子设备400也可以与一个或多个外部设备404（例如键盘、遥控器等）通信，还可以与一个或者多个使得用户能与电子设备400交互的设备通信（例如手机、电脑等），和/或，与使得电子设备400与一个或多个其它电子设备400进行通信的任何设备（例如路由器、调制解调器等）通信。这种通信可以通过输入/输出（Input /Output，I/O）接口403进行。并且，电子设备400还可以通过网络适配器406与一个或者多个网络（例如局域网（Local AreaNetwork，LAN），广域网（Wide Area Network，WAN）和/或公共网络，例如因特网）通信。如图10所示，网络适配器406通过总线403与电子设备400的其它模块通信。应当理解，尽管图10中未示出，可以结合电子设备400使用其它硬件和/或软件模块，包括但不限于：微代码、设备驱动器、冗余处理器、外部磁盘驱动阵列、磁盘阵列（Redundant Arrays of IndependentDisks，RAID）子系统、磁带驱动器以及数据备份存储子系统等。
需要说明的是，图10所示的电子设备400仅仅是一个示例，不应对本申请实施例的功能和使用范围带来任何限制。
实施例四：
本申请实施例还提供了一种计算机可读存储介质，该计算机可读存储介质存储有计算机指令，该计算机指令被处理器执行时实现本申请实施例提供的遥感影像变化检测方法。
应当注意，尽管在上文详细描述中提及了装置的若干单元或子单元，但是这种划分仅仅是示例性的并非强制性的。实际上，根据本申请的实施方式，上文描述的两个或更多单元的特征和功能可以在一个单元中具体化。反之，上文描述的一个单元的特征和功能可以进一步划分为由多个单元来具体化。
此外，尽管在附图中以特定顺序描述了本申请方法的操作，但是，这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。附加地或备选地，可以省略某些步骤，将多个步骤合并为一个步骤执行，和/或将一个步骤分解为多个步骤执行。
尽管已描述了本申请的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本申请范围的所有变更和修改。
最后应说明的是：以上各实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述各实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。

Claims (10)

1.一种遥感影像变化检测方法，其特征在于，包括：
对两个不同时序的原始影像分别进行影像处理，生成第一地表反射率影像和第二地表反射率影像；
对第一地表反射率影像和第二地表反射率影像进行光谱特征分析，生成光谱差异性影像；
对第一地表反射率影像和第二地表反射率影像进行影像线性变换，生成变化差异影像；
基于光谱差异性影像和变换差异影像，生成基于光谱特征的变化影像；
利用动态变化阈值从变化影像中划分出变化区域。
2.根据权利要求1所述的遥感影像变化检测方法，其特征在于，所述对两个不同时序的原始影像分别进行影像处理，生成第一地表反射率影像和第二地表反射率影像；包括：
利用第一时序的原始影像包含的元数据对第一时序的原始影像进行大气校正，生成第一原始地表反射率影像；利用第二时序的原始数据包含的元数据对第二时序的原始影像进行大气校正，生成第二原始地表反射率影像；
利用第一时序的原始影像包含的质量控制信息对第一原始地表反射率影像进行有效区域掩膜处理，得到第一影像;根据第二时序的原始数据包含的质量控制信息对第二原始地表反射率影像进行有效区域掩膜处理，得到第二影像;
对第一影像和第二影像进行影像配准；
对配准后的第一影像进行重采样和裁切，生成第一地表反射率影像；对配准后的第二影像进行重采样和裁切，生成第二地表反射率影像。
3.根据权利要求1所述的遥感影像变化检测方法，其特征在于，所述对第一地表反射率影像和第二地表反射率影像进行光谱特征分析，生成光谱差异性影像；包括：
获取第一地表反射率影像位置(i,j)处像元光谱向量
Figure P_210918102803578_578782001
，和第二地表反射率影像位置(i,j)处像元光谱向量
Figure P_210918102803641_641305002
，N为波段数，xn为
Figure P_210918102803672_672757003
的第n波段光谱值，yn为
Figure P_210918102803816_816105004
第n波段光谱值；第一地表反射率影像和第二地表反射率影像的大小均为H*W，1≤i≤H，1≤j≤W；H为影像的高度，W为影像的宽度；
计算光谱向量
Figure P_210918102803977_977725001
和光谱向量
Figure P_210918102804008_008976002
之间的光谱角
Figure P_210918102804055_055863003
：
Figure P_210918102804090_090258001
计算光谱向量
Figure P_210918102804168_168188001
和光谱向量
Figure P_210918102804230_230783002
之间的距离
Figure P_210918102804261_261978003
：
Figure P_210918102804313_313163001
计算第一地表反射率影像和第二地表反射率影像的(i,j)位置处像元的光谱差异性值D1(i,j)：
Figure P_210918102804391_391292001
基于各个像元光谱差异性值，生成光谱差异性影像D1。
4.根据权利要求3所述的遥感影像变化检测方法，其特征在于，所述对第一地表反射率影像和第二地表反射率影像进行影像线性变换，生成变化差异影像；包括：
将第一地表反射率影像转换为大小为N*M，均值为0的矩阵FX，其中M=H×W，矩阵FX的每一列为第一地表反射率影像的像元光谱向量；
将第二地表反射率影像转换为大小为N*M，均值为0的矩阵FY，矩阵FY的每一列为第二地表反射率影像的像元光谱向量；
计算FX和FY的协方差
Figure P_210918102804422_422534001
、FY和FX的协方差
Figure P_210918102804532_532414002
、FX的方差
Figure P_210918102804579_579288003
以及FY的方差
Figure P_210918102804626_626134004
；SXY、SYX、SXX和SYY均为N*N的矩阵；
根据相关性系数公式，构建优化目标函数：
Figure P_210918102804673_673050001
其中，相关性系数向量
Figure P_210918102804709_709645001
，所对应的数值为两幅影像N个波段间对应的相关性系数，α为FX的线性变换系数矩阵，大小为N*N；b为FY的线性变换系数矩阵，大小为N*N；
上述目标函数可转换为：
Figure P_210918102804756_756826001
构造Lagrange等式后求导，得到下列公式：
Figure P_210918102804803_803396001
对
Figure P_210918102804834_834664001
进行特征分解，解算出相关性系数向量ρ和矩阵α；
对ρ的各个元素从小到大进行排序，按照该顺序调整α的列向量的顺序；
然后将重排后的ρ和α，代入
Figure P_210918102804882_882958001
中，解算出矩阵b；
计算矩阵X'和矩阵Y'：
Figure P_210918102804914_914727001
利用矩阵X'和矩阵Y'，计算FX和FY之间的变化差异矩阵
Figure P_210918102804946_946434001
：
Figure P_210918102805171_171548001
其中，
Figure P_210918102805218_218443001
为X'第n波段的像元值，大小为1*M;
Figure P_210918102805265_265364002
为Y'第n波段的像元值，大小为1*M；
Figure P_210918102805298_298994003
的大小为1*M；
将
Figure P_210918102805330_330254001
的各个元素按照H*W的影像格式进行重新排列，生成差异影像D2。
5.根据权利要求4所述的遥感影像变化检测方法，其特征在于，所述基于光谱差异性影像和变换差异影像，生成基于光谱特征的变化影像，包括：
对于光谱差异影像D1位置(i,j)处像元值
Figure P_210918102805392_392832001
，以及变化差异影像D2位置(i,j)处像元值
Figure P_210918102805439_439770002
，计算两者之间的乘积，作为变化影像D位置(i,j)处像元值
Figure P_210918102805473_473773003
：
Figure P_210918102805521_521272001
由此生成变化影像D。
6.根据权利要求1-5之一所述的遥感影像变化检测方法，其特征在于，利用动态变化阈值从基于光谱特征的变化影像中划分出变化区域，包括：
对于变化影像的各个像元，如果像元的像元值大于变化阈值，则判断该像元为变化像元，由此划分出变化区域。
7.根据权利要求1所述的遥感影像变化检测方法，其特征在于，所述方法还包括：消除变化区域的孤岛或填补变化区域的空洞。
8.一种遥感影像变化检测装置，其特征在于，包括：
地表反射率影像生成单元，用于对两个不同时序的原始影像分别进行影像处理，生成第一地表反射率影像和第二地表反射率影像；
光谱差异性影像生成单元，用于对第一地表反射率影像和第二地表反射率影像进行光谱特征分析，生成光谱差异性影像；
变化差异影像生成单元，用于对第一地表反射率影像和第二地表反射率影像进行影像线性变换，生成变化差异影像；
变化影像生成单元，用于基于光谱差异性影像和变换差异影像，生成基于光谱特征的变化影像；
变化区域监测单元，用于利用动态变化阈值从变化影像中划分出变化区域。
9.一种电子设备，其特征在于，包括：存储器、处理器和存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如权利要求1-7任一项所述的遥感影像变化检测方法。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器运行时实现如权利要求1-7任一项所述的遥感影像变化检测方法。
CN202111104059.1A 2021-09-22 2021-09-22 遥感影像变化检测方法、装置、设备及存储介质 Active CN113554011B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111104059.1A CN113554011B (zh) 2021-09-22 2021-09-22 遥感影像变化检测方法、装置、设备及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111104059.1A CN113554011B (zh) 2021-09-22 2021-09-22 遥感影像变化检测方法、装置、设备及存储介质

Publications (2)

Publication Number Publication Date
CN113554011A true CN113554011A (zh) 2021-10-26
CN113554011B CN113554011B (zh) 2022-01-04

Family

ID=78134538

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111104059.1A Active CN113554011B (zh) 2021-09-22 2021-09-22 遥感影像变化检测方法、装置、设备及存储介质

Country Status (1)

Country Link
CN (1) CN113554011B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115861823A (zh) * 2023-02-21 2023-03-28 航天宏图信息技术股份有限公司 一种基于自监督深度学习的遥感变化检测方法和装置

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150323666A1 (en) * 2014-05-09 2015-11-12 Nec Corporation Change detection device, change detection method and recording medium
CN106157292A (zh) * 2015-04-23 2016-11-23 尤淑撑 基于两时相遥感图像的国土资源变化监测
CN106548146A (zh) * 2016-11-01 2017-03-29 北京航天泰坦科技股份有限公司 基于时空分析的地表覆盖变化算法和系统
CN108827880A (zh) * 2018-04-23 2018-11-16 吉林大学 基于多光谱影像和ndvi时间序列的地表覆盖变化检测方法
US10325370B1 (en) * 2016-05-31 2019-06-18 University Of New Brunswick Method and system of coregistration of remote sensing images
CA3101681A1 (en) * 2018-06-29 2020-01-02 Farmers Edge Inc. Automatic crop health change detection & alerting system
CN110648347A (zh) * 2019-09-24 2020-01-03 北京航天宏图信息技术股份有限公司 一种基于遥感影像的海岸线提取方法及装置
CN111179230A (zh) * 2019-12-18 2020-05-19 星际空间（天津）科技发展有限公司 遥感影像对比变化检测方法、装置、存储介质及电子设备
CN111931744A (zh) * 2020-10-09 2020-11-13 航天宏图信息技术股份有限公司 一种遥感影像变化检测方法和装置
US20210033589A1 (en) * 2018-03-27 2021-02-04 Flying Gybe Inc. Hyperspectral sensing system and methods

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150323666A1 (en) * 2014-05-09 2015-11-12 Nec Corporation Change detection device, change detection method and recording medium
CN106157292A (zh) * 2015-04-23 2016-11-23 尤淑撑 基于两时相遥感图像的国土资源变化监测
US10325370B1 (en) * 2016-05-31 2019-06-18 University Of New Brunswick Method and system of coregistration of remote sensing images
CN106548146A (zh) * 2016-11-01 2017-03-29 北京航天泰坦科技股份有限公司 基于时空分析的地表覆盖变化算法和系统
US20210033589A1 (en) * 2018-03-27 2021-02-04 Flying Gybe Inc. Hyperspectral sensing system and methods
CN108827880A (zh) * 2018-04-23 2018-11-16 吉林大学 基于多光谱影像和ndvi时间序列的地表覆盖变化检测方法
CA3101681A1 (en) * 2018-06-29 2020-01-02 Farmers Edge Inc. Automatic crop health change detection & alerting system
CN110648347A (zh) * 2019-09-24 2020-01-03 北京航天宏图信息技术股份有限公司 一种基于遥感影像的海岸线提取方法及装置
CN111179230A (zh) * 2019-12-18 2020-05-19 星际空间（天津）科技发展有限公司 遥感影像对比变化检测方法、装置、存储介质及电子设备
CN111931744A (zh) * 2020-10-09 2020-11-13 航天宏图信息技术股份有限公司 一种遥感影像变化检测方法和装置

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
MICHELE VOLPI等: ""Spectral alignment of multi-temporal cross-sensor images with automated kernel canonical correlation analysis"", 《ISPRS JOURNAL OF PHOTOGRAMMETRY AND REMOTE SENSING》 *
TODD A.SCHROEDER等: ""Radiometric correction of multi-temporal Landsat data for characterization of early successional forest patterns in western Oregon"", 《REMOTE SENSING OF ENVIRONMENT》 *
张志新等: ""基于混合像元分解的南方地区植被覆盖度遥感监测——以广州市为例"", 《国土资源遥感》 *
陈曦东等: ""增强型多时相云检测"", 《遥感学报》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115861823A (zh) * 2023-02-21 2023-03-28 航天宏图信息技术股份有限公司 一种基于自监督深度学习的遥感变化检测方法和装置
CN115861823B (zh) * 2023-02-21 2023-05-09 航天宏图信息技术股份有限公司 一种基于自监督深度学习的遥感变化检测方法和装置

Also Published As

Publication number Publication date
CN113554011B (zh) 2022-01-04

Similar Documents

Publication Publication Date Title
Broadwater et al. Hybrid detectors for subpixel targets
Nielsen The regularized iteratively reweighted MAD method for change detection in multi-and hyperspectral data
Wang et al. Low rank constraint and spatial spectral total variation for hyperspectral image mixed denoising
Wu et al. An error-bound-regularized sparse coding for spatiotemporal reflectance fusion
Canty et al. Linear and kernel methods for multivariate change detection
CN115236655B (zh) 基于全极化sar的滑坡识别方法、系统、设备及介质
Li et al. Hyperspectral anomaly detection via image super-resolution processing and spatial correlation
CN113554011B (zh) 遥感影像变化检测方法、装置、设备及存储介质
CN107680081B (zh) 基于卷积神经网络的高光谱图像解混方法
Li et al. Variational Bayesian independent component analysis-support vector machine for remote sensing classification
CN116363521A (zh) 一种遥感影像语义预测方法
Ni et al. Hurricane eye morphology extraction from SAR images by texture analysis
CN112818920B (zh) 一种双时相高光谱图像空谱联合变化检测方法
CN112819769A (zh) 基于核函数和联合字典的非线性高光谱图像异常探测算法
CN113705538A (zh) 基于深度学习的高分遥感影像道路变化检测装置及方法
CN107358625B (zh) 基于SPP Net和感兴趣区域检测的SAR图像变化检测方法
CN110084203B (zh) 基于上下文关联的全卷积网络飞机层级检测方法
CN116758419A (zh) 针对遥感图像的多尺度目标检测方法、装置和设备
CN116310851A (zh) 遥感图像变化检测方法
Wang et al. Multiple Spectral-Spatial Representation Based on Tensor Decomposition for HSI Anomaly Detection
CN113887656B (zh) 一种联合深度学习与稀疏表示的高光谱图像分类方法
CN115730718A (zh) 结合超光谱卫星与人工智能的大气no2时空预测算法
CN113962897B (zh) 基于序列遥感影像的调制传递函数补偿方法及装置
CN113743373A (zh) 基于深度学习的高分遥感影像耕地变化检测装置及方法
Xu et al. The regularized iteratively reweighted object-based MAD method for change detection in bi-temporal, multispectral 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
TR01 Transfer of patent right

Effective date of registration: 20221020

Address after: 100094 building 4, yard 2, Cuihu North Ring Road, Haidian District, Beijing

Patentee after: Aerospace Hongtu Information Technology Co.,Ltd.

Patentee after: Henan Aerospace Hongtu Information Technology Co.,Ltd.

Address before: 100195 Room 301, 3 / F, building 5, zone 4, Xishan Creative Park, Haidian District, Beijing

Patentee before: Aerospace Hongtu Information Technology Co.,Ltd.

TR01 Transfer of patent right