基于网络结构优化的SAR图像变化区域检测方法

摘要

本发明公开了一种基于网络结构优化的合成孔径雷达SAR图像变化区域检测方法，主要解决现有技术无法提取合成孔径雷达SAR图像的局部特征，需要人工输入学习误差值进行测试的问题。本发明的具体步骤如下：(1)读入合成孔径雷达SAR图像；(2)归一化处理；(3)构建变化检测矩阵；(4)更新变化检测矩阵；(5)选择训练样本集；(6)搭建堆栈自动编码器网络；(7)训练堆栈自动编码器网络；(8)优化堆栈自动编码器网络；(9)构建概率矩阵；(10)更新概率矩阵；(11)获得变化检测图像。本发明提取了合成孔径雷达SAR图像的局部特征，提高了变化检测的精度。

权利要求

1.一种基于网络结构优化的合成孔径雷达SAR图像变化区域检测方法，其特征在于，该方法是构建一个堆栈自动编码器网络，将合成孔径雷达SAR图像的每一个像素点与其周围邻域的像素点输入堆栈自动编码器网络，使用粒子群算法，优化堆栈自动编码器网络的结构，得到图像的局部特征，用局部特征对所有像素点进行分类，得到最终的变化检测结果图，该方法的具体步骤包括如下：(1)读入合成孔径雷达SAR图像：

读入同一地区不同时相的两幅已配准且校正的合成孔径雷达SAR图像I1和I2；

(2)归一化处理：

使用归一化公式，对读入的合成孔径雷达SAR图像I1和I2进行归一化处理，得到归一化后的合成孔径雷达SAR图像I1'和I2'；

(3)构建变化检测矩阵：

(3a)将归一化后的合成孔径雷达SAR图像I1'的每一个像素点与其周围5×5邻域的25个像素点，作为检测子矩阵M1的行向量；

(3b)将归一化后的合成孔径雷达SAR图像I2'的每一个像素点与其周围5×5邻域的25个像素点，作为检测子矩阵M2的行向量；

(3c)构建一个变化检测矩阵M，该变化检测矩阵M的行数等于检测子矩阵M1的行数和检测子矩阵M2的行数之和，该变化检测矩阵M的列数等于检测子矩阵M1的列数，将该变化检测矩阵M的每一个元素的初始值设置为0；

(4)更新变化检测矩阵：

将检测子矩阵M1与检测子矩阵M2的所有行向量，依次存放到变化检测矩阵M的行向量中，组成更新后的变化检测矩阵；

(5)选择训练样本集：

(5a)从更新后的变化检测矩阵中随机选择一个行向量，作为一个训练样本；

(5b)从更新后的变化检测矩阵中随机选择行向量总数的75％的行向量，组成一个训练样本集；

(6)搭建堆栈自动编码器网络：

(6a)搭建一个5层堆栈自动编码器网络，其结构为：输入层→第一个隐藏层→第二个隐藏层→第三个隐藏层→输出层，其中，堆栈自动编码器网络的每一个隐藏层的结构为一个自动编码器；

(6b)设置堆栈自动编码器网络的每一层参数如下：将输入层的节点总数设置为50，将第一个隐藏层的节点总数设置100，将第二个隐藏层的节点总数设置为50，将第三个隐藏层的节点总数设置为25，将输出层的节点数设置为1；

(7)训练堆栈自动编码器网络：

将训练样本集的样本，依次输入到堆栈自动编码器网络中，逐层训练堆栈自动编码器网络；

(8)使用粒子群算法，逐层优化训练好的堆栈自动编码器网络：(8a)将每一个权重粒子变量的初始值设置为0；

(8b)用每一个自动编码器的每一个权重值w1，更新每一个权重粒子变量的值；

(8c)将每一个偏差粒子变量的初始值设置为0；

(8d)用每一个自动编码器的每一个偏差值b1，更新每一个偏差粒子变量的值；

(8e)按照下式，计算一个权重粒子变量的适应度最小值：p1＝min{f11,f12,...,f1N}

其中，p1表示一个权重粒子变量的适应度最小值，min表示取最小值操作，f11表示第1个权重粒子值适应度，f12表示第2个权重粒子值的适应度，N表示一个自动编码器权重值的总数，f1N表示第N个权重粒子值的适应度；

(8f)按照下式，计算一个偏差粒子变量的适应度最小值：p2＝min{f21,f22,...,f2M}

其中，p2表示一个偏差粒子变量的适应度最小值，f21表示第1个偏差粒子值的适应度，f22表示第2个偏差粒子值的适应度，M表示一个自动编码器偏差值的总数，f2M表示第M个偏差粒子值的适应度；

(8g)将一个优化代数变量的初始值设置为0；

(8h)将每一个权重粒子值的更新前速度的初始值设置为0；

(8i)将每一个权重粒子值的更新后速度的初始值设置为0；

(8j)将每一个权重粒子值更新前的初始值设置为0；

(8k)将每一个权重粒子值的更新后的初始值设置为0；

(8l)将每一个偏差粒子值的更新前速度的初始值设置为0；

(8m)将每一个偏差粒子值的更新后速度的初始值设置为0；

(8n)将每一个偏差粒子值更新前的初始值设置为0；

(8o)将每一个偏差粒子值的更新后的初始值设置为0；

(8p)按照下式，计算每一个权重粒子值的更新后速度值：v2i＝v1i+c1×r1×(p1-w1i)+c2×r2×w1i,其中，v2i表示第i个权重粒子值的更新后速度值，i∈{1,2,...,N}，∈表示属于符号，{}表示集合符号，v1i表示第i个权重粒子值的更新前速度值，c1表示取值为2的局部学习速率值，r1表示一个介于0到1之间的随机数，w1i表示第i个权重粒子值更新前的值，c2表示取值为2的全局学习速率值，r2表示一个介于0到1之间的随机数；

(8q)将第j个权重粒子值更新前的值w1j加上第j个权重粒子值的更新后速度值v2j，j∈{1,2,...,N}，得到第j个更新值h2j；

(8r)利用第j个更新值h2j更新第j个权重粒子值w2j；

(8s)按照下式，计算每一个偏差粒子值的更新后速度值：v4t＝v3t+c3×r3×(p2-b1t)+c4×r4×b1t,其中，v4t表示第t个偏差粒子值的更新后速度值，t∈{1,2,...,M}，v3t表示第t个偏差粒子值更新前速度值，c3表示局部学习速率的值，r3表示一个介于0到1之间的随机数，b1t表示第t个偏差粒子值更新前的值，c4表示全局学习速率的值，r4表示一个介于0到1之间的随机数；

(8t)将第k个偏差粒子值更新前的值b1k加上第k个偏差粒子值的更新后速度值v4k，k∈{1,2,...,M}，得到第k个更新值h2k；

(8u)利用第k个更新值h2k更新第k个权重粒子值b2k；

(9)构建概率矩阵：

构建一个概率矩阵，概率矩阵的行数等于变化检测矩阵的行数，概率矩阵的列数等于变化检测矩阵的列数，该概率矩阵的每一个元素的初始值是0；

(10)更新概率矩阵：

(10a)将更新后的变化检测矩阵的每一个元素，依次输入到优化好的堆栈自动编码器网络中，得到每一个元素的变化概率，将变化概率保存到概率矩阵中；

(10b)遍历概率矩阵中的每一个元素，若元素的值大于等于0.5，则执行步骤(10c)，否则，执行步骤(10d)；

(10c)将元素的值更新为255；

(10d)将元素的值更新为0；

(10e)概率矩阵中的每一个元素更新完后，得到更新后的概率矩阵；

(11)获得变化检测图像：

将更新后的概率矩阵的所有元素值作为像素灰度值，利用Matlab中的imshow函数，将像素灰度值转化为变化检测图像。

2.根据权利要求1所述的基于网络结构优化的合成孔径雷达SAR图像变化区域检测方法，其特征在于：步骤(2)中所述的归一化公式如下：其中，I1'表示合成孔径雷达SAR图像I1归一化后的合成孔径雷达SAR图像，min表示取最小值操作，max表示取最大值操作，I2'表示合成孔径雷达SAR图像I2归一化后的合成孔径雷达SAR图像。

3.根据权利要求1所述的基于网络结构优化的合成孔径雷达SAR图像变化区域检测方法，其特征在于：步骤(7)中所述的逐层训练堆栈自动编码器网络的具体步骤如下：第一步，用-5到5以内的浮点数随机初始化每一层自动编码器的每一个权重值和每一个偏置值；

第二步，用所选择的训练样本集的样本单独地无监督训练每一个自动编码器，得到每一个自动编码器的最优权重值和最优偏差值；

第三步，使用基于最小交叉熵的反向传播算法，对堆栈自动编码器网络进行微调，得到训练好的堆栈自动编码器网络。

说明书

基于网络结构优化的SAR图像变化区域检测方法

技术领域

[0001] 本发明属于图像处理技术领域，更进一步涉及遥感图像变化检测技术领域中的一种基于网络结构优化的合成孔径雷达SAR(Synthetic Aperture Radar,SAR)图像变化区域检测方法。本发明可用于提取同一地区不同时段的两幅合成孔径雷达SAR图像的邻域像素信息，并用堆栈自动编码器网络对提取的像素信息进行学习，得到最终的变化检测图。

背景技术

[0002] 目前合成孔径雷达SAR图像变化区域检测的分析步骤大致分为两步：(1)生成差异图。此步骤是初步区分2幅合成孔径雷达SAR图像中的未变化类和变化类，并为差异图分析提供基础。其典型方法包括：差值法、均值比值法和对数比值法。这些方法存在的缺点是对相干斑噪声敏感，检测变化精度不高；(2)分析差异图。该步骤是对步骤(1)获得的差异图进行分析，提取变化信息，最终得到两幅图像的变化区域和非变化区域。其典型方法包括阈值法、聚类法。这些方法存在的缺点是不能有效地提取图像的局部特征。

[0003] 西安电子科技大学在其拥有的专利技术“基于稀疏表示的合成孔径雷达SAR图像相干斑抑制方法”(专利申请号：201110346349.7，授权公告号：CN102346908B)中提出了一种基于稀疏表示的合成孔径雷达SAR图像变化区域检测方法。该方法的步骤是，首先对原始合成孔径雷达SAR图像进行对数变换。然后将对数变换后的图像分割为的重叠分块，将图像块的自身信息作为控制因子对其进行稀疏表示。然后应用近似KSVD(K-Singular Value Decomposition,KSVD)算法进行字典学习，需要人工输入学习误差值进行测试，得到自适应字典和更新后的稀疏表示系数。再利用自适应字典和更新后的稀疏表示系数得到图像W，对W进行指数变换得到图像R。最终对由原始合成孔径雷达SAR图像Y和指数变换后的图像R得到的差值图像V进行非线性各向异性扩散，得到最终变化检测图像。该方法存在的不足之处是，需要人工输入学习误差值进行测试，易造成图像的部分纹理信息丢失，增加后期变化检测的误检率。

[0004] Zheng等人在其发表的论文“Using combined difference image and k-means clustering for合成孔径雷达SAR image change detection”(IEEE Geoscience and Remote Sensing Letters,2014,11(3)：691-695)中提出了一种简单实用的基于差异图融合的合成孔径雷达SAR图像变化区域检测方法。该方法的步骤是，首先用差值算子和对数比值算子分别处理合成孔径雷达SAR图像单个像素的值，图像的所有像素处理完毕后，得到合成孔径雷达SAR图像的差值差异图和对数比值差异图。然后对上述差异图进行均值滤波和中值滤波，初步去除噪声干扰和野点，获得融合差异图。最终使用K-means算法分析融合差异图。该方法存在的不足之处是，只考虑合成孔径雷达SAR图像的单个像素点信息，无法提取合成孔径雷达SAR图像的局部特征，降低后期变化检测的准确度。

发明内容

[0005] 本发明的目的在于针对上述已有技术的缺点，提出了一种基于网络结构优化的合成孔径雷达SAR图像变化区域检测方法，以实现对合成孔径雷达SAR图像变化区域的准确检测。该方法从变化检测矩阵中随机选择不同的训练样本，用选择的训练样本逐层训练堆栈自动编码器网络，使用粒子群算法逐层优化堆栈自动编码器网络，提升了堆栈自动编码器网络的自学习能力，提高了变化检测矩阵中每一个元素被划分为变化类的概率的准确度。该方法思路简单明确，通过有效提取合成孔径雷达SAR图像的局部特征提高了变化检测的精度。

[0006] 本发明实现上述目的思路是：首先构建归一化合成孔径雷达SAR图像的变化检测矩阵，从变化检测矩阵中随机选择不同的训练样本，逐层训练一个堆栈自动编码器网络。然后从变化检测矩阵中随机选择不同的训练样本，逐层优化堆栈自动编码器网络。再把变化检测矩阵所有元素作为测试数据，依次输入到优化好的堆栈自动编码器网络里进行测试。最后将堆栈自动编码器网络的输出结果依次保存到概率矩阵中，更新概率矩阵，最终输出变化检测图像。

[0007] 本发明的具体步骤包括如下：

[0008] (1)读入合成孔径雷达SAR图像：

[0009] 读入同一地区不同时相的两幅已配准且校正的合成孔径雷达SAR图像I1和I2；

[0010] (2)归一化处理：

[0011] 使用归一化公式，对读入的合成孔径雷达SAR图像I1和I2进行归一化处理，得到归一化后的合成孔径雷达SAR图像I1'和I2'；

[0012] (3)构建变化检测矩阵：

[0013] (3a)将归一化后的合成孔径雷达SAR图像I1'的每一个像素点与其周围5×5邻域的25个像素点，作为检测子矩阵M1的行向量；

[0014] (3b)将归一化后的合成孔径雷达SAR图像I2'的每一个像素点与其周围5×5邻域的25个像素点，作为检测子矩阵M2的行向量；

[0015] (3c)构建一个变化检测矩阵M，该变化检测矩阵M的行数等于检测子矩阵M1的行数和检测子矩阵M2的行数之和，该变化检测矩阵M的列数等于检测子矩阵M1的列数，将该变化检测矩阵M的每一个元素的初始值设置为0；

[0016] (4)更新变化检测矩阵：

[0017] 将检测子矩阵M1与检测子矩阵M2的所有行向量，依次存放到变化检测矩阵M的行向量中，组成更新后的变化检测矩阵；

[0018] (5)选择训练样本集：

[0019] (5a)从更新后的变化检测矩阵中随机选择一个行向量，作为一个训练样本；

[0020] (5b)从更新后的变化检测矩阵中随机选择行向量总数的75％的行向量，组成一个训练样本集；

[0021] (6)搭建堆栈自动编码器网络：

[0022] (6a)搭建一个5层堆栈自动编码器网络，其结构为：输入层→第一个隐藏层→第二个隐藏层→第三个隐藏层→输出层，其中，堆栈自动编码器网络的每一个隐藏层的结构为一个自动编码器；

[0023] (6b)设置堆栈自动编码器网络的每一层参数如下：将输入层的节点总数设置为50，将第一个隐藏层的节点总数设置100，将第二个隐藏层的节点总数设置为50，将第三个隐藏层的节点总数设置为25，将输出层的节点数设置为1；

[0024] (7)训练堆栈自动编码器网络：

[0025] 将训练样本集的样本，依次输入到堆栈自动编码器网络中，逐层训练堆栈自动编码器网络；

[0026] (8)使用粒子群算法，逐层优化训练好的堆栈自动编码器网络：

[0027] (8a)将每一个权重粒子变量的初始值设置为0；

[0028] (8b)用每一个自动编码器的每一个权重值w1，更新每一个权重粒子变量的值；

[0029] (8c)将每一个偏差粒子变量的初始值设置为0；

[0030] (8d)用每一个自动编码器的每一个偏差值b1，更新每一个偏差粒子变量的值；

[0031] (8e)按照下式，计算一个权重粒子变量的适应度最小值：

[0032] p1＝min{f11,f12,...,f1N}

[0033] 其中，p1表示一个权重粒子变量的适应度最小值，min表示取最小值操作，f11表示第1个权重粒子值适应度，f12表示第2个权重粒子值的适应度，N表示一个自动编码器权重值的总数，f1N表示第N个权重粒子值的适应度；

[0034] (8f)按照下式，计算一个偏差粒子变量的适应度最小值：

[0035] p2＝min{f21,f22,...,f2M}

[0036] 其中，p2表示一个偏差粒子变量的适应度最小值，f21表示第1个偏差粒子值的适应度，f22表示第2个偏差粒子值的适应度，M表示一个自动编码器偏差值的总数，f2M表示第M个偏差粒子值的适应度；

[0037] (8g)将一个优化代数变量的初始值设置为0；

[0038] (8h)将每一个权重粒子值的更新前速度的初始值设置为0；

[0039] (8i)将每一个权重粒子值的更新后速度的初始值设置为0；

[0040] (8j)将每一个权重粒子值更新前的初始值设置为0；

[0041] (8k)将每一个权重粒子值的更新后的初始值设置为0；

[0042] (8l)将每一个偏差粒子值的更新前速度的初始值设置为0；

[0043] (8m)将每一个偏差粒子值的更新后速度的初始值设置为0；

[0044] (8n)将每一个偏差粒子值更新前的初始值设置为0；

[0045] (8o)将每一个偏差粒子值的更新后的初始值设置为0；

[0046] (8p)按照下式，计算每一个权重粒子值的更新后速度值：

[0047] v2i＝v1i+c1×r1×(p1-w1i)+c2×r2×w1i,

[0048] 其中，v2i表示第i个权重粒子值的更新后速度值，i∈{1,2,...,N}，∈表示属于符号，{}表示集合符号，v1i表示第i个权重粒子值的更新前速度值，c1表示取值为2的局部学习速率值，r1表示一个介于0到1之间的随机数，w1i表示第i个权重粒子值更新前的值，c2表示取值为2的全局学习速率值，r2表示一个介于0到1之间的随机数；

[0049] (8q)将第j个权重粒子值更新前的值w1j加上第j个权重粒子值的更新后速度值v2j，j∈{1,2,...,N}，得到第j个更新值h2j；

[0050] (8r)利用第j个更新值h2j更新第j个权重粒子值w2j；

[0051] (8s)按照下式，计算每一个偏差粒子值的更新后速度值：

[0052] v4t＝v3t+c3×r3×(p2-b1t)+c4×r4×b1t,

[0053] 其中，v4t表示第t个偏差粒子值的更新后速度值，t∈{1,2,...,M}，v3t表示第t个偏差粒子值更新前速度值，c3表示取值为2的局部学习速率的值，r3表示一个介于0到1之间的随机数，b1t表示第t个偏差粒子值更新前的值，c4表示取值为2的全局学习速率的值，r4表示一个介于0到1之间的随机数；

[0054] (8t)将第k个偏差粒子值更新前的值b1k加上第k个偏差粒子值的更新后速度值v4k，k∈{1,2,...,M}，得到第k个更新值h2k；

[0055] (8u)利用第k个更新值h2k更新第k个权重粒子值b2k；

[0056] (9)构建概率矩阵：

[0057] 构建一个概率矩阵，概率矩阵的行数等于变化检测矩阵的行数，概率矩阵的列数等于变化检测矩阵的列数，该概率矩阵的每一个元素的初始值是0；

[0058] (10)更新概率矩阵：

[0059] (10a)将更新后的变化检测矩阵的每一个元素，依次输入到优化好的堆栈自动编码器网络中，得到每一个元素的变化概率，将变化概率保存到概率矩阵中；

[0060] (10b)遍历概率矩阵中的每一个元素，若元素的值大于等于0.5，则执行步骤(10c)，否则，执行步骤(10d)；

[0061] (10c)将元素的值更新为255；

[0062] (10d)将元素的值更新为0；

[0063] (10e)概率矩阵中的每一个元素更新完后，得到更新后的概率矩阵；

[0064] (11)获得变化检测图像：

[0065] 将更新后的概率矩阵的所有元素值作为像素灰度值，利用Matlab中的imshow函数，将像素灰度值转化为变化检测图像。

[0066] 本发明与现有技术相比具有以下优点：

[0067] 第一，由于本发明搭建了一个堆栈自动编码器网络，该堆栈自动编码器网络充分得利用了合成孔径雷达SAR图像的每一个像素点的邻域像素信息，提取了合成孔径雷达SAR图像的局部特征，克服了现有技术中无法提取合成孔径雷达SAR图像的局部特征的问题，使得本发明在进行合成孔径雷达SAR图像变化检测时，提高了变化检测的精度。

[0068] 第二，由于本发明采用了粒子群算法，逐层优化堆栈自动编码器网络，自动地获得了堆栈自动编码器网络的最优权重值和最优偏差值，克服了现有技术中需要人工地输入学习误差值的问题，使得本发明在进行合成孔径雷达SAR图像变化检测时，提高了变化检测的精度。

附图说明

[0069] 图1是本发明的流程图；

[0070] 图2是本发明对Bern地区合成孔径雷达SAR图像的变化检测的仿真图；

[0071] 图3是本发明对Ottawa地区合成孔径雷达SAR图像的变化检测的仿真图；

[0072] 图4是本发明对Mulargia地区合成孔径雷达SAR图像的变化检测的仿真图。

具体实施方式

[0073] 下面结合附图对本发明作进一步描述。

[0074] 参照图1，本发明的具体实现步骤如下：

[0075] 步骤1，读入合成孔径雷达SAR图像。

[0076] 读入同一地区不同时相的两幅已配准且校正的合成孔径雷达SAR图像I1和I2。

[0077] 步骤2，归一化处理。

[0078] 使用归一化公式，对读入的合成孔径雷达SAR图像I1和I2进行归一化处理，得到归一化后的合成孔径雷达SAR图像I1'和I2':

[0079]

[0080]

[0081] 其中，I1'表示合成孔径雷达SAR图像I1归一化后的合成孔径雷达SAR图像，min表示取最小值操作，max表示取最大值操作，I2'表示合成孔径雷达SAR图像I2归一化后的合成孔径雷达SAR图像。

[0082] 步骤3，构建变化检测矩阵。

[0083] 将归一化后的合成孔径雷达SAR图像I1'的每一个像素点与其周围5×5邻域的25个像素点，作为检测子矩阵M1的行向量。

[0084] 将归一化后的合成孔径雷达SAR图像I2'的每一个像素点与其周围5×5邻域的25个像素点，作为检测子矩阵M2的行向量。

[0085] 构建一个变化检测矩阵M，该变化检测矩阵M的行数等于检测子矩阵M1的行数和检测子矩阵M2的行数之和，该变化检测矩阵M的列数等于检测子矩阵M1的列数，将该变化检测矩阵M的每一个元素的初始值设置为0。

[0086] 步骤4，更新变化检测矩阵。

[0087] 将检测子矩阵M1与检测子矩阵M2的所有行向量，依次存放到变化检测矩阵M的行向量中，组成更新后的变化检测矩阵。

[0088] 步骤5，选择训练样本集。

[0089] 从更新后的变化检测矩阵中随机选择一个行向量，作为一个训练样本。

[0090] 从更新后的变化检测矩阵中随机选择行向量总数的75％的行向量，组成一个训练样本集。

[0091] 步骤6，搭建堆栈自动编码器网络。

[0092] 搭建一个5层堆栈自动编码器网络，其结构为：输入层→第一个隐藏层→第二个隐藏层→第三个隐藏层→输出层，其中，堆栈自动编码器网络的每一个隐藏层的结构为一个自动编码器。

[0093] 设置堆栈自动编码器网络的每一层参数如下：将输入层的节点总数设置为50，将第一个隐藏层的节点总数设置100，将第二个隐藏层的节点总数设置为50，将第三个隐藏层的节点总数设置为25，将输出层的节点数设置为1。

[0094] 步骤7，训练堆栈自动编码器网络。

[0095] 将训练样本集的样本，依次输入到堆栈自动编码器网络中，逐层训练堆栈自动编码器网络。

[0096] 步骤8，使用粒子群算法，逐层优化训练好的堆栈自动编码器网络。

[0097] 将每一个权重粒子变量的初始值设置为0。

[0098] 用每一个自动编码器的每一个权重值w1，更新每一个权重粒子变量的值。

[0099] 将每一个偏差粒子变量的初始值设置为0。

[0100] 用每一个自动编码器的每一个偏差值b1，更新每一个偏差粒子变量的值。

[0101] 按照下式，计算一个权重粒子变量的适应度最小值：

[0102] p1＝min{f11,f12,...,f1N}

[0103] 其中，p1表示一个权重粒子变量的适应度最小值，min表示取最小值操作，f11表示第1个权重粒子值适应度，f12表示第2个权重粒子值的适应度，N表示一个自动编码器权重值的总数，f1N表示第N个权重粒子值的适应度。

[0104] 按照下式，计算一个偏差粒子变量的适应度最小值：

[0105] p2＝min{f21,f22,...,f2M}

[0106] 其中，p2表示一个偏差粒子变量的适应度最小值，f21表示第1个偏差粒子值的适应度，f22表示第2个偏差粒子值的适应度，M表示一个自动编码器偏差值的总数，f2M表示第M个偏差粒子值的适应度。

[0107] 将一个优化代数变量的初始值设置为0。

[0108] 将每一个权重粒子值的更新前速度的初始值设置为0。

[0109] 将每一个权重粒子值的更新后速度的初始值设置为0。

[0110] 将每一个权重粒子值更新前的初始值设置为0。

[0111] 将每一个权重粒子值的更新后的初始值设置为0。

[0112] 将每一个偏差粒子值的更新前速度的初始值设置为0。

[0113] 将每一个偏差粒子值的更新后速度的初始值设置为0。

[0114] 将每一个偏差粒子值更新前的初始值设置为0。

[0115] 将每一个偏差粒子值的更新后的初始值设置为0。

[0116] 按照下式，计算每一个权重粒子值的更新后速度值：

[0117] v2i＝v1i+c1×r1×(p1-w1i)+c2×r2×w1i,

[0118] 其中，v2i表示第i个权重粒子值的更新后速度值，i∈{1,2,...,N}，∈表示属于符号，{}表示集合符号，v1i表示第i个权重粒子值的更新前速度值，c1表示取值为2的局部学习速率值，r1表示一个介于0到1之间的随机数，w1i表示第i个权重粒子值更新前的值，c2表示取值为2的全局学习速率值，r2表示一个介于0到1之间的随机数。

[0119] 将第j个权重粒子值更新前的值w1j加上第j个权重粒子值的更新后速度值v2j，j∈{1,2,...,N}，得到第j个更新值h2j。

[0120] 利用第j个更新值h2j更新第j个权重粒子值w2j。

[0121] 按照下式，计算每一个偏差粒子值的更新后速度值：

[0122] v4t＝v3t+c3×r3×(p2-b1t)+c4×r4×b1t,

[0123] 其中，v4t表示第t个偏差粒子值的更新后速度值，t∈{1,2,...,M}，v3t表示第t个偏差粒子值更新前速度值，c3表示取值为2的局部学习速率的值，r3表示一个介于0到1之间的随机数，b1t表示第t个偏差粒子值更新前的值，c4表示取值为2的全局学习速率的值，r4表示一个介于0到1之间的随机数。

[0124] 将第k个偏差粒子值更新前的值b1k加上第k个偏差粒子值的更新后速度值v4k，k∈{1,2,...,M}，得到第k个更新值h2k。

[0125] 利用第k个更新值h2k更新第k个权重粒子值b2k。

[0126] 步骤9，构建概率矩阵。

[0127] 构建一个概率矩阵，概率矩阵的行数等于变化检测矩阵的行数，概率矩阵的列数等于变化检测矩阵的列数，该概率矩阵的每一个元素的初始值是0。

[0128] 步骤10，更新概率矩阵。

[0129] 将更新后的变化检测矩阵的每一个元素，依次输入到优化好的堆栈自动编码器网络中，得到每一个元素的变化概率，将变化概率保存到概率矩阵中。

[0130] 遍历概率矩阵中的每一个元素，若元素的值大于等于0.5，则执行第1步，否则，执行第2步。

[0131] 第1步，将元素的值更新为255。

[0132] 第2步，将元素的值更新为0。

[0133] 概率矩阵中的每一个元素更新完后，得到更新后的概率矩阵。

[0134] 步骤11，获得变化检测图像。

[0135] 将更新后的概率矩阵的所有元素值作为像素灰度值，利用Matlab中的imshow函数，将像素灰度值转化为变化检测图像。

[0136] 下面结合仿真实验对本发明的效果做进一步的说明。

[0137] 1、仿真条件：

[0138] 本发明的仿真实验是在主频2.30GHz的Intel Pentium(R)Dual-Core CPU、内存5GB的硬件环境和MATLAB R2016的软件环境下进行的。

[0139] 本发明仿真实验所使用的仿真参数如下：

[0140] 正检变化类像素数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中发生变化且实验结果图中检测为变化的像素个数，称为正检变化类像素数TP。

[0141] 正检未变化类像素数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中未发生变化且实验结果图中检测为未变化的像素个数，称为正检未变化类像素数TN。

[0142] 误检数：统计实验结果图中未发生变化区域的像素个数，与参考图中未变化区域的像素个数进行对比，把参考图中未发生变化但实验结果图中检测为变化的像素个数，称为误检数FP。

[0143] 漏检数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中发生变化但实验结果图中检测为未变化的像素个数，称为漏检数FN。

[0144] 总错误像素数＝误检数+漏检数。

[0145] 正确率PCC：PCC＝1-总错误像素数/总像素数。

[0146] 衡量检测结果图与参考图一致性的Kappa系数：

[0147]

[0148]

[0149] 其中，PRE表示理论检测精度的一致率,Nc表示实验结果图中像素值是255的像素的总个数，Nu表示实验结果图中像素值是0的像素的总个数，Mc为标准参考图像的像素点的总个数，卡帕系数Kappa表示变化检测的精度。

[0150] 2、仿真内容与结果分析：

[0151] 本发明的仿真实验共有3个。

[0152] (1)仿真实验1.

[0153] 本发明的仿真实验1是，使用本发明以及两个现有技术(基于邻域比值差异图NR(Neighborhood-based Ratio)方法、基于多目标优化模糊聚类MOFCM(Multi-objective Fuzzy Clustering Method)方法)，分别对Bern地区的合成孔径雷达SAR图像进行仿真，得到了Bern地区的合成孔径雷达SAR图像的变化检测仿真图。

[0154] 图2是仿真实验1对Bern地区合成孔径雷达SAR图像的变化检测仿真图。图2(a)是仿真实验1中使用的选自于Bern地区1999年4月的一幅合成孔径雷达SAR图。图2(b)是仿真实验1中使用的选自于Bern地区1999年5月的一幅合成孔径雷达SAR图。图2(c)是仿真实验1中使用的Bern地区的变化检测参考图。图2(d)是仿真实验1中采用NR方法，获得的Bern地区的合成孔径雷达SAR图像的变化检测仿真图。图2(e)是仿真实验1中采用MOFCM方法，获得的Bern地区的合成孔径雷达SAR图像的变化检测仿真图。图2(f)是采用本发明技术，获得的Bern地区的合成孔径雷达SAR图像的变化检测仿真图。

[0155] 图2(c)，图2(d)，图2(e)和图2(f)中的黑色部分表示仿真实验后检测到的变化区域，白色部分表示仿真实验后检测到的非变化区域。由图2(d)和图2(e)可见，两图中的黑色部分比图2(c)中的黑色部分少，两图中白色部分比图2(c)中的白色部分多，黑色部分少表示可以检测到的非变化区域小，非变化区域小表示变化检测精度低，白色部分多表示检测到的变化区域大，变化区域大表示变化检测精度低。由图2(f)可见，图中黑色部分比图2(c)中的黑色部分多，图2(f)中白色部分比图2(c)中的白色部分少，黑色部分多表示可以检测到的非变化区域大，非变化区域大表示变化检测精度低，白色部分少表示检测到的变化区域小，变化区域小表示变化检测精度高。

[0156] 对图2(d)进行参数分析，得到图2(d)的5个参数值：漏检数、误检数、总错误像素数、正确率、Kappa系数。对图2(e)进行参数分析，得到图2(e)的5个参数值：漏检数、误检数、总错误像素数、正确率、Kappa系数。对图2(f)进行参数分析，得到图2(f)的5个参数值：漏检数、误检数、总错误像素数、正确率、Kappa系数。将图2(d)、图2(e)及图2(f)对应的5个参数值列出表1。

[0157] 表1 Bern地区变化检测结果一览表

[0158] 方法 漏检数 误检数 总错误像素数 正确率 Kappa系数NR 414 83 497 0.9945 0.809
MOFCM 1762 868 2630 0.9290 0.764
本发明 196 247 443 0.9953 0.837

[0159] 由表1可见，本发明方法的总错误像素数比NR方法的总错误像素数少了54个，总错误像素数少表示变化检测精度高。本发明方法的总错误像素数比MOFCM方法的总错误像素数少了2178个，总错误像素数少表示变化检测精度高。本发明方法的Kappa系数比NR方法的Kappa系数大0.028，Kappa系数大表示变化检测精度高。本发明方法的Kappa系数比MOFCM方法的Kappa系数大0.073，Kappa系数大表示变化检测精度高。

[0160] 综上所述，采用NR方法和MOFCM方法误检测的非变化区域多，误检测的非变化区域多表示变化检测精度低。采用本发明的方法误检测的非变化区域少，误检测的非变化区域少表示变化检测精度高。

[0161] (2)仿真实验2.

[0162] 本发明的仿真实验2是，使用本发明以及两个现有技术(NR方法、MOFCM方法)，分别对Ottawa地区的合成孔径雷达SAR图像进行仿真，得到了Ottawa地区的合成孔径雷达SAR图像的变化检测仿真图。

[0163] 图3是仿真实验2对Ottawa地区合成孔径雷达SAR图像的变化检测仿真图。图3(a)是仿真实验2中使用的选自于Ottawa地区1997年5月的一幅合成孔径雷达SAR图。图3(b)是仿真实验2中使用的选自于Ottawa地区1997年8月的一幅合成孔径雷达SAR图。图3(c)是仿真实验2中使用的Ottawa地区的变化检测参考图。图3(d)是仿真实验2中采用NR方法，获得的Ottawa地区的合成孔径雷达SAR图像的变化检测仿真图。图3(e)是仿真实验2中采用MOFCM方法，获得的Ottawa地区的合成孔径雷达SAR图像的变化检测仿真图。图3(f)采用本发明技术，获得的Ottawa地区的合成孔径雷达SAR图像的变化检测仿真图。

[0164] 图3(c)，图3(d)，图3(e)和图3(f)中的黑色部分表示仿真实验后检测到的变化区域，白色部分表示仿真实验后检测到的非变化区域。由图3(d)和图3(e)可见，两图中的黑色部分比图3(c)中的黑色部分少，两图中白色部分比图3(c)中的白色部分多，黑色部分少表示可以检测到的非变化区域小，非变化区域小表示变化检测精度低，白色部分多表示检测到的变化区域大，变化区域大表示变化检测精度低。由图3(f)可见，图中黑色部分比图3(c)中的黑色部分多，图3(f)中白色部分比图3(c)中的白色部分少。黑色部分多表示可以检测到的非变化区域大，非变化区域大表示变化检测精度低，白色部分少表示检测到的变化区域小，变化区域小表示变化检测精度高。

[0165] 对图3(d)进行参数分析，得到图3(d)的5个参数值：漏检数、误检数、总错误像素数、正确率、Kappa系数。对图3(e)进行参数分析，得到图3(e)的5个参数值：漏检数、误检数、总错误像素数、正确率、Kappa系数。对图3(f)进行参数分析，得到图3(f)的5个参数值：漏检数、误检数、总错误像素数、正确率、Kappa系数。将图3(d)、图3(e)及图3(f)对应的5个参数值列出表2。

[0166] 表2 Ottawa地区变化检测结果一览表

[0167]方法 漏检数 错检数 总错误像素数 正确率 Kappa系数
NR 2473 1272 3745 0.9753 0.866
MOFCM 2521 3017 5538 0.9420 0.785
本发明 635 2032 2667 0.9731 0.894

[0168] 由表2可见，本发明方法的总错误像素数比NR方法的总错误像素数少了1078个，总错误像素数少表示变化检测精度高。本发明方法的总错误像素数比MOFCM方法的总错误像素数少了2871个，总错误像素数少表示变化检测精度高。本发明方法的Kappa系数比NR方法的Kappa系数大0.028，Kappa系数大表示变化检测精度高。本发明方法的Kappa系数比MOFCM方法的Kappa系数大0.109，Kappa系数大表示变化检测精度高。

[0169] 综上所述，采用NR方法和MOFCM方法误检测的非变化区域多，误检测的非变化区域多表示变化检测精度低。采用本发明的方法误检测的非变化区域少，误检测的非变化区域少表示变化检测精度高。

[0170] (3)仿真实验3.

[0171] 本发明的仿真实验3是，使用本发明以及两个现有技术(NR方法、MOFCM方法)，分别对意大利撒丁岛Mulargia湖泊区域的合成孔径雷达SAR图像进行仿真，得到了意大利撒丁岛Mulargia湖泊区域的合成孔径雷达SAR图像的变化检测仿真图。

[0172] 图4是仿真实验3对意大利撒丁岛Mulargia湖泊区域合成孔径雷达SAR图像的变化检测仿真图。图4(a)是仿真实验3中使用的选自于意大利撒丁岛Mulargia湖泊区域1996年7月的一幅合成孔径雷达SAR图。图4(b)是仿真实验3中使用的选自于意大利撒丁岛Mulargia湖泊区域的1996年9月的一幅合成孔径雷达SAR图像。图4(c)是仿真实验3中使用的意大利撒丁岛Mulargia湖泊区域的变化检测参考图。图4(d)是仿真实验3中采用NR方法，获得的意大利撒丁岛Mulargia湖泊区域的合成孔径雷达SAR图像的变化检测仿真图。图4(e)是本发明仿真实验3中采用MOFCM方法，获得的意大利撒丁岛Mulargia湖泊区域的合成孔径雷达SAR图像的变化检测仿真图。图4(f)是采用本发明技术，获得的意大利撒丁岛Mulargia湖泊区域的合成孔径雷达SAR图像的变化检测仿真图。

[0173] 图4(c)，图4(d)，图4(e)和图4(f)中的黑色部分表示仿真实验后检测到的变化区域，白色部分表示仿真实验后检测到的非变化区域。由图4(d)和图4(e)可见，两图中的黑色部分比图4(c)中的黑色部分少，两图中白色部分比图4(c)中的白色部分多，黑色部分少表示可以检测到的非变化区域小，非变化区域小表示变化检测精度低，白色部分多表示检测到的变化区域大，变化区域大表示变化检测精度低。由图4(f)可见，图中黑色部分比图4(c)中的黑色部分多，图4(f)中白色部分比图4(c)中的白色部分少。黑色部分多表示可以检测到的非变化区域大，非变化区域大表示变化检测精度低，白色部分少表示检测到的变化区域小，变化区域小表示变化检测精度高。

[0174] 对图4(d)进行参数分析，得到图4(d)的5个参数值：漏检数、误检数、总错误像素数、正确率、Kappa系数。对图4(e)进行参数分析，得到图4(e)的5个参数值：漏检数、误检数、总错误像素数、正确率、Kappa系数。对图4(f)进行参数分析，得到图4(f)的5个参数值：漏检数、误检数、总错误像素数、正确率、Kappa系数。将图4(d)、图4(e)及图4(f)对应的5个参数值列出表3。

[0175] 表3意大利撒丁岛Mulargia湖泊区域变化检测结果一览表

[0176]方法 漏检数 错检数 总错误像素数 正确率 Kappa系数
NR 6036 784 6820 0.9430 0.665
MOFCM 6541 3952 10493 0.9015 0.542
本发明 5615 247 5862 0.9522 0.679

[0177] 由表3可见，本发明方法的总错误像素数比NR方法的总错误像素数少了958个，总错误像素数少表示变化检测精度高。本发明方法的总错误像素数比MOFCM方法的总错误像素数少了4631个，总错误像素数少表示变化检测精度高。本发明方法的Kappa系数比NR方法的Kappa系数大0.014，Kappa系数大表示变化检测精度高。本发明方法的Kappa系数比MOFCM方法的Kappa系数大0.137，Kappa系数大表示变化检测精度高。

[0178] 综上所述，采用NR方法和MOFCM方法误检测的非变化区域多，误检测的非变化区域多表示变化检测精度低。采用本发明的方法误检测的非变化区域少，误检测的非变化区域少表示变化检测精度高。