具体实施方式

下面将结合实施例对本发明的实施方案进行详细描述，但是本领域的技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围。

参考图1，本发明提供的一种基于自适应参数遗传算法的森林变化遥感检测方法，包括以下步骤：

步骤1，获取不同时刻的两幅森林遥感影像I₁和I₂，对每幅所述森林遥感影像分别进行预处理，得到两幅预处理后的森林遥感影像I₁”和I₂”；图像大小为M×N×B；M、N分别表示I₁和I₂的行数和列数，B表示图像的光谱波段数目，其中，M≥20,N≥20,B≥4；

步骤1.1，读取两个不同时刻的哨兵2A卫星森林遥感影像，其中M＝350，N＝200，B＝4；

步骤1.2，由于遥感影像数据在获取和传输过程中易产生系统的、随机的辐射失真或畸变，从而造成遥感图像的失真，影响遥感图像的判读和解译。为消除或改正因辐射误差而引起影像畸变，本发明将森林遥感图像I₁和I₂输入遥感影像处理软件ENVI 5.2中，使用ENVI的radiometric calibration工具包对I₁和I₂进行辐射校正，获取辐射校正后的森林遥感影像I₁'、I₂‘；

由于大气和光照等因素会影响地物反射率，为了获取地物真实的反射率，本发明利用ENVI 5.2中的FLASH Atmospheric Correction工具包对辐射校正后的森林遥感影像I₁'、I₂‘进行大气校正,获得大气校正后的森林遥感影像I₁”、I₂”，即为预处理后的图像。

步骤2，利用光谱波段的反射率分别计算两幅预处理后的森林遥感影像I₁”、I₂”的归一化植被指数矩阵NDVI₁、NDVI₂，并令其为对应的植被指数特征图SF₁、SF₂；

植被指数可以根据植物的光谱特性，将卫星可见光和近红外波段进行组合,实现对地表植被状况的有效度量，有效区分植被与土壤，在森林变化的检测领域被广泛应用。归一化植被指数NDVI是植被覆盖的最佳指示因子，对土壤背景的变化较敏感，在很大程度上消除地形和阴影的影响，对植被提取效果较好。本发明利用遥感影像光谱波段的反射率计算大气校正后的森林遥感影像的归一化植被指数矩阵NDVI，计算公式为：

其中，ρ_NIR和ρ_RED分别表示森林遥感影像的红外光谱波段反射率和红光谱波段反射率，令每幅预处理后的森林遥感影像I₁”、I₂”求得的NDVI矩阵为对应森林影像的植被指数特征图SF₁、SF₂，图像大小为M×N，具体为350×200。

步骤3，对植被指数特征图SF₁和SF₂分别进行中值滤波，对应得到去噪后的特征图SF₁'、SF₂'；使用差值法对去噪后的特征图SF₁'、SF₂'构造植被指数差异图DI；

3.1，使用中值滤波算法对两幅植被指数特征图SF₁、SF₂进行去噪处理，得到去噪后的植被指数特征图SF₁'、SF₂’；

3.2，使用差值法对去噪后的特征图SF₁'、SF₂'构造植被指数差异图DI，即将两幅植被指数特征图SF₁、SF₂对应位置的像素值做差，将差值作为DI(x,y)对应位置的像素值，遍历所有像素点即可得到差异图DI。

DI(x,y)＝SF₁'(x,y)-SF₂'(x,y)。

步骤4，对植被指数差异图DI进行显著性检测，得到显著性差异图DS；

4.1，为了减少分类计算的搜索空间，本发明对植被指数差异图DI进行了一种基于全局对比度的自下而上的图像显著性检测。该方法计算每个像素在整个图像上的全局对比度，把该像素与图像中其他所有像素在颜色上的欧式距离之和作为该像素的显著性特征，计算公式为：

其中，Sam(x,y)是像素(x,y)的显著性特征，且(x,y)≠(i,j)；然后将Sam(x,y)归一化为[0,255]。因为植被指数差异图DI是单通道图像，所以显著性特征值就是像素的灰度值。

4.2，对显著性图Sam使用OTSU(最大类间方差法)图像分割算法，生成二值化的显著性差异图S，其中非显著区域的像素值为0，显著性区域的像素值为1；然后利用二值化的显著性图S和植被指数差异图DI生成显著性差异图DS：

DS(x,y)＝S(x,y)·DI(x,y)。

步骤5，使用模糊C均值聚类算法对显著性差异图DS进行预分类，将DS中的像素划分为变化类、未变化类和不确定类，对应得到变化类像素集、未变化类像素集和不确定类像素集；

5.1，设置模糊C均值聚类算法的目标函数为：

其中，c为类别数，n为显著性差异图DS的像素总数，u_ij代表像素x_j属于第i聚类的隶属度，取值为[0,1]，且m表示模糊权重系数，d_ij＝||x_j-v_i||，表示像素x_j与聚类中心v_i的欧氏距离。

具体操作中，本发明设c＝3，m＝2，阈值ε＝10^-5，最大迭代次数T＝10³；并初始化迭代次数t＝1，随机初始化三个聚类中心v₁、v₂、v₃；

5.2，计算新的聚类中心和隶属度函数：

若当前迭代代数t＝1，则直接根据聚类中心的像素值更新每个像素的隶属度；若当前迭代代数t≥2，则先用上一代求得的隶属度矩阵更新聚类中心，再根据聚类中心更新隶属度矩阵。具体公式为：

其中，k＝1,2,3；表示第t次迭代中第i聚类中心，表示第t次迭代中像素x_j属于第i聚类的隶属度；

5.3，若且t≤T，则令t＝t+1，返回步骤5.2；否则，迭代终止，根据每个像素最终的隶属度分配其类别标签：对于像素x_j，若arg_k{max{u_kj}}＝arg_k{max{v_k}}，且max{u_kj}≥0.9，则该像素属于变化类；arg_k{max{v_k}}表示具有最大聚类中心像素值的类别k，arg_k{max{u_kj}}表示像素x_j的最大隶属度所对应的类别；若arg_k{max{u_kj}}＝arg_k{min{v_k}}，且max{u_kj}≥0.9，则该像素属于未变化类；其他情况下，该像素属于不确定类。在步骤6之后的算法中不再对已经确定类别的像素即变化类和未变化类的像素进行计算。

步骤6，使用基于自适应参数的遗传算法对不确定类像素集中像素的类别标签进行优化，得到不确定类对应的优化后的类别标签；根据优化后的类别标签生成最终的变化图。

6.1，种群初始化：

设置种群大小为Num＝40，即种群中有40个编码为二进制二维矩阵的个体，个体矩阵Ind的大小与显著性差异图DS大小相同，即350×200；根据步骤5所得的类别标签，对所有个体矩阵Ind进行初始化，即将变化类像素集对应位置的像素值修改为1，未变化类像素集对应位置的像素值修改为0，不确定类像素集对应位置的像素值随机编码为0或1；并建立不确定类像素在个体矩阵中的位置索引，设迭代次数k'＝1，阈值ε'＝10^-2，最大迭代次数K＝10⁵；

6.2，计算个体适应度：

根据个体与目标函数的适应度值确定选择该个体执行遗传操作的概率。与传统遗传算法不同，本发明所使用的改进遗传算法引入了像素的空间信息，设计带有邻域因子的目标函数OF，计算公式为：

其中，r为类别标识，r＝0,1，r＝0对应未变化类，r＝1对应变化类；c_r表示显著性区域中对应类别的所有像素，n_r和m_r表示对应类别像素的统计个数和平均灰度值，N_q表示中心像素x_j的3×3邻域N_j中的第q个像素，其中，q＝1,2,3...,9，l_q为N_q与中心像素x_j的欧氏距离，L是像素x_j与其周围像素之间的距离权重；DS_j表示显著性差异图DS中与中心像素x_j位置对应的像素。

在得到每个个体的目标函数后，对其取倒数求得每个个体的适应度Fitness：

6.3，对种群进行精英选择和交叉操作：

若当前迭代代数k'＝1，则对所有个体进行遗传操作；若当前迭代代数k'≥2，则当前迭代的个体为子代个体，k'-1代个体为父代个体，那么，根据子代和父代的所有个体的适应度值从高到低选择Num个个体执行遗传操作，即选择40个个体，挑选出适应度值最高的个体作为当前迭代代数的最佳个体，并将最佳适应度值保存在矩阵Maxfit中，未被选中的个体在进化中被淘汰。然后对被选中执行遗传操作的个体进行交叉操作，具体做法为：首先，将执行遗传操作的所有个体进行两两分组，形成多组待交叉个体；再随机在(0,1)之间初始化一个与个体大小一致的矩阵A，遍历矩阵A的元素，当某元素数值小于预先设置的交叉概率p_c＝0.8，则标记该元素的位置；然后将每组待交叉个体中与标记位置相同的元素进行交换，其他位置不变，形成交叉后的种群；

6.4，对交叉后的种群中的每个个体逐像素计算其变异概率：

对于交叉后的种群中的每个个体，遍历个体矩阵，对目标像素x_j的3×3邻域L_j中的每个邻域像素L_j(s)，s＝1,2,3...,9，计算其属于类别r的隶属度u_sr：

其中，r＝0,1：m_r表示步骤6.2)求得的各类别像素的平均灰度值；

根据隶属度u_sr，为邻域像素L_j(s)分配类别标签，若u_s0≥u_s1，则c_s＝0，否则，c_s＝1；因此，目标像素x_j的变异概率为：

其中，d_s是邻域像素L_j(s)到中心像素点x_j的空间距离，B_j为x_j的二进制编码值；

若求得的变异概率p(x_j)大于预先设置的标准变异概率p＝0.2，则该位置的编码改变，即0变为1或者1变为0，否则编码值不变；遍历完每个个体矩阵则得到变异后的种群，并将其作为上一代的子代种群；计算子代种群中每个个体的适应度，并选出适应度最大的作为该子代的最佳适应度；

6.5，若相邻两次迭代代数求得的最佳适应度的差大于所设定的阈值，即|Maxfit_k'+1-Maxfit_k'|≥ε'，且k'≤K，则令k'＝k'+1，返回步骤6.3，否则输出本次迭代的适应度最高的个体；

6.6，在输出的最佳个体上补充步骤5中得到的变化类和未变化类对应的像素的类别标签，输出最终的带有标签的二进制变化图。

仿真实验

下面通过仿真数据处理结果进一步说明本发明的正确性和有效性。

仿真内容：使用原始遗传算法与本发明算法对图2(a)和图2(b)所示的两幅多光谱森林图像进行变化检测，结果如图2(c)和图2(d)所示，其中，图2(c)是原始遗传算法所得的变化图，图2(d)是本发明算法求得的变化图。

最后，统计其平均分类精度，结果如表1所示。

表1原始遗传算法与本发明算法的检测精度

从图2(c)和2(d)可以看出，本发明方法能够对遥感图像实现准确的地物分类，对于森林地物中的干扰因素能够有限排除，具有较强的鲁棒性，实现准确检测。且从表1可以看出，本发明具有更高的检测精度。

本发明通过使用森林遥感影像的NDVI光谱特征图作为变化检测待分图像，对其进行去噪和显著性检测预处理。为了进一步减少变化检测算法的搜索空间，本发明使用了模糊C均值聚类，根据灰度值给隶属度大于阈值的像素分配类别标签，对隶属度小于阈值的像素进行基于自适应参数的遗传算法分类。

本发明与传统遗传算法不同的是，在目标函数和变异概率中引入图像的空间邻域信息，自适应地修改遗传算子参数，避免了传统遗传算法分类模型性抗噪性差和搜索空间巨大的问题，有效地平衡了噪声抑制和细节保留，使进化朝着最优解的方向快速收敛，大大提高了算法的检测精度，从而更准确地获取森林地区的变化信息，增强了对林地变化的检测能力。

虽然，本说明书中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。