本发明涉及遥感图像应用领域，更具体的说，涉及一种基于云分布特点的遥感图像地物完整性分析方法。
背景技术：
近年来，卫星遥感技术进入了一个能快速及时地获取对地观测海量信息的新阶段。卫星遥感图像被越来越多地应用在地球资源调查，自然灾害预测预报，环境污染监测，气象卫星云图处理，以及地面军事目标识别等各个方面。在实际情况下，由于天气情况复杂多变，遥感卫星获取卫星图像，尤其是获取大范围的遥感影像时，多数情况下都会有云的存在，真正的无云晴空地面目标物体影像很难获取。遥感图像上，云的存在会遮蔽或者覆盖云层区域下面的真实的地面目标物体的信息。为了能够有效的从遥感图像数据中提取地面目标物体的信息，提高遥感卫星图像数据的可用性和利用率，必须对遥感图像进行分析，通过分析确定其是否可用。
由于云的种类繁多，纹理呈现随机分布特性，其形态、特性等随时间和空间发生变化，同时还受季节和气候等多种条件的综合影响，即不同季节、不同空间高度、不同天气条件下的云的反射率和亮度都有所不同，这对分析遥感图像是否可用提出了新的挑战。
卫星遥感图像地物完整性分析的前提是假定所获的图像质量是符合标准的，虽然已有一些图像质量评估软件问世，将他们直接用于图像分析的效果仍不尽人意。卫星遥感图像地物完整性分析是根据图像用户使用需求，对图像中地物信息的完整性进行分析的过程。不管是国内还是国外，主观分析方法在使用时，对图像分析人员必须进行专门训练，既需要具有丰富的地学知识和目视判读经验，又需要花费大量的时间去目视判读，劳动强度大，信息获取周期长，遥感图像判读质量受目视判读者的个体差异、经验以及对解译区域的熟悉程度等各种因素制约。
在实际应用场合中，由于卫星拍摄的实时性，系统会获取海量的照片，其中很大部分图片因受到云分布的干扰而丧失实用价值，因此对获取的卫星遥感数据进行质量评估具有重要的意义。而目前对遥感图像的研究重点主要局限在云检测算法方面，对于图像使用价值的评估缺少理论一句，还依赖专业的图像分析人员进行目视判读，既耗费人力，降低效率，又无法保证判别标准的一致性。地物完整性判别方法的提出弥补了这一缺陷，为图像自身包含的信息量多少提供了科学的计算方法，提高了卫星遥感数据的处理效率，具有重要的理论意义和应用价值。
由于云的分布和下垫面地物的复杂性，导致目前云检测算法的性能、运算的自适应性、判别的灵活性、运行的实时性等多方面远不能满足可用度评估的要求。长期以来，云分布的定量统计一直以人工测试为主，受人为因素影响较大，不同研究者往往因其对云区域的结构、形态辨别能力不完全相同，对云分布的定量结果也有所差异。虽然国内外已经出现一些云检测算法，但是，这些方法的侧重点主要集中于云区域的检测，关于遥感图像地物完整性的分析探讨甚少。
技术实现要素：
为解决现有技术中存在的问题，本发明旨在提供一种基于云区域间拓扑结构与云区域内部信息熵的遥感图像地物信息完整性分析方法，在给出方法的理论定义的同时，提供详细的量化计算说明，以对卫星拍摄的遥感图像进行有依据的筛选甄别，筛选出地物完整性较高的可用图像，剔除地物完整性较低的无用图像，达到有效衡量遥感图像所包含的地物信息量多少的目的。
本发明所述地物完整性是指：在遥感图像中，地物受到不同程度的云遮挡后所保持的信息连续性。
为此，本发明提供如下技术方案：一种基于云分布特点的遥感图像地物完整性分析方法，其特征在于，包括如下步骤：
s1：计算云破碎度；
s2：计算信息损失抵抗度；
s3：根据步骤s1计算的云破碎度和步骤s2计算的信息损失抵抗度计算遥感图像的信息损失度；
s4：根据所述遥感图像的信息损失度分析遥感图像的地物完整性；
其中，云破碎度是指单位面积的云像素在图像空间被分割的程度，根据连通的云区域数量和大小来描述；
信息损失抵抗度是指被云覆盖的遥感图像对云覆盖信息损失的抵抗能力，其与可用面积比正相关，与对应地物类型的平均熵负相关；
遥感图像的信息损失度是指在云分布的影响下，遥感图像地物信息的完整度受到何种程度的损失；
遥感图像的信息损失度根据公式(1)确定：
其中，ild为信息损失度；cf为云破碎度，cf01为云破碎度归一化后的结果；dr为信息损失抵抗度。
优选地，所述步骤s1中云破碎度根据公式(2)确定：
其中，pd为云区域密度，指单位面积内云斑块连通域数量；sv为云区域面积方差。
优选地，对于同样面积的云像素，图像完整性被破坏程度随着它在图像中表现的斑块数量和斑块面积而变化，具体地说，
独立云团斑块数量越大，面积越小，则云遮挡像素分布得越离散，对图像中地物完整性破坏越大，反之，云遮挡像素分布得越密集连续，图像地物完整性保留程度增加；对于同样云遮挡面积，由于云分布不同，图像包含的地物信息多少不同；
所述云破碎度的计算公式中，考虑极限破碎的情况，定义图像尺寸为hsize×wsize，云区域数最多为[(hsize+1)/2]×[(wsize+1)/2]，其中每个云区域面积均为1，设图像i中云区域为ic，ic中云区域数量为rn(i)，sc为云覆盖总面积；pdmax为破碎到极限时单位面积内云斑块连通域数量；s总为图像总面积hsize×wsize；
将云破碎度经过归一化后可得：
优选地，所述i中云区域数量的统计方法包括如下步骤：
1)对检测出云的图像做区域分割，每个独立的云斑块属于一个独立区域，形成图像阵列a，每个区域中像素取相同的值，和其它区域的像素取值相异；
2)扫描图像阵列a并在每一个像素点[i,j]完成下列各步；
3)让r1＝a[i,j]；
4)查看像素在[i,j]的邻接像素[k,l]，对每一个邻接像素，进行下一步；
5)让r2＝a[k,l],如果r1≠r2，在图像的节点r1和r2之间增加一条弧线；
6)根据所述弧线和节点生成区域邻接图；
7)对于区域邻接图中的每一个云区域，完成下列步骤：
a)查看是否与邻接云区域相似；
b)合并相似的云区域，并修改生成的区域邻接图；
8)重复步骤7)，直到没有区域可以合并；
9)统计合并后的云区域数量。
优选地，所述步骤7)中查看是否与邻接云区域相似的步骤具体包括：
假定两个相邻云斑块区域r1和r2分别包含有m1和m2个像素，那么，有两种可能的假设：
h0：两个区域属于同一云斑块，此时两个区域的灰度值都服从单一高斯分布(μ0,σ02)；
h1：属于不同云斑块的区域，此时每一个区域的灰度值服从不同的高斯分布(μ1,σ12)和(μ2,σ22)；
对上面所述未知参数μ1,μ2,σ1,σ2使用样本估计，当区域包含有n个像素，每个像素灰度值为gi,i＝1,2,...,n,服从正态分布：
这些参数的最大似然估计方程为：
在假设h0下，所有的像素独立服从同一个分布n(u0,σ02)，在h0下的联合概率密度是：
在假设h1下，属于区域r1的m1个像素服从分布n(u1,σ12)，属于区域r2的m2个像素服从分布n(u2,σ22)，在h1下，联合分布密度函数为：
将这两种假设下的概率密度之比为似然比l：
计算得到似然比后分析：
如果似然比低于设定的阈值，则将这两个区域合并为一个区域来统计；
如果似然比高于设定阈值，则将这两个区域仍作为两个区域来统计；
在上述方程(4)-(9)中，参数μ1,μ2,σ1,σ2可由方程(5)和(6)通过使用m1+m2个像素估计得到。
应当注意，μ0,σ0只是作为泛化的参数说明所使用的概率公式原型，实际要求的参数只有上段所述的四个。
优选地，所述云区域面积方差根据公式(10)确定：
其中，si为第i个云区域面积，n为云区域数量。
优选地，所述信息损失抵抗度根据公式(11)确定：
dr＝aap/ae(11)
其中，aap为可用面积比；ae为遥感图像地物类型的平均熵；
所述可用面积比受基于地物类型的最小有效尺度影响，所述基于地物类型的最小有效尺度是指：针对某种地物类型的遥感图像，存在一个以像素为单位的最小的尺度，由它组成的矩形的区域刚好可以被某种识别算法完全识别出来，而小于这个尺度将无法识别的尺度；
所述可用面积比是指由全部由基于地物类型的最小有效尺度组成的面积占图像总面积的比例。
优选地，所述可用面积比的计算方法包括如下步骤：
step1：确认图像i的地物类型的mes；
step2：扫描i中未被云覆盖的区域，将其标记为g；
step3：扫描g中的每一个像素单元(i,j)，并执行以下步骤；
step4：将(i,j)标记为边长为mes的的矩形的一个顶点，并将这个矩形标记为mes(i,j)；
step5：若mes(i,j)∈g，存储mes(i,j)包含的所有像素，并返回step3；若
step6：处理完全部像素后，用所有存储的像素总数除以总像素数，计算得到aap的数值；
其中，mes是指上述基于地物类型的最小有效尺度。
优选地，所述遥感图像地物类型的平均熵是指遥感图像中灰度分布的聚集特征所包含的信息量，其根据公式(12)确定：
其中，pi是指图像中灰度值为i的像素所占的比例，即灰度值i的出现频率。
优选地，所述步骤s4具体包括如下内容：
分析计算得到的信息损失度，若信息损失度较高，则遥感图像地物完整性较差，不适合用于后续的陆地观测分析；
若信息损失度较低，则遥感图像地物完整性较好，适合用于后续的陆地观测分析。
与现有技术相对比，本发明产生的有益效果是：
(1)本发明提供的分析方法，为衡量遥感图像的地物信息量多少提供了描述依据，在云检测算法和图像的可用度评价标准之间架设了一座桥梁，使得遥感图像的分析系统能够更加完善。本发明提供的分析方法适用于各种类型的遥感图像分析，能够高效筛选出地物完整度强的遥感图像，可有效用于高精度地形测绘、地物提取与分类、森林防火、灾害评估、资源观测、和地质灾害灾后灾情分析等领域。
(2)本发明提供的分析方法中提出云破碎度的概念，通过对云破碎度进行分析，能够分析出云层的破碎程度对图像地物完整性保留程度的影响，云破碎度越高，说明云斑块越分散，图像地物完整性保留程度越低；反之，云破碎度越高，说明云斑块越大，分散程度越低，此时图像地物完整性保留程度就较高。
(3)本发明提供的分析方法中，提出了信息损失抵抗度的概念，通过对信息损失抵抗度进行分析，能够描述不同类型地物的遥感图像所具有的对信息破坏的抵抗能力，本发明研究了影响它大小的两方面因素，分别为信息保存因素(aap)及信息丢失因素(ae)，并将它们融合，能够较为准确地描述地物对信息破坏的抵抗能力，对遥感图像的地物完整性分析具有明显帮助。
(4)本发明创造性地提出遥感图像信息损失度的概念，通过云破碎度和地物对信息损失的抵抗度计算出信息损失度的值，能够分析遥感图像被云遮挡造成的损失有多大，进而模拟人类能够从图形中获取信息的感知能力，分析遥感图像的地物完整性。
附图说明
图1为五种典型地物样本图；
图2为云破碎度计算流程图；
图3为用于图像分割的区域邻接图；
图4为地物信息损失抵抗度计算流程图；
图5为地物完整性描述模型图。
具体实施方式
下面结合附图，对本发明的具体实施方式作详细的说明。
图1为五种典型地物样本图，在遥感图像中，一般将地物类型分为五种：海洋、沙漠、山脉、绿地和城市，各类样本图像如附图中中图1所示。分析过程中我们需要选用合适的特征可以对地物类型进行简单判别，本发明中利用分形维数作为地物类型的判别标准。
分形维数：反映了复杂形体占有空间的有效性，它是复杂形体不规则性的量度，代表图片的复杂程度。对均匀的、无缝隙的目标体，分形维数与拓扑维数相同，但是，分形维数是个分数值，而不是整数值。它的计算公式如下：
式中ε是立方体边长，n(ε)是用此立方体覆盖被测形体所得的数目，维数公式意味着通过用边长为ε的立方体覆盖被测形体来确定形体的维数。
从纹理上看，同种地物的图像之间具有某种一致性，而不同地物之间表现出明显不同的复杂性，例如海洋和沙漠的图像比较平缓，复杂性明显低于山脉、绿地和城市图像，分形维数作为一种特征方法正好以量化的方式反映出这一复杂性差异，因此本发明借用这一特征来区分不同地物类型。
图2为本发明提供的云破碎度计算流程图。本发明提出了一种基于云区域间拓扑结构的云破碎度描述方法，用以度量遥感图像中反映地物信息完整性的技术指标。该方法描述一幅遥感图像中云的空间分布特性。本发明中定义云破碎度(cloudfragment,cf)为单位面积的云像素在图像空间被分割的程度。根据连通的云区域数量和大小来反映。计算公式为：
其中，pd为云区域密度，指单位面积内云连通域数量；sv为云区域面积方差。
对于同样面积的云像素，图像完整性被破坏程度随着它在图像中表现的斑块数量和斑块面积而变化，具体来说，独立云团斑块数量越大，面积越小，则云遮挡像素分布得越离散，对图像中地物完整性破坏越大，反之，云遮挡像素分布得越密集连续，图像地物完整性保留程度增加。对于同样云遮挡面积，由于云分布不同，图像包含的地物信息多少不同。
云破碎度(cloudfragment,cf)计算公式(2)，考虑极限破碎的情况。图像尺寸为hsize×wsize，云区域数最多为[(hsize+1)/2]×[(wsize+1)/2]，其中每个云区域面积均为1。设图像中云区域为i，i中云区域数量为rn(i)，sc为云覆盖总面积；pdmax为破碎到极限时单位面积内云斑块连通域数量；s总为图像总面积hsize×wsize。经过归一化后可得：
为了准确计算归一化云破碎度cf01值，要解决两个主要问题。第一是准确统计云斑块连通域数量rn(i)(即云区域数量)，第二是准确计算云区域面积之间的方差sv。
本发明中提出基于云区域连通域拓扑关联方法进行云区域块数量统计，其中采用了区域邻接图生成算法和区域合并法，其中区域合并时考虑两个相邻云区域的统计特性，根据相邻区域是否具有相同的概率分布函数以便决定是否合并它们。
区域邻接图(regionadjacencygraphs,rag)表示图像中区域与区域之间的拓扑关系，它主要强调由区域构成图像的划分和每一个划分的特性。区域的不同特性可以存贮在不同的节点数据结构中。rag中的节点表示一个云区域，节点之间的弧线表示云区域之间的连接。
i中云区域数量的统计方法包括如下步骤：
1)对检测出云的图像做区域分割，每个独立的云斑块属于一个独立区域，形成图像阵列a，每个区域中像素取相同的值，和其它区域的像素取值相异；
2)扫描图像阵列a并在每一个像素点[i,j]完成下列各步；
3)让r1＝a[i,j]；
4)查看像素在[i,j]的邻接像素[k,l]，对每一个邻接像素，进行下一步；
5)让r2＝a[k,l],如果r1≠r2，在图像的节点r1和r2之间增加一条弧线；
6)根据所述弧线和节点生成区域邻接图，区域邻接图的样式见附图3，图3中(a)为步骤1)中已分割的图像，(b)为该图像的区域邻接图；
7)对于区域邻接图中的每一个云区域，完成下列步骤：
a)查看是否与邻接云区域相似；
b)合并相似的云区域，并修改生成的区域邻接图；
8)重复7)，直到没有区域可以合并；
9)统计合并后的云区域数量。
经过上面算法处理得到的已经标记的云连通域，可能存在很多相邻的区域，它们有可能属于同一片连续的云区域，当物理形态具有一定程度的相似性时，应被视为是一个云区域，否则会产生虚假的破碎度描述，故需进行合并。
上述步骤7)中查看是否与邻接云区域相似的步骤具体包括：
假定两个相邻云斑块区域r1和r2分别包含有m1和m2个像素，那么，有两种可能的假设：
h0：两个区域属于同一云斑块，此时两个区域的灰度值都服从单一高斯分布(μ0,σ02)。
h1：属于不同云斑块的区域，此时每一个区域的灰度值服从不同的高斯分布(μ1,σ12)和(μ2,σ22)。
对上面所述未知参数使用样本估计，当区域包含有n个像素，每个像素灰度值为gi,i＝1,2,...,n,服从正态分布：
这些参数的最大似然估计方程为：
在假设h0下，所有的像素独立服从同一个分布n(u0,σ02)。在h0下的联合概率密度是：
在假设h1下，属于区域r1的m1个像素服从分布n(u1,σ12)，属于区域r2的m2个像素服从分布n(u2,σ22)，在h1下，联合分布密度函数为：
将这两种假设下的概率密度之比为似然比l：
计算得到似然比后分析：
如果似然比低于设定的阈值，则将这两个区域合并为一个区域来统计；
如果似然比高于设定阈值，则将这两个区域仍作为两个区域来统计；
其中，上文所述阈值可以根据实际需要认为设定，因此本发明提供的方法具有通用性。
在上述方程(4)-(9)中，参数μ1,μ2,σ1,σ2可由方程(5)和(6)通过使用m1+m2个像素估计得到。
本发明提出独立连通区域的大小差异对破碎度的影响。采用云独立连通域的面积方差来衡量该差异大小，云区域之间面积差异越大，存在大面积云团的可能性越大；大块连续区域云对下垫面地物完整性破坏程度小于云区域数目多且面积小而碎的云区域对图像的破坏程度。
云区域面积方差根据公式(10)确定：
其中si为第i个云区域面积，n为云区域数量。
云区域之间面积差异越大，存在大面积云团的可能性越大，可以理解为，大块连续云对下垫面地物景象的破坏程度要小于斑块数目多且分布小而碎的云团对图像的破坏度。
图4为本发明提出的地物信息损失抵抗度计算流程图。
为了描述不同类型地物的遥感图像所具有的对信息破坏的抵抗能力，本发明研究了影响它大小的两方面因素，并分别称为信息保存因素及信息丢失因素，并将它们融合以描述地物对信息破坏的抵抗能力。
信息损失抵抗度根据公式(11)确定：
dr＝aap/ae(11)
其中，aap为可用面积比；ae为遥感图像地物类型的平均熵；
可用面积比：对于地物图像信息的保存因素，即图像中地物的可用面积，值得注意的是并非所有的地物区域块都能够有效地表示它们的信息，地物区域需要达到一定的尺寸才可能包含有效的、可用的、不被破坏的完整信息。由于不同类型地物的纹理特征不同，地物实际信息可分辨的区域尺寸与它们的类型强烈相关。为了描述这个尺寸的下限，并用于获取全部含有有效信息的地物区域，本发明提出了地物类型的最小有效尺度(minimaleffectivescale，mes)的概念：针对某种地物类型的遥感图像，存在一个最小的尺度(以像素为单位)，由它组成的矩形的区域刚好可以被某种识别算法完全识别出来，而小于这个尺度将无法识别，这个尺度称为最小有效尺度。使用mes来得出可用面积比(availableareapercentage，aap)，即由全部由mes组成的面积占图像总面积的比例。
可用面积比aap的计算方法包括如下步骤：
step1：确认图像i的地物类型的mes；
step2：扫描i中未被云覆盖的区域，将其标记为g；
step3：扫描g中的每一个像素单元(i,j)，并执行以下步骤；
step4：将(i,j)标记为边长为mes的的矩形的一个顶点，并将这个矩形标记为mes(i,j)；
step5：若mes(i,j)∈g，存储mes(i,j)包含的所有像素，并返回step3；若
step6：处理完全部像素后，用所有存储的像素总数除以总像素数，计算得到aap的数值；
图像的信息熵：图像熵表示为图像灰度级集合的比特平均数，单位比特/像素，也描述了图像信源的平均信息量。图像的信息熵表示遥感图像中灰度分布的聚集特征所包含的信息量，令pi表示图像中灰度值为i的像素所占的比例，即灰度值i的出现频率，则定义灰度图象的信息熵为：
其中，pi是指图像中灰度值为i的像素所占的比例，即灰度值i的出现频率。
图5为地物完整性描述模型图，本发明根据云破碎度对不同地物的影响因素，以及被云团覆盖遮挡的地物的自身性质，提出了遥感图像信息损失度(informationlossdegree，ild)概念用来测量cf和dr的相互作用。通过结合云破碎度的影响以及图像地物对本身的对这种影响的抵抗程度，ild值能够模拟人类从图像中获取信息的感知能力。
通过上文部分已经计算出cf01和dr，故此处即可通过公式(1)计算出图像信息损失度ild。当图像信息丢失过多时，可以从公式(1)中推导出，可能是地物抵消破碎云的影响的能力太弱，抑或是云的分布过分离散。
参考图5可知，本发明提供的基于云分布特点的遥感图像地物完整性分析方法具体工作过程如下：
获得原始图像后，首先进行云检测，区分出云区域和地物区域，然后分别对这两者进行分析。首先进行云区域分析，结合统计特性计算云破碎度cf；然后进行地物区域分析，在用分形维数等特征完成地物类型判别的基础上，计算地物的信息损失抵抗度dr；根据计算出的云破碎度cf和地物的信息损失抵抗度dr计算出图像信息损失度ild，并且分析遥感图像的地物完整性。若信息损失度ild较高，则遥感图像地物完整性较差，不适合后续的陆地分析应用；若信息损失度ild较低，则遥感图像地物完整性较好，适合用于多种陆地分析应用的场合，如森林防火、土地资源监测、震后救灾等。因此，本发明提供的分析方法具有能够准确智能描述遥感图像质量的能力。
上文所述的一系列的详细说明仅仅是针对本发明的可行性实施方式的具体说明，它们并不是用以限制本发明的保护范围，在所属技术领域普通技术人员所具备的知识范围内，在不脱离本发明宗旨的前提下作出的各种变化均属于本发明的保护范围。