CN107392940A - 一种基于堆栈半监督自适应去噪自编码器的sar图像变化检测方法 - Google Patents

一种基于堆栈半监督自适应去噪自编码器的sar图像变化检测方法 Download PDF

Info

Publication number
CN107392940A
CN107392940A CN201710437638.5A CN201710437638A CN107392940A CN 107392940 A CN107392940 A CN 107392940A CN 201710437638 A CN201710437638 A CN 201710437638A CN 107392940 A CN107392940 A CN 107392940A
Authority
CN
China
Prior art keywords
mrow
msub
layer
image
network
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
CN201710437638.5A
Other languages
English (en)
Other versions
CN107392940B (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.)
Xidian University
Original Assignee
Xidian 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 Xidian University filed Critical Xidian University
Priority to CN201710437638.5A priority Critical patent/CN107392940B/zh
Publication of CN107392940A publication Critical patent/CN107392940A/zh
Application granted granted Critical
Publication of CN107392940B publication Critical patent/CN107392940B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/254—Analysis of motion involving subtraction of images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10044—Radar image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20081—Training; Learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20084—Artificial neural networks [ANN]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20224—Image subtraction

Abstract

本发明公开了一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法，解决了现有方法对相干斑噪声点和较多边缘的变化区域检测精度不高的问题。其实现步骤为：首先生成多尺度差异指导图；将时相1图像作为输入来训练SDAE；将多尺度差异指导图、时相1和时相2图像作为输入来训练SSADAE，SSADAE自适应误差函数中使用了SDAE训练得到的权重；然后用SSADAE计算时相1和时相2图像的特征向量；将两者相减得到差异向量，再对其进行FCM分类，得到变化检测结果图。本发明首先提出了多尺度差异指导图，能够突出差异图中的变化区域；之后提出的SSADAE能够利用图像中少量的标记样本，进一步提高变化检测准确度。

Description

一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检 测方法
技术领域
本发明属于图像处理技术领域，涉及SAR图像的变化检测，具体是一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法。该方法可用在SAR图像的变化检测中。
背景技术
变化检测是遥感领域的关键技术之一，它通过检测同一成像场景中不同时期图像灰度值或局部纹理之间的变化，获取感兴趣的地表或地物在形状、位置、数量及其它属性的变化信息。在社会、环境及军事等领域有广泛的应用。
在多时相SAR图像变化检测方法中，主要的路线有两种，一种是分类后比较(PostClassification Comparison,PCC)，另一种是比较后分类。前一种方法直接对每个时相的SAR图像进行分类，然后对分类的结果进行比较。该方法可以不需要对不同传感器和不同天气条件下拍摄的两幅图像进行辐射校正，但是对不同时相图像像素的分类差异会造成比较结果的巨大差异，且分类方法对变化检测结果影响较大。目前绝大多数的方法都采用先比较后分类的策略。
先比较后分类就是先生成差异图，然后对差异图进行分析。目前差异图生成方法有：
2012年Maoguo Gong和Yu Cao等在IEEE Geoscience and Remote SensingLetters,9卷第2期307-311页发表的文章A Neighborhood-Based Ratio Approach forChange Detection in SAR Images提出了基于邻域的比值算子(NR)，NR算子加入了图像的异质性/同质性算子，但是当变化前后的两幅SAR图像噪声分布不一致时，该方法检测效果不够准确。Jordi Inglada和Grégoire Mercier在IEEE Transactions on Geoscience andRemote Sensing,45卷第5期1432-1445页发表的文章A New Statistical SimilarityMeasure for Change Detection in Multitemporal SAR Images and Its Extension toMultiscale Change Analysis使用KL散度来度量变化前后两幅图像的相似性，具体做法是使用KL散度度量两个独立随机变量间的概率密度函数之间的距离，如果两个变量概率密度值接近KL散度值较小，相反KL散度值会更大。由于这种方法需要顾及变化前后两幅图像的概率密度函数，所以需要较大的计算量来估计概率密度函数的参数，并且需要对图像的分布进行合理的假设。
近年来，深度学习由于能够提取出更加抽象的特征，因此它作为一种新的机器学习方法在语音识别，目标识别,目标检测等应用上都显著的提高了检测结果。而SAR图像由于其乘性噪声的复杂性，基于深度学习的方法能够提取抽象特征的优势给SAR图像变化检测提供了新的思路。目前基于深度学习的变化检测方法有：
2016年Maoguo Gong和Jiaojiao Zhao等在IEEE Transactions on NeuralNetworks and Learning Systems,27卷第1期125-138页发表的文章Change Detection inSynthetic Aperture Radar Images Based on Deep Neural Networks提出了一种基于RBM的无监督变化检测方法。该方法使用联合分类得到初始的变化检测二值图，然后使用样本选择策略对初始的结果进行进一步筛选剔除掉噪声等样本，最后使用该学习样本训练RBM网络。训练后的RBM网络只需要输入变化前后的两时相SAR图像，就可以得到变化检测结果图。该方法省略了生成差异图的步骤，并且由于样本选择策略剔除掉了噪声样本，因此训练后的RBM抗噪性能好；但正因为这种剔除策略同时也剔除了一些边缘信息等有效样本，所以变化检测的漏检率会升高。2016年Puzhao Zhang和Maoguo Gong等在JournalofPhotogrammetry and Remote Sensing,116卷24-41页发表的文章Change detectionbased on deep feature representation and mapping transformation for multi-spatial-resolution remote sensing images结合SDAE特征表示和基于映射的特征变化分析来解决不同分辨率SAR图像的变化检测问题。该方法首先生成了一副粗糙的初始变化掩膜，其中对于同质区使用先比较后分类方法，对于异质区使用先分类后比较方法，然后挑选非变化类作为后续映射神经网络的学习样本；之后使用SDAE分别对两幅不同分辨率图像进行特征学习；将不同分辨率的两幅图像的特征分别作为映射神经网络的输入和输出来学习映射函数，映射函数是这个方法的关键，因为有了映射函数才能真正建立起不同分辨率图像之间的联系，将低分辨率图像映射后的图像就可以直接和高分辨率图像比较，最后使用特征变化分析和FLICM得到最终的变化检测结果。该方法的缺点：一是对光学图像和SAR图像生成初始变化掩膜会引入较大的误差，这个误差较大的学习样本对结果的影响也会较大；二是SDAE是用无监督的方式来提取特征的，所以提取的特征有一定的随意性，而作者提出的基于映射的特征变化分析又是完全依赖于SDAE提取的特征，会进一步引入误差。
综上所述，当变化前后图像噪声分布不一致，上述方法对边缘细节检测效果不够好，变化检测整体错误率较高。
发明内容
本发明为了克服上述已有技术的缺点，提高变化检测的边缘细节等检测效果，有效抑制相干斑噪声，本发明提出了一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法。
本发明的技术方案是：一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法，包括如下步骤：
(1)输入时相1图像I和时相2图像J，图像大小均为U行V列，图像的灰度级范围为0～255；
(2)由时相1图像I和时相2图像J计算多尺度差异指导图MGI；
(3)建立堆栈去噪自编码器SDAE网络，该网络由两层去噪自编码器DAE网络堆栈而成，堆栈结构为：输入层、第一隐含层、第二隐含层、第三隐含层、输出层神经元数目分别为NN1、NN2、NN3、NN2、NN1；在单层预训练阶段，第一层DAE网络的结构为：输入层、隐含层、输出层神经元数目分别为NN1、NN2、NN1；第二层DAE网络的结构为：输入层、隐含层、输出层神经元数目分别为NN2、NN3、NN2；
(4)将时相1图像I作为SDAE网络的输入，进行SDAE网络训练，得到输入层与第一隐含层之间的权重第一隐含层与第二隐含层之间的权重第二隐含层与第三隐含层之间的权重第三隐含层与输出层之间的权重
(5)建立半监督去噪自编码器Semi-DAE网络，网络结构为：输入层、隐含层、输出层神经元数目分别为NN1、NN2、NN1；
(6)建立堆栈半监督自适应去噪自编码器SSADAE网络，该网络由两层网络堆栈而成，堆栈结构与SDAE网络相同，在单层预训练阶段，第一层网络为Semi-DAE网络，第二层网络为DAE网络；
(7)将时相1图像I、时相2图像J和多尺度差异指导图MGI作为SSADAE网络输入，进行SSADAE网络训练，得到SSADAE网络权重，多层微调阶段的误差计算采用自适应误差函数；
(8)将时相1图像I和时相2图像J输入训练好的SSADAE网络得到时相1图像特征和时相2图像特征；将这两个图像特征相减，使用模糊C均值聚类算法FCM分为两类，得到变化检测结果。
步骤(2)中所述的计算多尺度差异指导图MGI的方法包括如下步骤：
(2a)对时相1图像I和时相2图像J中位置(u,v)处像素点的3×3邻域分别计算该3×3邻域中共9个像素值的均值，分别记为μN3(I(u,v))和μN3(J(u,v))，再根据下式计算(u,v)处的3×3邻域均值差异值IS(u,v)，
遍历所有像素点，得到3×3邻域均值差异图像IS＝{IS(u,v)|1≤u≤U,1≤v≤V}，即小尺度差异图IS，这里u和v分别为图像的行序号和列序号。
再根据下式计算(u,v)处的3×3邻域均值差异对数值LIS(u,v)，
LIS(u,v)＝log(IS(u,v)) (2)
遍历所有像素点，得到3×3邻域均值差异对数图像LIS＝{LIS(u,v)|1≤u≤U,1≤v≤V}；
(2b)对时相1图像I和时相2图像J中位置(u,v)处像素点的7×7邻域分别计算该7×7邻域中共49个像素值的均值，分别记为μN7(I(u,v))和μN7(J(u,v))，再根据下式计算(u,v)处两时相之间的7×7邻域均值差异值IM(u,v)，
遍历所有像素点，得到7×7邻域均值差异图像IM＝{IM(u,v)|1≤u≤U,1≤v≤V}，即中尺度差异图IM；
(2c)对时相1图像I和时相2图像J中位置(u,v)处像素点的11×11邻域分别计算该邻域块共121个像素值的均值，分别记为μN11(I(u,v))和μN11(J(u,v))，再根据下式计算(u,v)处两时相之间的11×11邻域均值差异值IL(u,v)，
遍历所有像素点，得到11×11邻域均值差异图像IL＝{IL(u,v)|1≤u≤U,1≤v≤V}，即大尺度差异图IL；
(2d)对步骤(2a)得到的3×3邻域均值差异对数图像LIS，计算其像素点(u,v)的5×5邻域的标准差σN5(u,v)和均值μN5(u,v)的比值r(u,v)，
遍历所有像素点，得到5×5邻域标准差均值比图像R＝{r(u,v)|1≤u≤U,1≤v≤V}；
(2e)按照下式计算(u,v)处两时相之间的多尺度差异参考值MGI(u,v)，
遍历所有像素点，得到多尺度差异指导图像，简称MGI图。
步骤(5)中所述的建立半监督去噪自编码器Semi-DAE网络的具体步骤如下：
(5a)建立Semi-DAE网络，其中输入层神经元到隐单元的权重记为WS1，维数为NN2行NN1列，隐单元的偏置记为bS1，维数为NN2行1列；隐单元到输出层的权重记为WS2，维数为NN1行NN2列，WS2是WS1的转置，输出层神经元的偏置记为bS2，维数为NN1行1列；神经元个数NN1的取值为9，这由输入特征的维数决定的，神经元个数NN2的取值范围是[50,300]，具体NN2取值为100；
对Semi-DAE网络输入样本xS加入噪声后的加噪输入使用如下公式计算隐单元输出hS：
其中，sigm(·)为sigmoid型激励函数，计算公式为sigm(·)＝1/(1+exp(·))。
传统的自编码器直接由隐单元输出hS得到网络输出，为了更有利于突出判别作用较大的神经元，弱化对判别作用较少的神经元，将隐单元输出hS使用转换矩阵A进行映射，映射结果记为
再使用如下公式，计算网络输出zS
其中，转换矩阵A为由步骤(5b)用LPP算法构造得出；
其中LPP算法具体过程如下：
LPP算法目标函数的最小化问题可转变成求解特征值λ和特征向量A的问题：
MLMTA＝λMDMTA (10)
对应于特征值最小的特征向量A就是映射隐单元输出所需的转换矩阵A；
其中，M＝[m1,...,mt,...,mn]，这里mt为输入样本集Φ中具有对应类别标记的输入样本xS对应的隐单元输出，这里M为所有mt拼接成的矩阵，n为对应有类别标记的隐单元输出hS的总数；具有对应类别标记的隐单元输出mt与其对应的类别标记δt构成训练样本集ψ＝{(mt,δt)|1≤t≤n}，输入样本集Φ＝{xS}，输入样本xS可以是具有类别标记的，也可以是不具有类别标记的；δt∈{0,1}，0表示非变化类，1表示变化类；
其中L是拉普拉斯矩阵，L＝D-S；
D是一个对角矩阵，其第i行第i列元素Dii对应于相似性矩阵S中的第j列之和，即Dii＝∑jSij，Sij为相似性矩阵S中第i行第j列的值，其计算公式可简化为如下直接得到图像标记信息的公式：
其中，当第i个和第j个输入样本的类别标记δi和δj一致时，相似性矩阵对应位置取1，其它取0；特别地，当i＝j时，Sij＝1。
步骤(6)中所述的建立堆栈半监督自适应去噪自编码器SSADAE网络的方法具体包括如下步骤：
(6a)SSADAE网络的单层预训练阶段中，第一层网络使用的是步骤(5a)的Semi-DAE网络结构和映射关系，第二层网络使用的是DAE网络；其中DAE网络的网络结构为：输入层神经元个数为NN2，隐单元个数为NN3，输出层神经元个数为NN2；输入层神经元到隐单元的权重记为WD3，维数为NN3行NN2列，隐单元的偏置记为bD3，维数为NN3行1列，隐含层到输出层神经元的权重记为WD4，维数为NN2行NN3列，这里WD4是WD3的转置，输出层神经元的偏置记为bD4，维数为NN2行1列；神经元个数NN3的的取值范围是[10,50]，NN3具体取值为20；
第二层DAE网络的输入为加入噪声后的第一层Semi-DAE网络隐含层输出；对第一层Semi-DAE网络的隐单元输出hS加入噪声后的结果记为本层DAE网络输出记为h′D1；
第二层DAE网络的映射关系如下：
对网络的加噪输入使用如下公式，计算隐单元输出hD2：
再使用如下公式，计算网络输出h′D1：
h′D1＝sigm(WD4·hD2+bD4) (13)
(6b)SSADAE网络的多层微调阶段所建立的堆栈网络：输入层神经元到第一层隐单元的权重记为维数为NN2行NN1列，第一层隐单元的偏置记为维数为NN2行1列，第一层隐单元的输出记为hDD1；第一层隐单元到第二层隐单元的权重记为维数为NN3行NN2列，第二层隐单元的偏置记为维数为NN3行1列，第二层隐单元的输出记为hDD2；第二层隐单元到第三层隐单元的权重记为维数为NN2行NN3列，第三层隐单元的偏置记为维数为NN2行1列，第三层隐单元的输出记为hDD3；第三层隐单元到输出层神经元的权重记为维数为NN1行NN2列，输出层神经元的偏置记为维数为NN1行1列，网络输出记为zD；其中，是的转置，是的转置；
多层微调阶段的网络映射关系如下：
多层微调阶段的网络输入同单层预训练阶段第一层网络的输入相同，对网络加噪后的输入使用如下公式计算第一层隐单元输出hDD1：
第二层隐单元的输出hDD2的计算公式如下：
第三层隐单元的输出hDD3的计算公式如下：
SSADAE网络的输出zS的计算公式如下：
步骤(7)中所述的进行SSADAE网络训练的具体步骤如下：
(7a)SSADAE网络的无标记训练样本集的生成：对时相1图像I、时相2图像J和步骤(2)中得到的MGI图，分别取各图像中对应像素(u,v)的3×3的邻域，将各图像邻域中的9个像素按照先从左到右后从上到下的顺序排成一行，得到各图像的(u,v)像素值向量，分别记为ILu,v、JLu,v和MGILu,v，分别遍历图像I、J和MGI的所有像素点，得到3个9行U×V列的矩阵；将这3个9行U×V列的矩阵沿水平方向拼接，得到9行U×V×3列的矩阵Ξ，即为训练样本集；
(7b)SSADAE网络中Semi-DAE的有标记训练样本集的生成：对步骤(2a)得到的3×3邻域均值差异对数图像LIS，使用FCM分类为变化类和非变化类，从变化类中选择隶属度较高的300个像素位置，从非变化类中选择隶属度较高的300个像素位置；将这600个像素对应的时相1图像I的(u,v)像素值向量ILu,v与时相2图像J的(u,v)像素值向量JLu,v相减，得到差异向量DLu,v＝ILu,v-JLu,v，差异向量DLu,v所对应的类别标记为Lu,v；
(7c)SSADAE网络单层预训练阶段的各层网络权重初始化方法如下：
连接网络的第K层神经元到第K+1层神经元的权重WK,K+1的初始化范围如下所示：
WK,K+1～Q[-ε,ε] (18)
其中
Q[-ε,ε]是在区间(-ε,ε)内均匀分布的随机数，符号～表示WK,K+1初始化为在区间(-ε,ε)内的随机数；HK表示第K层神经元个数，HK+1表示第K+1层神经元个数；
SSADAE网络单层预训练阶段各层的神经元偏置均初始化为0；
(7d)将所有差异向量DLu,v的集合加入masking corruption噪声(即将输入元素随机进行删除)，作为Semi-DAE网络的加噪输入，与权重WS1、偏置bS1代入公式(7)中，计算Semi-DAE网络的隐单元输出，记为{HLu,v}；将差异向量DLu,v对应位置(u,v)的类别标记Lu,v与隐单元输出HLu,v构成的有标记训练样本对记为(HLu,v,Lu,v)，将所有600对有标记训练样本构成的集合{(HLu,v,Lu,v)}(u,v)∈B作为有标记的训练样本集ψ；按照步骤(5b)计算转换矩阵A；
(7e)对步骤(7a)得到的训练样本集Ξ，随机抽取一批100个，抽取后的样本不再重复抽取，所有批次的输入样本xS的集合{xS}＝Ξ；输入样本xS加入masking corruption加性噪声后得到加噪输入
(7f)将加噪输入权重WS1和偏置bS1代入公式(7)中，计算Semi-DAE网络的隐单元输出hS；
(7g)然后将隐单元输出hS和步骤(7d)得到的转换矩阵A代入公式(8)中计算转换后的隐单元输出再使用公式(9)计算Semi-DAE网络输出zS，计算网络输出zS和原始信号xS的误差函数J(xS,zS)，采用欧式距离进行计算，计算公式如下：
使用BP算法对误差函数J(xS,zS)进行最小化，更新权重和偏置；
(7h)采用步骤(7g)更新后的权重和偏置，重复步骤(7d)到步骤(7g)直至全部训练样本三次迭代，通过Semi-DAE网络训练后，得到最终更新的网络权重WS1、WS2和偏置bS1、bS2；
(7i)将随机更换的不同批次的加噪输入步骤(7h)所得的最终的网络权重WS1和偏置bS1代入公式(7)中，计算第一层Semi-DAE网络的隐单元输出hS，并将其加入maskingcorruption噪声，得到加噪的隐单元输出记为
(7j)将加噪的隐单元输出权重WD3和WD4、偏置bD3和bD4代入公式(12)和公式(13)中，计算第二层的DAE网络输出h′D1，然后计算网络输出h′D1和原始信号hS的误差函数J(hS，h′D1)，公式如下：
其中，hS为步骤(7i)计算得到的第一层Semi-DAE网络的隐单元输出；
使用BP算法对误差函数J(hS，h′D1)进行最小化，更新权重WD3、WD4和偏置bD3、bD4；
(7k)重复步骤(7i)和(7j)直至训练样本集Ξ中的全部输入样本三次迭代通过SDAE的网络训练后，得到第二层DAE最终的网络权重WD3、WD4和偏置bD3、bD4；
(7l)SSADAE多层微调阶段的权重和偏置初始化为步骤(7h)和步骤(7k)得到的网络权重和偏置，即
(7m)SSADAE多层微调阶段的输入和第一层Semi-SDAE网络输入相同，为按照步骤(6b)计算多层微调阶段网络输出，记为zS；
(7n)对网络输入样本xS和网络输出zS使用下面的自适应的误差函数:
使用BP算法和梯度下降法最小化误差函数J(xS,zS)，更新多层微调阶段网络权重和偏置；
其中，λ是使网络的误差函数自适应的正则项参数，β是用于保证平稳过渡权重并防止误差函数J(xS,zS)过拟合，其中λ＝0.0001，β＝0.05；
(7o)随机更换不同批次的输入样本的加噪输入，采用(7n)更新后的权重和偏置，重复步骤(7m)到(7n)直至训练样本集ψ中的全部训练样本三次迭代通过SDAE网络多层微调阶段训练后，得到最终的网络权重偏置
本发明的有益效果：本发明提出了一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法，首先对时相1图像和时相2图像生成的不同邻域大小的均值比图像，加入差异图像的异质性和同质性信息，得到多尺度差异指导图，接着对变化前的图像输入堆栈去噪自编码器网络得到自适应权重，然后将时相1图像、时相2图像和多尺度图像加入到本专利提出的堆栈半监督自适应去噪自编码器训练网络权重，最后将时相1图像和时相2图像输入训练好的堆栈半监督自适应去噪自编码器得到时相1图像特征和时相2图像特征，将这两个特征相减，使用FCM分为两类，得到变化检测结果。本发明由于构造了多尺度差异指导图，并建立了堆栈半监督自适应去噪自编码器网络来提取图像特征，所以具有如下优点：
(A)当时相1图像和时相2图像噪声分布不一致时，能有效的检测出变化区域。
(B)能够充分利用少量的标记样本信息，提高变化检测准确度。
(C)对SAR图像存在的相干斑噪声有良好的抑制效果，并能有效的检测边缘等细节信息。
实验证明，本发明能够有效的检测SAR图像边缘等细节信息，对相干斑噪声抑制良好，变化检测具有较低的漏检率和误检率。
附图说明
图1本发明的整体实现流程图；
图1(a)网络训练流程图；
图1(b)变化检测流程图；
图2堆栈去噪自编码器网络结构；
图2(a)SDAE网络单层预训练中第一层训练的网络结构；
图2(b)SDAE网络单层预训练中第二层训练的网络结构；
图2(c)SDAE网络多层微调的网络结构；
图3半监督自编码器网络结构；
图4 Ottawa数据集对应图像；
图5 Bern数据集对应图像；
图6 Ottawa数据集变化检测结果；
图7 Bern数据集变化检测结果；
具体实施方式
参见图1，本发明提供了一种基于无监督的深度神经网络SAR图像变化检测方法，具体包括如下步骤：
步骤1：输入时相1图像I和时相2图像J，I＝{I(u,v)|1≤u≤U,1≤v≤V}，J＝{J(u,v)|1≤u≤U,1≤v≤V}，其中I(u,v)和J(u,v)分别为图像I和图像J在像素点(u,v)的灰度值，这里u和v分别为图像的行序号和列序号，最大行序号为U，最大列序号为V。
步骤2：计算多尺度差异指导图
(2a)对时相1图像I和时相2图像J中位置(u,v)处像素点的3×3邻域，分别计算该3×3邻域中共9个像素值的均值，分别记为μN3(I(u,v))和μN3(J(u,v))，再根据下式计算(u,v)处的3×3邻域均值差异值IS(u,v)，
遍历所有像素点，得到3×3邻域均值差异图像IS＝{IS(u,v)|1≤u≤U,1≤v≤V}，即小尺度差异图IS。
再根据下式计算(u,v)处的3×3邻域均值差异对数值LIS(u,v)，
LIS(u,v)＝log(IS(u,v)) (2)
遍历所有像素点，得到3×3邻域均值差异对数图像LIS＝{LIS(u,v)|1≤u≤U,1≤v≤V}。
(2b)对时相1图像I和时相2图像J中位置(u,v)处像素点的7×7邻域，分别计算该7×7邻域中共49个像素值的均值，分别记为μN7(I(u,v))和μN7(J(u,v))，再根据下式计算(u,v)处两时相之间的7×7邻域均值差异值IM(u,v)，
遍历所有像素点，得到7×7邻域均值差异图像IM＝{IM(u,v)|1≤u≤U,1≤v≤V}，即中尺度差异图IM。
(2c)对时相1图像I和时相2图像J中位置(u,v)处像素点的11×11邻域，分别计算该11×11邻域块共121个像素值的均值，分别记为μN11(I(u,v))和μN11(J(u,v))，再根据下式计算(u,v)处两时相之间的11×11邻域均值差异值IL(u,v)，
遍历所有像素点，得到11×11邻域均值差异图像IL＝{IL(u,v)|1≤u≤U,1≤v≤V}，即大尺度差异图IL。
(2d)对步骤(2a)得到的3×3邻域均值差异对数图像LIS，计算其像素点(u,v)的5×5邻域的标准差σN5(u,v)和均值μN5(u,v)的比值r(u,v)，
遍历所有像素点，得到5×5邻域标准差均值比图像R＝{r(u,v)|1≤u≤U,1≤v≤V}。
(2e)按照下式计算(u,v)处两时相之间的多尺度差异参考值MGI(u,v)，
遍历所有像素点，得到多尺度差异指导图像，简称MGI图。
步骤3：建立堆栈去噪自编码器网络Stacked Denoising Autoencoders，简称SDAE。(参照图2)
(3a)SDAE网络的每一层都是一个去噪自编码器DAE。第一层DAE网络结构为：输入层神经元个数为NN1，隐单元个数为NN2，输出层神经元个数为NN1。输入层神经元到隐单元的权重记为WD1，维数为NN2行NN1列，隐单元的偏置记为bD1，维数为NN2行1列；隐单元到输出层神经元的权重记为WD2，维数为NN1行NN2列，这里WD2是WD1的转置；输出层神经元的偏置记为bD2，维数为NN1行1列。神经元个数NN1的取值为9，这由输入特征的维数决定的，神经元个数NN2的取值范围是[50,300]，本发明实例中，NN2取值为100。
第一层DAE网络的映射关系如下：
对网络输入样本xD加入噪声后的加噪输入使用如下公式，计算隐单元输出hD1：
其中，sigm(·)为sigmoid型激励函数，计算公式为sigm(·)＝1/(1+exp(·))。
再使用如下公式，计算网络输出x′D：
x′D＝sigm(WD2·hD1+bD2) (8)
(3b)第二层DAE的网络结构为：输入层神经元个数为NN2，隐单元个数为NN3，输出层神经元个数为NN2。输入层神经元到隐单元的权重记为WD3，维数为NN3行NN2列，隐单元的偏置记为bD3，维数为NN3行1列，隐含层到输出层神经元的权重记为WD4，维数为NN2行NN3列，这里WD4是WD3的转置，输出层神经元的偏置记为bD4，维数为NN2行1列。神经元个数NN3的的取值范围是[10,50]，本发明实例中，NN3取值为20。
对第一层DAE网络的隐单元输出hD1加入噪声后的结果记为本层DAE网络输出记为h′D1。
第二层DAE网络的映射关系如下：
对网络的加噪输入使用如下公式，计算隐单元输出hD2：
再使用如下公式，计算网络输出h′D1：
h′D1＝sigm(WD4·hD2+bD4) (10)
(3c)对应SDAE网络的多层微调阶段所建立的堆栈网络结构为：输入层神经元个数为NN1，第一层隐单元个数为NN2，第二层隐单元个数为NN3，第三层隐单元个数为NN2，输出层神经元个数为NN1。输入层神经元到第一层隐单元的权重记为WDD1，维数为NN2行NN1列，第一层隐单元的偏置记为bDD1，维数为NN2行1列，第一层隐单元的输出记为hDD1；第一层隐单元到第二层隐单元的权重记为WDD2，维数为NN3行NN2列，第二层隐单元的偏置记为bDD2，维数为NN3行1列，第二层隐单元的输出记为hDD2；第二层隐单元到第三层隐单元的权重记为WDD3，维数为NN2行NN3列，第三层隐单元的偏置记为bDD3，维数为NN2行1列，第三层隐单元的输出记为hDD3；第三层隐单元到输出层神经元的权重记为WDD4，维数为NN1行NN2列，输出层神经元的偏置记为bDD4，维数为NN1行1列，网络输出记为zD。其中，WDD4是WDD1的转置，WDD3是WDD2的转置。
多层微调阶段的网络映射关系如下：
多层微调阶段的网络输入同单层预训练阶段第一层网络的输入相同，即对xD加噪所得使用如下公式计算第一层隐单元输出hDD1：
第二层隐单元的输出hDD2的计算公式如下：
hDD2＝sigm(WDD2·hDD1+bDD2) (12)
第三层隐单元的输出hDD3的计算公式如下：
hDD3＝sigm(WDD3·hDD2+bDD3) (13)
SDAE网络的输出zD的计算公式如下：
zD＝sigm(WDD4·hDD3+bDD4) (14)
步骤4：SDAE网络的训练。分为两个阶段，分别是单层预训练阶段和多层微调阶段，其中单层预训练阶段分第一层训练和第二层训练，就是分别对步骤(3a)和步骤(3b)所建立的DAE网络进行训练。
(4a)SDAE真实训练的样本集生成。取时相1图像I中像素(u,v)的3×3的邻域，将邻域中的9个像素按照先从左到右后从上到下的顺序排成一行，记为ILu,v，称做时相1图像I的(u,v)像素值向量；遍历所有像素点，得到9行U×V列的矩阵，即构成了训练样本集Ж。
(4b)将训练样本集Ж中的样本随机抽取一批100个组成SDAE每次训练的输入样本，记为xD。抽取后的这一批样本不再重复抽取，所有批次的输入样本xD的全体即全部训练样本集Ж＝{xD}。加噪输入是输入样本xD加入加性噪声，如高斯噪声、椒盐噪声等，本发明实例中添加的噪声是masking corruption，即将输入元素随机进行删除。
(4c)权重初始化公式为2010年Glorot Xavier和Yoshua Bengio在Aistats.2010,9:249-256发表的论文Understanding the difficulty of training deep feedforwardneural networks中提出的。连接网络的第K层神经元到第K+1层神经元的权重WK,K+1的初始化范围如下所示：
WK,K+1～Q[-ε,ε] (15)
其中，Q[-ε,ε]是在区间(-ε,ε)内均匀分布的随机数，符号～表示WK,K+1初始化为在区间(-ε,ε)内的随机数。HK表示第K层神经元个数，HK+1表示第K+1层神经元个数。
所有层的神经元偏置均初始化为0。
(4d)将网络加噪输入权重WD1和WD2、偏置bD1和bD2代入步骤(3a)中，计算第一层的去噪自编码器的隐单元输出hD1和网络输出x′D。然后计算网络输出x′D和原始信号xD的误差函数J(xD,x′D)，一般采用交叉熵和欧式距离，本发明专利使用欧式距离的计算公式如下：
使用BP算法对误差函数J(xD,x′D)进行最小化，更新权重WD1、WD2和偏置bD1、bD2。
(4e)随机更换不同批次的输入样本的加噪输入采用步骤(4d)更新后的权重和偏置，重复步骤(4d)直至训练样本集Ж中的全部输入样本迭代三次通过SDAE的第一层网络训练后，得到第一层DAE最终的网络权重WD1、WD2和偏置bD1、bD2。
(4f)将随机更换的不同批次的加噪输入步骤(4e)所得的最终的网络权重WD1、WD2和偏置bD1、bD2代入步骤(3a)的公式(7)中，计算第一层DAE的隐单元输出hD1，并将其加入masking corruption噪声，得到加噪的隐单元输出
(4g)将加噪的隐单元输出权重WD3和WD4、偏置bD3和bD4代入步骤(3b)中，计算第二层的DAE网络输出h′D1。然后计算网络输出h′D1和原始信号hD1的误差函数J(hD1，h′D1)，公式如下：
其中，hD1为步骤(4f)计算得到的第一层DAE的隐单元输出。
使用BP算法对误差函数J(hD1，h′D1)进行最小化，更新权重WD3、WD4和偏置bD3、bD4。
(4h)重复步骤(4f)和(4g)直至训练样本集Ж中的全部输入样本三次迭代通过SDAE的网络训练后，得到第二层DAE最终的网络权重WD3、WD4和偏置bD3、bD4。
(4i)多层微调阶段的堆栈网络的权重和偏置初始化为步骤(4e)和步骤(4h)得到的网络权重和偏置，即WDD1＝WD1，bDD1＝bD1，WDD2＝WD3，bDD2＝bD3，WDD3＝WD4，bDD3＝bD4，WDD4＝WD2，bDD4＝bD2。
(4j)将加噪输入堆栈网络的权重和偏置代入步骤(3c)中，计算堆栈网络输出zD。然后计算堆栈网络输出zD和原始信号xD的误差函数J(xD,zD)，计算公式如下：
使用BP算法对误差函数J(xD,zD)进行最小化，更新多层微调阶段网络权重WDD1、WDD2、WDD3、WDD4和偏置bDD1、bDD2、bDD3、bDD4。
(4k)随机更换不同批次的输入样本的加噪输入采用步骤(4j)更新后的权重和偏置，重复步骤(4j)直至训练样本集Ж中的全部输入样本三次迭代通过SDAE的多层微调阶段训练后，得到最终的网络权重和偏置，分别记为和偏置
步骤5：建立半监督去噪自编码器网络Semi-Supervised DenoisingAutoencoders，简称Semi-DAE(参照图3)。
(5a)单层Semi-DAE网络结构为：输入层神经元个数为NN1，隐单元个数为NN2，输出层神经元个数为NN1。输入层神经元到隐单元的权重记为WS1，维数为NN2行NN1列，隐单元的偏置记为bS1，维数为NN2行1列；隐单元到输出层的权重记为WS2，维数为NN1行NN2列，WS2是WS1的转置，输出层神经元的偏置记为bS2，维数为NN1行1列。
Semi-DAE网络的映射关系如下：
对Semi-DAE网络输入样本xS加入噪声后的加噪输入使用如下公式计算隐单元输出hS：
传统的自编码器直接由隐单元输出hS得到网络输出，本专利为了更有利于突出判别作用较大的神经元，弱化对判别作用较少的神经元，将隐单元输出hS使用转换矩阵A进行映射，映射结果记为
再使用如下公式，计算网络输出zS
其中，转换矩阵A为由步骤(5b)用LPP算法构造的。
(5b)本专利所采用LPP算法为2003年Xiaofei He和Partha Niyogi在NIPS,vol.16发表的文章Locality Preserving Projections里提出的局部保留投影算法。具体如下：
LPP算法目标函数的最小化问题可以转变成求解特征值λ和特征向量A的问题：
MLMTA＝λMDMTA (22)
对应于特征值最小的特征向量A就是映射隐单元输出所需的转换矩阵A。
其中，M＝[m1,...,mt,...,mn]，这里mt为输入样本集Φ中具有对应类别标记的输入样本xS对应的隐单元输出，这里M为所有mt拼接成的矩阵，n为对应有类别标记的隐单元输出hS的总数；具有对应类别标记的隐单元输出mt与其对应的类别标记δt构成有标记训练样本集ψ＝{(mt,δt)|1≤t≤n}；输入样本集Φ＝{xS}，输入样本xS可以是具有类别标记的，也可以是不具有类别标记的。δt∈{0,1}，0表示非变化类，1表示变化类。
L是拉普拉斯矩阵，L＝D-S。
D是一个对角矩阵，其第i行第i列元素Dii对应于相似性矩阵S中的第j列之和，即Dii＝∑jSij，Sij为相似性矩阵S中第i行第j列的值，其计算公式可简化为如下直接得到图像标记信息的公式：
其中，当第i个和第j个输入样本的类别标记δi和δj一致时，相似性矩阵对应位置取1，其它取0。特别地，当i＝j时，Sij＝1。
步骤6：建立堆栈半监督自适应去噪自编码器网络Stacked Semi-SupervisedAdaptive Denoising Autoencoders，简称SSADAE。
(6a)SSADAE网络的单层预训练阶段中，第一层网络使用的是步骤(5a)的Semi-DAE网络结构和映射关系，第二层网络使用的是步骤(3b)的SDAE网络中第二层DAE网络结构和映射关系。
其中，第二层DAE网络的输入为加入噪声后的第一层Semi-DAE网络隐含层输出。
(6b)SSADAE网络的多层微调阶段使用步骤(3c)的网络结构和映射关系。
步骤7：SSADAE网络的真实单层预训练。
(7a)SSADAE网络的无标记训练样本集的生成。对时相1图像I、时相2图像J和步骤2中得到的MGI图，分别取各图像中对应像素(u,v)的3×3的邻域，将各图像邻域中的9个像素按照先从左到右后从上到下的顺序排成一行，得到各图像的(u,v)像素值向量，分别记为ILu,v、JLu,v和MGILu,v，分别遍历图像I、J和MGI的所有像素点，得到3个9行U×V列的矩阵；将这3个9行U×V列的矩阵沿水平方向拼接，得到9行U×V×3列的矩阵Ξ，即为训练样本集；
(7b)SSADAE网络中Semi-DAE的有标记训练样本集的生成。对步骤(2a)得到的3×3邻域均值差异对数图像LIS，使用FCM分类为变化类和非变化类，从变化类中选择隶属度较高的300个像素位置，从非变化类中选择隶属度较高的300个像素位置，将这600个像素对应的时相1图像I的(u,v)像素值向量ILu,v与时相2图像J的(u,v)像素值向量JLu,v相减，得到差异向量DLu,v＝ILu,v-JLu,v，差异向量DLu,v所对应的类别标记为Lu,v；
(7c)按照步骤(4c)来初始化Semi-DAE网络的权重WS1和WS2、偏置bS1和bS2；
(7d)将所有差异向量DLu,v的集合加入masking corruption噪声作为Semi-DAE网络的加噪输入，与权重WS1和偏置bS1代入公式(19)中，计算Semi-DAE网络的隐单元输出，记为{HLu,v}；将差异向量DLu,v对应位置(u,v)的类别标记Lu,v与隐单元输出HLu,v构成的有标记训练样本对记为(HLu,v,Lu,v)，将所有600对有标记训练样本构成的集合{(HLu,v,Lu,v)}(u,v)∈B作为有标记的训练样本集ψ；按照步骤(5b)计算转换矩阵A；
(7e)对步骤(7a)得到的训练样本集Ξ，随机抽取一批100个，抽取后的样本不再重复抽取，所有批次的输入样本xS的集合{xS}＝Ξ；输入样本xS加入masking corruption加性噪声后得到加噪输入
(7f)将加噪输入权重WS1和偏置bS1代入公式(19)中，计算Semi-DAE网络的隐单元输出hS；
(7g)然后将隐单元输出hS和步骤(7d)得到的转换矩阵A代入公式(20)中计算转换后的隐单元输出再使用公式(21)计算Semi-DAE网络输出zS，计算网络输出zS和原始信号xS的误差函数J(xS,zS)，一般采用交叉熵和欧式距离，本发明专利使用欧式距离的计算公式如下：
使用BP算法对误差函数J(xS,zS)最小化，更新权重WS1和WS2、偏置bS1和bS2；
(7h)采用步骤(7g)更新后的权重和偏置，重复步骤(7d)到步骤(7g)直至全部训练样本三次迭代，通过Semi-DAE网络训练后，得到最终更新的网络权重WS1、WS2和偏置bS1、bS2。
(7i)将步骤(7h)所得的最终的网络权重WSR1、WSR2和偏置bSR1、bSR2分别作为步骤(5a)中的权重WS1、WS2和偏置bS1、bS2，代入公式(19)中，计算隐单元输出hS；
(7j)SSADAE单层预训练阶段第二层DAE网络的权重和偏置的初始化公式按照步骤(4c)进行；
(7k)SSADAE单层预训练阶段第二层DAE网络训练如下：用训练样本集Ξ替代训练样本集Ж，并用加噪输入作为用步骤(7h)所得的最终的网络权重WSR1、WSR2和偏置bSR1、bSR2分别替代步骤(4e)中的权重WD1、WD2和偏置bD1、bD2，将步骤(7i)得到的隐单元输出hS加入masking corruption噪声，执行步骤(4g)到(4h)，得到第二层去噪自编码器最终的网络权重和偏置，分别记为WSR3、WSR4和bSR3、bSR4。
步骤8：SSADAE网络的真实多层微调训练。
(8a)SSADAE多层微调阶段的权重和偏置初始化为步骤(7h)和步骤(7k)得到的网络权重和偏置，即
(8b)SSADAE多层微调阶段的输入是将加噪输入作为将网络权重 分别作为WDD1、WDD2、WDD3、WDD4并将网络偏置 分别作为bDD1、bDD2、bDD3、bDD4代入步骤(3c)中计算多层微调阶段网络输出，记为zS。
(8c)对网络输入样本xS和网络输出zS使用下面的自适应的误差函数：
使用BP算法对误差函数J(xS,zS)最小化，更新多层微调阶段网络权重 和偏置
其中，λ是使网络的误差函数自适应的正则项参数，β是用于保证平稳过渡权重并防止误差函数J(xS,zS)过拟合，λ和β采用2014年Jun Deng和Zixing Zhang等在IEEE SignalProcessing Letters第21卷第9期1068-1072页发表的文章Autoencoder-basedUnsupervised Domain Adaptation for Speech Emotion Recognition里提出的值，即λ＝0.0001，β＝0.05。
(8d)随机更换不同批次的输入样本的加噪输入，采用(8c)更新后的权重和偏置，重复步骤(8b)到(8c)直至训练样本集ψ中的全部训练样本三次迭代通过SDAE网络多层微调阶段训练后，得到最终的网络权重偏置
步骤9：用完成训练的SSADAE网络计算时相1图像特征和时相2图像特征，生成差异图，得到变化检测结果。
(9a)计算时相1图像特征：
对步骤(4a)的时相1图像的ILu,v和步骤(8d)得到的权重和偏置采用如下公式计算SSADAE第一层的隐单元的输出，记为IHu,v：
SSADAE第一层隐单元的输出IHu,v作为第二层隐单元的输入，将由步骤(8d)得到的权重和偏置代入如下的公式中，计算SSADAE第二层隐单元的输出，记为IFu,v：
IFu,v为时相1图像对应位置(u,v)的特征；
(9b)计算时相2图像特征：
将步骤(7a)中的时相2图像的JLu,v和步骤(8d)权重和偏置采用如下公式计算SSADAE第一层隐单元的输出JHu,v：
SSADAE第一层隐单元的输出JHu,v作为第二层隐单元的输入，将由步骤(8d)得到的权重和偏置代入如下的公式中，计算SSADAE第二层隐单元的输出JFu,v：
JFu,v为时相2图像对应位置(u,v)的特征；
(9c)计算差异特征：
将时相1图像和时相2图像在位置(u,v)处的特征IFu,v和JFu,v相减，得差异特征DFu,v，即DFu,v＝IFu,v-JFu,v，遍历图像所有位置，得到差异特征矩阵DF。
(9d)对DF使用FCM分类，类别数为2，得到变化检测结果图。
本发明的效果可以通过仿真实验具体说明：
1.实验条件
实验所用微机CPU为Intel Core i7 2.5GHz内存4GB，编程平台是Matlab R2014a。
实验中采用的第一组数据为Ottawa数据集，Radarsat-1SAR图像，大小为350×290，时间分别为1997年5月和1997年8月，如图5所示，其中：
图4(a)是Ottawa数据集第一时相图像；
图4(b)是Ottawa数据集第二时相图像；
图4(c)是Ottawa数据集变化参考图。
第二组数据为Bern数据集，ERS-2SAR图像，大小为301×301，时间分别为2008年6月和2009年6月，如图6所示，其中：
图5(a)是Bern数据集第一时相图像；
图5(b)是Bern数据集第二时相图像；
图5(c)是Bern数据集变化参考图。
2.实验内容
本专利从变化检测结果图和定量指标两个角度来分析。
假设图像中所有的像素个数为ON，实际未发生变化的像素在变化检测后被当做变化类的像素个数为Uc，实际发生变化的像素在变化检测后被当做非变化类的像素个数为Cu，实际发生变化的像素在变化检测后被当做变化类的像素个数为Cc，实际未发生变化的像素在变化检测后被当做非变化类的像素个数为Uu，则定量指标可以表示如下：
(1)MA：漏检数。用于统计实际发生变化的像素在变化检测后被当做非变化类的像素个数。
MA＝Cu
(2)FA：虚警数。用于统计在真实结果图中属于非变化类的像素而在实验结果图中被检测为变化类的像素。
FA＝Uc
(3)OE：整体错误数。漏检数和错检数的求和，用于统计检测类别错误的像素数。
OE＝MA+FA
(4)PCC：统计正确分类所占百分比
(5)KC：KAPPA系数。这个指标常用来统计分类算法的指标，由于包含了更多的细节信息所以比PCC更加准确。图像中真实的非变化类像素个数记为Nu，真实的变化类像素个数记为Nc。
其中
由于本专利方法是基于深度学习的方法，所以本专利选择了三个对比算法，分别是2016年Maoguo Gong和Jiaojiao Zhao等在IEEE Transactions on Neural Networksand Learning Systems,27卷125-138页发表的文章Change Detection in SyntheticAperture Radar Images Based on Deep Neural Networks发表的基于深度学习RBM的算法，简记为RBM。2016年Yawei Li和Lizuo Jin等在International Joint Conference onNeural Networks上4032-4039页发表的文章Semi-supervised Auto-encoder Based onManifold Learning，简记为semAE。
此外，实验结果还比较了本文提出的SSADAE和原始的SDAE、自适应的SDAE(A-SDAE)方法比较，这三种网络每层神经元个数均相同，其中SDAE方法是步骤3介绍的堆栈去噪自编码器网络，由两个DAE网络堆栈组成；自适应的SDAE方法依然有两个DAE网络堆栈组成，但是目标函数采用步骤9公式(25)介绍自适应的误差函数。
图6为Ottawa数据集变化检测结果，其中图6(a)Ottawa数据集变化参考图；图6(b)Ottawa数据集使用RBM方法得到的变化检测结果；图6(c)Ottawa数据集使用semAE方法得到的变化检测结果；图6(d)Ottawa数据集使用SDAE方法得到的变化检测结果；图6(e)Ottawa数据集使用A-SDAE方法得到的变化检测结果；图6(f)Ottawa数据集使用本专利方法得到的变化检测结果。
图7为Bern数据集变化检测结果，其中图7(a)Bern数据集变化参考图；图7(b)Bern数据集使用RBM方法得到的变化检测结果；图7(c)Bern数据集使用semAE方法得到的变化检测结果；图7(d)Bern数据集使用SDAE方法得到的变化检测结果；图7(e)Bern数据集使用A-SDAE方法得到的变化检测结果；图7(f)Bern数据集使用本专利方法得到的变化检测结果。
3.实验结果
总结分析：从表1可以看到，对比RBM方法和semAE，本专利方法无论在漏检数、误检数和整体错误率方面表现都最好；表1还比较了本专利方法和改进之前的SDAE，和自适应的SDAE结果进行比较，可以看出来A-SDAE在漏检数方面比SDAE要好，但是A-SDAE方法误检数高，而本专利提出的SSADAE方法误检数比A-SDAE方法要低的多，此外本专利方法在整体错误率和Kappa系数方面表现也最好。
表1变化检测结果对比
从图6(f)和图6(b)、图6(c)、图6(d)、图6(e)比较，可以看到本专利方法对Ottawa左上角的水域检测的效果最好。从图7(f)和图7(b)、图7(c)、图7(d)、图7(e)比较，可以看到本专利方法不仅检测的噪声点少，此外对Bern边缘检测效果也好。
表2两组SAR图像有MGI图和没有MGI图结果比较
表2记录本专利提出的SSADAE迭代三次时，加入MGI图和不加入MGI图的结果比较。之所以选择迭代到第三次的时候，是因为迭代到三次时，SSADAE训练加入MGI和SDAE训练不加MGI基本都达到稳定，所以取达到稳定的结果比较更有说服力。从表2所示，可以看到无论从漏检率，整体错误率，还有KAPPA系数方面有MGI都比没有MGI表现要更好。
本实施方式中没有详细叙述的部分属本行业的公知的常用手段，这里不一一叙述。以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。

Claims (5)

1.一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法，其特征在于，包括如下步骤：
(1)输入时相1图像I和时相2图像J，图像大小均为U行V列，图像的灰度级范围为0～255；
(2)由时相1图像I和时相2图像J计算多尺度差异指导图MGI；
(3)建立堆栈去噪自编码器SDAE网络，该网络由两层去噪自编码器DAE网络堆栈而成，堆栈结构为：输入层、第一隐含层、第二隐含层、第三隐含层、输出层神经元数目分别为NN1、NN2、NN3、NN2、NN1；在单层预训练阶段，第一层DAE网络的结构为：输入层、隐含层、输出层神经元数目分别为NN1、NN2、NN1；第二层DAE网络的结构为：输入层、隐含层、输出层神经元数目分别为NN2、NN3、NN2；
(4)将时相1图像I作为SDAE网络的输入，进行SDAE网络训练，得到输入层与第一隐含层之间的权重第一隐含层与第二隐含层之间的权重第二隐含层与第三隐含层之间的权重第三隐含层与输出层之间的权重
(5)建立半监督去噪自编码器Semi-DAE网络，网络结构为：输入层、隐含层、输出层神经元数目分别为NN1、NN2、NN1；
(6)建立堆栈半监督自适应去噪自编码器SSADAE网络，该网络由两层网络堆栈而成，堆栈结构与SDAE网络相同，在单层预训练阶段，第一层网络为Semi-DAE网络，第二层网络为DAE网络；
(7)将时相1图像I、时相2图像J和多尺度差异指导图MGI作为SSADAE网络输入，进行SSADAE网络训练，得到SSADAE网络权重，多层微调阶段的误差计算采用自适应误差函数；
(8)将时相1图像I和时相2图像J输入训练好的SSADAE网络得到时相1图像特征和时相2图像特征；将这两个图像特征相减，使用模糊C均值聚类算法FCM分为两类，得到变化检测结果。
2.根据权利要求1所述的一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法，其特征在于，步骤(2)中所述的计算多尺度差异指导图MGI的方法包括如下步骤：
(2a)对时相1图像I和时相2图像J中位置(u,v)处像素点的3×3邻域分别计算该3×3邻域中共9个像素值的均值，分别记为μN3(I(u,v))和μN3(J(u,v))，再根据下式计算(u,v)处的3×3邻域均值差异值IS(u,v)，
<mrow> <msub> <mi>I</mi> <mi>S</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>-</mo> <mi>min</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>3</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>I</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>3</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>J</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>3</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>J</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>3</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>I</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
遍历所有像素点，得到3×3邻域均值差异图像IS＝{IS(u,v)|1≤u≤U,1≤v≤V}，即小尺度差异图IS，这里u和v分别为图像的行序号和列序号。
再根据下式计算(u,v)处的3×3邻域均值差异对数值LIS(u,v)，
LIS(u,v)＝log(IS(u,v)) (2)
遍历所有像素点，得到3×3邻域均值差异对数图像LIS＝{LIS(u,v)|1≤u≤U,1≤v≤V}；
(2b)对时相1图像I和时相2图像J中位置(u,v)处像素点的7×7邻域分别计算该7×7邻域中共49个像素值的均值，分别记为μN7(I(u,v))和μN7(J(u,v))，再根据下式计算(u,v)处两时相之间的7×7邻域均值差异值IM(u,v)，
<mrow> <msub> <mi>I</mi> <mi>M</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>-</mo> <mi>min</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>7</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>I</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>7</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>J</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>7</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>J</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>7</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>I</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
遍历所有像素点，得到7×7邻域均值差异图像IM＝{IM(u,v)|1≤u≤U,1≤v≤V}，即中尺度差异图IM；
(2c)对时相1图像I和时相2图像J中位置(u,v)处像素点的11×11邻域分别计算该邻域块共121个像素值的均值，分别记为μN11(I(u,v))和μN11(J(u,v))，再根据下式计算(u,v)处两时相之间的11×11邻域均值差异值IL(u,v)，
<mrow> <msub> <mi>I</mi> <mi>L</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>-</mo> <mi>min</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>11</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>I</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>11</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>J</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> <mo>,</mo> <mfrac> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>11</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>J</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>11</mn> </mrow> </msub> <mrow> <mo>(</mo> <mi>I</mi> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
遍历所有像素点，得到11×11邻域均值差异图像IL＝{IL(u,v)|1≤u≤U,1≤v≤V}，即大尺度差异图IL；
(2d)对步骤(2a)得到的3×3邻域均值差异对数图像LIS，计算其像素点(u,v)的5×5邻域的标准差σN5(u,v)和均值μN5(u,v)的比值r(u,v)，
<mrow> <mi>r</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&amp;sigma;</mi> <mrow> <mi>N</mi> <mn>5</mn> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>LI</mi> <mi>S</mi> </msub> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msub> <mi>&amp;mu;</mi> <mrow> <mi>N</mi> <mn>5</mn> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>LI</mi> <mi>S</mi> </msub> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
遍历所有像素点，得到5×5邻域标准差均值比图像R＝{r(u,v)|1≤u≤U,1≤v≤V}；
(2e)按照下式计算(u,v)处两时相之间的多尺度差异参考值MGI(u,v)，
<mrow> <mtable> <mtr> <mtd> <mrow> <mi>M</mi> <mi>G</mi> <mi>I</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>&amp;lsqb;</mo> <msub> <mi>I</mi> <mi>S</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>I</mi> <mi>M</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>&amp;times;</mo> <mi>r</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>+</mo> <mo>&amp;lsqb;</mo> <msub> <mi>I</mi> <mi>M</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>I</mi> <mi>L</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>&amp;lsqb;</mo> <mn>1</mn> <mo>-</mo> <mi>r</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>=</mo> <mo>&amp;lsqb;</mo> <msub> <mi>I</mi> <mi>S</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>I</mi> <mi>L</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>&amp;times;</mo> <mi>r</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>I</mi> <mi>M</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>I</mi> <mi>L</mi> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
遍历所有像素点，得到多尺度差异指导图像，简称MGI图。
3.根据权利要求2所述的一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法，其特征在于，步骤(5)中所述的建立半监督去噪自编码器Semi-DAE网络的具体步骤如下：
(5a)建立Semi-DAE网络，其中输入层神经元到隐单元的权重记为WS1，维数为NN2行NN1列，隐单元的偏置记为bS1，维数为NN2行1列；隐单元到输出层的权重记为WS2，维数为NN1行NN2列，WS2是WS1的转置，输出层神经元的偏置记为bS2，维数为NN1行1列；神经元个数NN1的取值为9，这由输入特征的维数决定的，神经元个数NN2的取值范围是[50,300]，具体NN2取值为100；
对Semi-DAE网络输入样本xS加入噪声后的加噪输入使用如下公式计算隐单元输出hS：
<mrow> <msub> <mi>h</mi> <mi>S</mi> </msub> <mo>=</mo> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>m</mi> <mrow> <mo>(</mo> <msub> <mi>W</mi> <mrow> <mi>S</mi> <mn>1</mn> </mrow> </msub> <mo>&amp;CenterDot;</mo> <msub> <mover> <mi>x</mi> <mo>~</mo> </mover> <mi>S</mi> </msub> <mo>+</mo> <msub> <mi>b</mi> <mrow> <mi>S</mi> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
其中，sigm(·)为sigmoid型激励函数，计算公式为sigm(·)＝1/(1+exp(·))。
传统的自编码器直接由隐单元输出hS得到网络输出，为了更有利于突出判别作用较大的神经元，弱化对判别作用较少的神经元，将隐单元输出hS使用转换矩阵A进行映射，映射结果记为
<mrow> <msub> <mover> <mi>h</mi> <mo>^</mo> </mover> <mi>S</mi> </msub> <mo>=</mo> <mi>A</mi> <mo>&amp;CenterDot;</mo> <msub> <mi>h</mi> <mi>S</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>
再使用如下公式，计算网络输出zS
<mrow> <msub> <mi>z</mi> <mi>S</mi> </msub> <mo>=</mo> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>m</mi> <mrow> <mo>(</mo> <msub> <mi>W</mi> <mrow> <mi>S</mi> <mn>2</mn> </mrow> </msub> <mo>&amp;CenterDot;</mo> <msub> <mover> <mi>h</mi> <mo>^</mo> </mover> <mi>S</mi> </msub> <mo>+</mo> <msub> <mi>b</mi> <mrow> <mi>S</mi> <mn>2</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
其中，转换矩阵A为由步骤(5b)用LPP算法构造得出；
其中LPP算法具体过程如下：
LPP算法目标函数的最小化问题可转变成求解特征值λ和特征向量A的问题：
MLMTA＝λMDMTA (10)
对应于特征值最小的特征向量A就是映射隐单元输出所需的转换矩阵A；
其中，M＝[m1,...,mt,...,mn]，这里mt为输入样本集Φ中具有对应类别标记的输入样本xS对应的隐单元输出，这里M为所有mt拼接成的矩阵，n为对应有类别标记的隐单元输出hS的总数；具有对应类别标记的隐单元输出mt与其对应的类别标记δt构成训练样本集ψ＝{(mt,δt)|1≤t≤n}，输入样本集Φ＝{xS}，输入样本xS可以是具有类别标记的，也可以是不具有类别标记的；δt∈{0,1}，0表示非变化类，1表示变化类；
其中L是拉普拉斯矩阵，L＝D-S；
D是一个对角矩阵，其第i行第i列元素Dii对应于相似性矩阵S中的第j列之和，即Dii＝∑jSij，Sij为相似性矩阵S中第i行第j列的值，其计算公式可简化为如下直接得到图像标记信息的公式：
<mrow> <msub> <mi>S</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mn>1</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <msub> <mi>&amp;delta;</mi> <mi>i</mi> </msub> <mo>=</mo> <msub> <mi>&amp;delta;</mi> <mi>j</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>0</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>e</mi> <mi>l</mi> <mi>s</mi> <mi>e</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
其中，当第i个和第j个输入样本的类别标记δi和δj一致时，相似性矩阵对应位置取1，其它取0；特别地，当i＝j时，Sij＝1。
4.根据权利要求3所述的一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法，其特征在于，步骤(6)中所述的建立堆栈半监督自适应去噪自编码器SSADAE网络的方法具体包括如下步骤：
(6a)SSADAE网络的单层预训练阶段中，第一层网络使用的是步骤(5a)的Semi-DAE网络结构和映射关系，第二层网络使用的是DAE网络；其中DAE网络的网络结构为：输入层神经元个数为NN2，隐单元个数为NN3，输出层神经元个数为NN2；输入层神经元到隐单元的权重记为WD3，维数为NN3行NN2列，隐单元的偏置记为bD3，维数为NN3行1列，隐含层到输出层神经元的权重记为WD4，维数为NN2行NN3列，这里WD4是WD3的转置，输出层神经元的偏置记为bD4，维数为NN2行1列；神经元个数NN3的的取值范围是[10,50]，NN3具体取值为20；
第二层DAE网络的输入为加入噪声后的第一层Semi-DAE网络隐含层输出；对第一层Semi-DAE网络的隐单元输出hS加入噪声后的结果记为本层DAE网络输出记为h′D1；
第二层DAE网络的映射关系如下：
对网络的加噪输入使用如下公式，计算隐单元输出hD2：
<mrow> <msub> <mi>h</mi> <mrow> <mi>D</mi> <mn>2</mn> </mrow> </msub> <mo>=</mo> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>m</mi> <mrow> <mo>(</mo> <msub> <mi>W</mi> <mrow> <mi>D</mi> <mn>3</mn> </mrow> </msub> <mo>&amp;CenterDot;</mo> <msub> <mover> <mi>h</mi> <mo>~</mo> </mover> <mi>S</mi> </msub> <mo>+</mo> <msub> <mi>b</mi> <mrow> <mi>D</mi> <mn>3</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>
再使用如下公式，计算网络输出h′D1：
h′D1＝sigm(WD4·hD2+bD4) (13)
(6b)SSADAE网络的多层微调阶段所建立的堆栈网络：输入层神经元到第一层隐单元的权重记为维数为NN2行NN1列，第一层隐单元的偏置记为维数为NN2行1列，第一层隐单元的输出记为hDD1；第一层隐单元到第二层隐单元的权重记为维数为NN3行NN2列，第二层隐单元的偏置记为维数为NN3行1列，第二层隐单元的输出记为hDD2；第二层隐单元到第三层隐单元的权重记为维数为NN2行NN3列，第三层隐单元的偏置记为维数为NN2行1列，第三层隐单元的输出记为hDD3；第三层隐单元到输出层神经元的权重记为维数为NN1行NN2列，输出层神经元的偏置记为维数为NN1行1列，网络输出记为zD；其中，是的转置，是的转置；
多层微调阶段的网络映射关系如下：
多层微调阶段的网络输入同单层预训练阶段第一层网络的输入相同，对网络加噪后的输入使用如下公式计算第一层隐单元输出hDD1：
<mrow> <msub> <mi>h</mi> <mrow> <mi>D</mi> <mi>D</mi> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>m</mi> <mrow> <mo>(</mo> <msubsup> <mi>W</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>&amp;CenterDot;</mo> <msub> <mover> <mi>x</mi> <mo>~</mo> </mover> <mi>S</mi> </msub> <mo>+</mo> <msubsup> <mi>b</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>
第二层隐单元的输出hDD2的计算公式如下：
<mrow> <msub> <mi>h</mi> <mrow> <mi>D</mi> <mi>D</mi> <mn>2</mn> </mrow> </msub> <mo>=</mo> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>m</mi> <mrow> <mo>(</mo> <msubsup> <mi>W</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>&amp;CenterDot;</mo> <msub> <mi>h</mi> <mrow> <mi>D</mi> <mi>D</mi> <mn>1</mn> </mrow> </msub> <mo>+</mo> <msubsup> <mi>b</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>
第三层隐单元的输出hDD3的计算公式如下：
<mrow> <msub> <mi>h</mi> <mrow> <mi>D</mi> <mi>D</mi> <mn>3</mn> </mrow> </msub> <mo>=</mo> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>m</mi> <mrow> <mo>(</mo> <msubsup> <mi>W</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>&amp;CenterDot;</mo> <msub> <mi>h</mi> <mrow> <mi>D</mi> <mi>D</mi> <mn>2</mn> </mrow> </msub> <mo>+</mo> <msubsup> <mi>b</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow>
SSADAE网络的输出zS的计算公式如下：
<mrow> <msub> <mi>z</mi> <mi>S</mi> </msub> <mo>=</mo> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>m</mi> <mrow> <mo>(</mo> <msubsup> <mi>W</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>&amp;CenterDot;</mo> <msub> <mi>h</mi> <mrow> <mi>D</mi> <mi>D</mi> <mn>3</mn> </mrow> </msub> <mo>+</mo> <msubsup> <mi>b</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>17</mn> <mo>)</mo> </mrow> <mo>.</mo> </mrow>
5.根据权利要求4所述的一种基于堆栈半监督自适应去噪自编码器的SAR图像变化检测方法，其特征在于，步骤(7)中所述的进行SSADAE网络训练的具体步骤如下：
(7a)SSADAE网络的无标记训练样本集的生成：对时相1图像I、时相2图像J和步骤(2)中得到的MGI图，分别取各图像中对应像素(u,v)的3×3的邻域，将各图像邻域中的9个像素按照先从左到右后从上到下的顺序排成一行，得到各图像的(u,v)像素值向量，分别记为ILu,v、JLu,v和MGILu,v，分别遍历图像I、J和MGI的所有像素点，得到3个9行U×V列的矩阵；将这3个9行U×V列的矩阵沿水平方向拼接，得到9行U×V×3列的矩阵Ξ，即为训练样本集；
(7b)SSADAE网络中Semi-DAE的有标记训练样本集的生成：对步骤(2a)得到的3×3邻域均值差异对数图像LIS，使用FCM分类为变化类和非变化类，从变化类中选择隶属度较高的300个像素位置，从非变化类中选择隶属度较高的300个像素位置；将这600个像素对应的时相1图像I的(u,v)像素值向量ILu,v与时相2图像J的(u,v)像素值向量JLu,v相减，得到差异向量DLu,v＝ILu,v-JLu,v，差异向量DLu,v所对应的类别标记为Lu,v；
(7c)SSADAE网络单层预训练阶段的各层网络权重初始化方法如下：
连接网络的第K层神经元到第K+1层神经元的权重WK,K+1的初始化范围如下所示：
WK,K+1～Q[-ε,ε] (18)
其中
Q[-ε,ε]是在区间(-ε,ε)内均匀分布的随机数，符号～表示WK,K+1初始化为在区间(-ε,ε)内的随机数；HK表示第K层神经元个数，HK+1表示第K+1层神经元个数；
SSADAE网络单层预训练阶段各层的神经元偏置均初始化为0；
(7d)将所有差异向量DLu,v的集合加入masking corruption噪声(即将输入元素随机进行删除)，作为Semi-DAE网络的加噪输入，与权重WS1、偏置bS1代入公式(7)中，计算Semi-DAE网络的隐单元输出，记为{HLu,v}；将差异向量DLu,v对应位置(u,v)的类别标记Lu,v与隐单元输出HLu,v构成的有标记训练样本对记为(HLu,v,Lu,v)，将所有600对有标记训练样本构成的集合{(HLu,v,Lu,v)}(u,v)∈B作为有标记的训练样本集ψ；按照步骤(5b)计算转换矩阵A；
(7e)对步骤(7a)得到的训练样本集Ξ，随机抽取一批100个，抽取后的样本不再重复抽取，所有批次的输入样本xS的集合{xS}＝Ξ；输入样本xS加入masking corruption加性噪声后得到加噪输入
(7f)将加噪输入权重WS1和偏置bS1代入公式(7)中，计算Semi-DAE网络的隐单元输出hS；
(7g)然后将隐单元输出hS和步骤(7d)得到的转换矩阵A代入公式(8)中计算转换后的隐单元输出再使用公式(9)计算Semi-DAE网络输出zS，计算网络输出zS和原始信号xS的误差函数J(xS,zS)，采用欧式距离进行计算，计算公式如下：
<mrow> <mi>J</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>S</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>S</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <munder> <mo>&amp;Sigma;</mo> <mrow> <msub> <mi>x</mi> <mi>S</mi> </msub> <mo>&amp;Element;</mo> <mi>&amp;Xi;</mi> </mrow> </munder> <mo>|</mo> <mo>|</mo> <msub> <mi>x</mi> <mi>S</mi> </msub> <mo>-</mo> <msub> <mi>z</mi> <mi>S</mi> </msub> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>19</mn> <mo>)</mo> </mrow> </mrow>
使用BP算法对误差函数J(xS,zS)进行最小化，更新权重和偏置；
(7h)采用步骤(7g)更新后的权重和偏置，重复步骤(7d)到步骤(7g)直至全部训练样本三次迭代，通过Semi-DAE网络训练后，得到最终更新的网络权重WS1、WS2和偏置bS1、bS2；
(7i)将随机更换的不同批次的加噪输入步骤(7h)所得的最终的网络权重WS1和偏置bS1代入公式(7)中，计算第一层Semi-DAE网络的隐单元输出hS，并将其加入maskingcorruption噪声，得到加噪的隐单元输出记为
(7j)将加噪的隐单元输出权重WD3和WD4、偏置bD3和bD4代入公式(12)和公式(13)中，计算第二层的DAE网络输出h′D1，然后计算网络输出h′D1和原始信号hS的误差函数J(hS，h′D1)，公式如下：
<mrow> <mi>J</mi> <mrow> <mo>(</mo> <msub> <mi>h</mi> <mi>S</mi> </msub> <mo>,</mo> <msubsup> <mi>h</mi> <mrow> <mi>D</mi> <mn>1</mn> </mrow> <mo>&amp;prime;</mo> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mi>&amp;Sigma;</mi> <mo>|</mo> <mo>|</mo> <msub> <mi>h</mi> <mi>S</mi> </msub> <mo>-</mo> <msubsup> <mi>h</mi> <mrow> <mi>D</mi> <mn>1</mn> </mrow> <mo>&amp;prime;</mo> </msubsup> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>
其中，hS为步骤(7i)计算得到的第一层Semi-DAE网络的隐单元输出；
使用BP算法对误差函数J(hS，h′D1)进行最小化，更新权重WD3、WD4和偏置bD3、bD4；
(7k)重复步骤(7i)和(7j)直至训练样本集Ξ中的全部输入样本三次迭代通过SDAE的网络训练后，得到第二层DAE最终的网络权重WD3、WD4和偏置bD3、bD4；
(7l)SSADAE多层微调阶段的权重和偏置初始化为步骤(7h)和步骤(7k)得到的网络权重和偏置，即
(7m)SSADAE多层微调阶段的输入和第一层Semi-SDAE网络输入相同，为按照步骤(6b)计算多层微调阶段网络输出，记为zS；
(7n)对网络输入样本xS和网络输出zS使用下面的自适应的误差函数:
<mrow> <mi>J</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>S</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>S</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mi>&amp;lambda;</mi> <mn>2</mn> </mfrac> <mrow> <mo>(</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>4</mn> </munderover> <mo>|</mo> <mo>|</mo> <msubsup> <mi>W</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;beta;W</mi> <mi>D</mi> <mrow> <mi>t</mi> <mi>e</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> </mrow> </msubsup> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <munder> <mo>&amp;Sigma;</mo> <mrow> <msub> <mi>x</mi> <mi>S</mi> </msub> <mo>&amp;Element;</mo> <mi>&amp;Xi;</mi> </mrow> </munder> <mo>|</mo> <mo>|</mo> <msub> <mi>x</mi> <mi>S</mi> </msub> <mo>-</mo> <msub> <mi>z</mi> <mi>S</mi> </msub> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>21</mn> <mo>)</mo> </mrow> </mrow>
使用BP算法和梯度下降法最小化误差函数J(xS,zS)，更新多层微调阶段网络权重和偏置；
其中，λ是使网络的误差函数自适应的正则项参数，β是用于保证平稳过渡权重并防止误差函数J(xS,zS)过拟合，其中λ＝0.0001，β＝0.05；
(7o)随机更换不同批次的输入样本的加噪输入，采用(7n)更新后的权重和偏置，重复步骤(7m)到(7n)直至训练样本集ψ中的全部训练样本三次迭代通过SDAE网络多层微调阶段训练后，得到最终的网络权重偏置
CN201710437638.5A 2017-06-12 2017-06-12 一种基于堆栈半监督自适应去噪自编码器的sar图像变化检测方法 Active CN107392940B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710437638.5A CN107392940B (zh) 2017-06-12 2017-06-12 一种基于堆栈半监督自适应去噪自编码器的sar图像变化检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710437638.5A CN107392940B (zh) 2017-06-12 2017-06-12 一种基于堆栈半监督自适应去噪自编码器的sar图像变化检测方法

Publications (2)

Publication Number Publication Date
CN107392940A true CN107392940A (zh) 2017-11-24
CN107392940B CN107392940B (zh) 2020-10-30

Family

ID=60332466

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710437638.5A Active CN107392940B (zh) 2017-06-12 2017-06-12 一种基于堆栈半监督自适应去噪自编码器的sar图像变化检测方法

Country Status (1)

Country Link
CN (1) CN107392940B (zh)

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108363382A (zh) * 2018-02-09 2018-08-03 哈尔滨工业大学（威海） 一种复杂装备故障诊断方法及系统
CN108492298A (zh) * 2018-04-13 2018-09-04 西安电子科技大学 基于生成对抗网络的多光谱图像变化检测方法
CN108596204A (zh) * 2018-03-15 2018-09-28 西安电子科技大学 一种基于改进型scdae的半监督调制方式分类模型的方法
CN109191418A (zh) * 2018-06-22 2019-01-11 西安电子科技大学 一种基于收缩自编码器特征学习的遥感图像变化检测方法
CN109377480A (zh) * 2018-09-27 2019-02-22 中国电子科技集团公司第五十四研究所 基于深度学习的耕地变化检测方法
CN109753896A (zh) * 2018-12-20 2019-05-14 西安电子科技大学 一种基于共性自动编码器的无监督异源遥感图像变化检测方法
CN110766084A (zh) * 2019-10-28 2020-02-07 北京理工大学 基于cae与hl-cnn的小样本sar目标识别方法
CN111046887A (zh) * 2018-10-15 2020-04-21 华北电力大学（保定） 一种带噪图像特征提取的方法
CN111931035A (zh) * 2019-05-13 2020-11-13 中国移动通信集团湖北有限公司 业务推荐方法、装置及设备
CN112836736A (zh) * 2021-01-28 2021-05-25 哈尔滨理工大学 一种基于深度自编码器构图的高光谱图像半监督分类方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103810699A (zh) * 2013-12-24 2014-05-21 西安电子科技大学 基于无监督深度神经网络的sar图像变化检测方法
CN104778717A (zh) * 2015-05-05 2015-07-15 西安电子科技大学 基于导向差异图的sar图像变化检测方法
CN105844279A (zh) * 2016-03-22 2016-08-10 西安电子科技大学 基于深度学习和sift特征的sar图像变化检测方法
CN105913090A (zh) * 2016-04-14 2016-08-31 西安电子科技大学 基于sdae-svm的sar图像目标分类方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103810699A (zh) * 2013-12-24 2014-05-21 西安电子科技大学 基于无监督深度神经网络的sar图像变化检测方法
CN104778717A (zh) * 2015-05-05 2015-07-15 西安电子科技大学 基于导向差异图的sar图像变化检测方法
CN105844279A (zh) * 2016-03-22 2016-08-10 西安电子科技大学 基于深度学习和sift特征的sar图像变化检测方法
CN105913090A (zh) * 2016-04-14 2016-08-31 西安电子科技大学 基于sdae-svm的sar图像目标分类方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
GUITING WANG ET AL: "A multiscale region-based approach to automatic SAR image registration using CLPSO", 《2014 IEEE GEOSCIENCE AND REMOTE SENSING SYMPOSIUM》 *
PASCAL VINCENT ET AL: "Stacked Denoising Autoencoders:Learning Useful Representations in a Deep Network with a Local Denoising Criterion", 《JOURNAL OF MACHINE LEARNING RESEARCH》 *
YAWEI LI ET AL: "Semi-supervised auto-encoder based on manifold learning", 《2016 INTERNATIONAL JOINT CONFERENCE ON NEURAL NETWORKS (IJCNN)》 *
赵飞翔 等: "基于栈式降噪稀疏自动编码器的雷达目标识别方法", 《雷达学报》 *

Cited By (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108363382A (zh) * 2018-02-09 2018-08-03 哈尔滨工业大学（威海） 一种复杂装备故障诊断方法及系统
CN108596204A (zh) * 2018-03-15 2018-09-28 西安电子科技大学 一种基于改进型scdae的半监督调制方式分类模型的方法
CN108596204B (zh) * 2018-03-15 2021-11-09 西安电子科技大学 一种基于改进型scdae的半监督调制方式分类模型的方法
CN108492298A (zh) * 2018-04-13 2018-09-04 西安电子科技大学 基于生成对抗网络的多光谱图像变化检测方法
CN108492298B (zh) * 2018-04-13 2021-11-09 西安电子科技大学 基于生成对抗网络的多光谱图像变化检测方法
CN109191418B (zh) * 2018-06-22 2021-10-01 西安电子科技大学 一种基于收缩自编码器特征学习的遥感图像变化检测方法
CN109191418A (zh) * 2018-06-22 2019-01-11 西安电子科技大学 一种基于收缩自编码器特征学习的遥感图像变化检测方法
CN109377480A (zh) * 2018-09-27 2019-02-22 中国电子科技集团公司第五十四研究所 基于深度学习的耕地变化检测方法
CN111046887A (zh) * 2018-10-15 2020-04-21 华北电力大学（保定） 一种带噪图像特征提取的方法
CN109753896A (zh) * 2018-12-20 2019-05-14 西安电子科技大学 一种基于共性自动编码器的无监督异源遥感图像变化检测方法
CN111931035A (zh) * 2019-05-13 2020-11-13 中国移动通信集团湖北有限公司 业务推荐方法、装置及设备
CN111931035B (zh) * 2019-05-13 2023-11-24 中国移动通信集团湖北有限公司 业务推荐方法、装置及设备
CN110766084B (zh) * 2019-10-28 2021-06-15 北京理工大学 基于cae与hl-cnn的小样本sar目标识别方法
CN110766084A (zh) * 2019-10-28 2020-02-07 北京理工大学 基于cae与hl-cnn的小样本sar目标识别方法
CN112836736A (zh) * 2021-01-28 2021-05-25 哈尔滨理工大学 一种基于深度自编码器构图的高光谱图像半监督分类方法
CN112836736B (zh) * 2021-01-28 2022-12-30 哈尔滨理工大学 一种基于深度自编码器构图的高光谱图像半监督分类方法

Also Published As

Publication number Publication date
CN107392940B (zh) 2020-10-30

Similar Documents

Publication Publication Date Title
CN107392940A (zh) 一种基于堆栈半监督自适应去噪自编码器的sar图像变化检测方法
Giacinto et al. Design of effective neural network ensembles for image classification purposes
Dos Santos et al. A dynamic overproduce-and-choose strategy for the selection of classifier ensembles
CN108805188A (zh) 一种基于特征重标定生成对抗网络的图像分类方法
Balakrishnan et al. A study of the classification capabilities of neural networks using unsupervised learning: A comparison with K-means clustering
Song et al. A competitive pixel-object approach for land cover classification
CN104484681B (zh) 基于空间信息和集成学习的高光谱遥感影像分类方法
Malone et al. Data mining using rule extraction from Kohonen self-organising maps
Ghosh et al. Semi-supervised change detection using modified self-organizing feature map neural network
Chang et al. Change detection of land use and land cover in an urban region with SPOT-5 images and partial Lanczos extreme learning machine
Stavrakoudis et al. A genetic fuzzy-rule-based classifier for land cover classification from hyperspectral imagery
Stavrakoudis et al. A boosted genetic fuzzy classifier for land cover classification of remote sensing imagery
CN103942562B (zh) 基于多分类器组合的高光谱图像分类方法
CN104268570B (zh) 一种基于类内差异的层次化单分类舰船目标虚警剔除方法
Bragilevsky et al. Deep learning for Amazon satellite image analysis
Vatsavai High-resolution urban image classification using extended features
CN110837787B (zh) 一种三方生成对抗网络的多光谱遥感图像检测方法及系统
Stavrakoudis et al. A multistage genetic fuzzy classifier for land cover classification from satellite imagery
Abazar et al. An efficient ensemble of convolutional deep steganalysis based on clustering
Eeti et al. A modified class-specific weighted soft voting for bagging ensemble
Otero et al. Intercomparison of deep learning architectures for the prediction of precipitation fields with a focus on extremes
Kalita et al. Class-Wise Subspace Alignment-Based Unsupervised Adaptive Land Cover Classification in Scene-Level Using Deep Siamese Network
Vassilas Efficient neural network-based methodology for the design of multiple classifiers
Canuto et al. The use of confidence measures to enhance combination strategies in multi-network neuro-fuzzy systems
Pajares et al. Combining classifiers through fuzzy cognitive maps in natural images

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