一种基于随机森林的卫星遥感影像云量计算方法
【技术领域】
[0001] 本发明属于卫星遥感影像质量检查技术领域，具体涉及一种基于随机森林的卫星 遥感影像云量计算方法。
【背景技术】
[0002] 在卫星遥感影像中，云层区域的存在将对影像自身质量和后续信息处理带来极大 的不利影响，因此云量的检测与识别是卫星遥感影像应用领域的主要问题之一。遥感云检 测技术可以用于删除卫星遥感影像中云所在区域的数据，大幅减少数据量，避免云量过大 的无效数据占用系统的存储空间、处理能力和传输带宽，它有星上在轨应用和地面应用两 种应用形式。
[0003] 现有的云检测方法主要是有三种，分别是基于光谱阈值的方法、基于影像特征的 方法以及综合方法。基于光谱阈值的方法根据云层自身的反射特性和温度特性，利用云 在不同波段光谱下的反射率，人工设定光谱阈值来检测，但实际云层区域由于季节、大气环 境、地理位置等因素带来的不稳定性使得该方法复杂度过高、适应性不强；基于影像特征的 方法提取影像的灰度、频率、纹理等特征，通过云图所含特征分类进行云检测，但由于云和 地物在一些特征方面存在重叠现象，检测结果取决于所选取特征的有效性、弱相关性、完整 性；综合方法利用光谱阈值方法进行初检，筛选出备选云区域，再用特征提取的方法对这些 区域再次进行云检测。
[0004] 目前的云检测方法存在以下问题：第一、在轨云检测方法对设备体积、重量、功耗 有较大的约束，限制了算法的复杂度及适应性，不能保证较理想的云检测效果；第二、目前 已有的阈值法只针对某一颗卫星，缺乏一般系统方法，通用性差，且检测结果受到时空类型 的影响，可靠性不高；第三、目前已有的基于影像特征的方法，只应用了单独的纹理或亮度 或频率特征，在选取特征的完整性上有许多欠缺，导致检测方法的适应性不强，对厚云的检 测效果尚好，但对于薄云、低云的检测仍然存在难度。
【发明内容】
[0005] 为了解决上述技术问题，本发明提供了一种可以提高云检测精度的同时提升云检 测方法的实用性和通用性，使其可以应用到资源三号、天绘一号以及高分一号等国产卫星 影像产品质量检查系统中的卫星遥感影像云量计算方法。
[0006] 本发明所采用的技术方案是：一种基于随机森林的卫星遥感影像云量计算方法， 其特征在于，包括以下步骤：
[0007] 步骤1 :样本获取；
[0008] 收集不同类型的遥感云图影像以及不同类型的地物影像，切分遥感云图影像，获 得云、地物影像样本，将云影像样本和地物影像样本作为训练集；
[0009] 步骤2 :特征提取；
[0010] 计算所有样本影像的灰度特征、频率特征和纹理特征矢量值，形成特征矢量集 合；
[0011] 步骤3 :影像分类器训练；
[0012] 使用随机森林方法来训练样本影像的特征矢量集合，得到由决策树森林构成的影 像分类器；
[0013] 步骤4:待测影像切分；
[0014] 将待测卫星遥感影像的原始影像进行下采样以获取缩略图，对缩略图进行影像切 分得到子影像，计算所有子影像的灰度特征、频率特征和纹理特征矢量值；
[0015] 步骤5:影像分类；
[0016] 将单个子影像的特征矢量值输入影像分类器，影像分类器中的每个决策树都对该 特征矢量进行分类投票，最终按照它在"云类"和"无云类"得票的多少来判断对应的子影 像是否为含云区域；
[0017] 步骤6:云量计算；
[0018] 用步骤5所述的方法对所有子影像进行分类，分类完毕之后即可计算该卫星遥感 影像的云量百分比。
[0019] 作为优选，步骤1中所述的切分遥感云图影像，获得云、地物影像样本，是对卫星 遥感影像进行下采样获得缩略图，分别切分缩略图中的云图和无云图为32X32像素的样 本影像。
[0020] 作为优选，步骤2的具体实现包括以下子步骤：
[0021] 步骤2. 1 :计算样本影像的灰度特征；
[0022] 选用灰度均值、灰度方差、一阶差分、直方图信息熵作为灰度特征矢量；其具体实 现包括以下子步骤：
[0023] 步骤2. I. 1 :计算样本影像的灰度均值：
[0025] 其中，f (i, j)是影像（i, j)处的灰度值，M是影像的宽，N是影像的尚；
[0026] 步骤2. 1. 2 :计算样本影像的灰度方差：
[0028] 灰度方差反映了影像整体灰度的分布均匀程度；
[0029] 步骤2. 1. 3 :计算样本影像的一阶差分：
[0033] -阶差分表达了影像中灰度变化的剧烈程度；
[0034] 步骤2. 1. 4 :计算样本影像的直方图信息熵：
[0036] 其中，Hist[g]是影像的直方图，Hist[g] (i)是在某灰度级下的像素分布频率，直 方图信息熵综合反映了影像灰度分布以及影像的有序程度；
[0037] 步骤2. 2 :计算样本影像的频率特征；
[0038] 选择傅里叶变换高频系数以及小波变换高频系数作为频率特征矢量，具体实现包 括以下子步骤：
[0039] 步骤2. 2. 1 :计算样本影像的傅里叶变换高频系数，选择如下的傅里叶变换函数：
[0041] 当 u, V = 0 时，C(U)C(V) = 2 1/2;其他情况，C(U)C(V) = 1 ;
[0042] 步骤2. 2. 2 :计算样本影像的小波变换高频系数，使用多贝西小波中的哈尔小波 基函数对影像进行小波变换，其表达式如下：
[0044] 对应的尺度函数为：
[0046] 步骤2. 3 :计算样本影像的纹理特征；
[0047] 选择灰度梯度共生矩阵的二次统计量：梯度均方差、混合熵和逆差距，以及影像的 纹理分数维作为纹理特征矢量，具体实现包括以下子步骤：
[0048] 步骤2.3. 1 :首先计算出影像的灰度梯度共生矩阵H(i，j)，并对其进行归一化处 理，得到归一化后的灰度梯度共生矩阵,从而用于计算二次统计特征量；
[0049] 步骤2. 3. 2 :计算样本影像的梯度均方差，使用如下公式：
[0051] 其中，Tavg为梯度平均，其表达式如下：
[0053] Lg表示最大灰度级，Ls表示最大梯度值；表示归一化后的灰度梯度共生矩 阵；
[0054] 步骤2. 3. 3 :计算样本影像的混合熵，使用如下公式：
[0056] 步骤2. 3. 4 :计算样本影像的逆差矩，使用如下公式：
[0058] 步骤2. 3. 5 :计算样本影像的纹理分数维，使用分形布朗分维估计方法来求解影 像的纹理分数维值，该方法的数学描述如下：
[0059] 设X e Rn，f (X)是关于X的实随机函数，若存在常数H (0 < H < 1)，使得F⑴满
是一个与X，△ X无关的分布函数时，则f (X)称为分形布朗函 数；其中H称为自相似参数，则影像的分数维D的表达式为：
[0060] D = η+1-H 〇
[0061] 作为优选，步骤3中所述的随机森林方法，其具体实现包括以下子步骤：
[0062] 步骤3. 1 :将所有样本的特征矢量作为训练影像分类器的训练集S，特征维数为记 为F，使用到的决策树的数量记为t，每个节点使用到的特征数量记为f ;终止条件为：节点 上达到最少样本数sc或树达到最大深度d ;
[0063] 步骤3. 2 :从总的训练集S中有放回地多次随机采样，抽取大小和S -样的训练集 得到每个决策树的训练集S(i)，对于第x(x彡t)棵树，i = X ;以S(i)作为根节点的样本， 从根节点开始训练；
[0064] 步骤3. 3 :如果当前节点达到终止条件，则设置当前节点为叶子节点，该叶子节点 用于分类时的预测输出为当前节点样本集合中数量最多的那一类c(j)，其占当前节点总样 本的比例记为概率P (j);如果当前节点没有达到终止条件，则从F维特征中无放回地随机 选取f维特征；计算这f维特征各自的Gini系数，寻找最优的特征k和阈值th，之后将当 前节点上样本第k维特征小于th的样本被划分到左子节点，其余的被划分到右子节点；
[0065] Gini系数的计算公式为：
[0066] Gini = 1- Σ (p(j) · p(j));
[0067] 其中，p (j)为当前节点下第j类样本所占比例；
[0068] 判断标准的数学表达式为：
[0069] argmax (Gini-GiniLeft_GiniRight);
[0070] 其中，Gini是当前节点的Gini系数，Giniu3ft是左子节点的Gini系数，Gini [^^是 右子节点的Gini系数，argmax为取最大值；
[0071] 步骤3. 4 :重复执行步骤3. 2和步骤3. 3直到所有节点都训练过或者被标记为叶 子节点；
[0072] 步骤3. 5 :重复执行步骤3. 2、步骤3. 3、步骤3. 4直到所有决策树都被训练过。
[0073] 作为优选，步骤4中所述的计算所有子影像的灰度特征、频率特征和纹理特征矢 量值，与步骤2中所述的所有样本影像的灰度特征、频率特征和纹理特征矢量值的计算方 法相同。
[0074] 作为优选，步骤5中所述的将单个子影像的特征矢量值输入影像分类器，其具体 实现包括如下子步骤：
[0075] 步骤5.