CN105046087A - 一种遥感卫星多光谱影像的水体信息自动提取方法 - Google Patents

Abstract

一种遥感卫星多光谱影像的水体信息自动提取方法如下：(1)对国产遥感卫星多光谱影像的表观反射率数据进行归一化水体指数提取，形成归一化水体指数影像NDWI；(2)对归一化水体指数影像进行基于最大类间方差算法的全局阈值分割初步提取水体信息；(3)对初步的水体信息进行局部缓冲区扩张和阈值分割，获取稳定的水体单元结果集合；(4)根据高程数据获取坡度和阴影信息并对水体单元结果集合进行阴影移除，最终完成水体信息提取。本发明采用计算机自动解译的方法，利用从全局分割到局部迭代分割的思想，对归一化水体指数流程化地进行水体提取，提取精度高，工程容易。

Description

一种遥感卫星多光谱影像的水体信息自动提取方法

技术领域

本发明涉及一种遥感卫星多光谱影像的水体信息自动提取方法，适用于无云覆盖的国产遥感卫星多光谱影像的水体信息提取，属于遥感信息处理领域

背景技术

随着经济的发展，全球环境破坏也日趋严重，干旱洪涝、水环境污染等问题频繁发生严重影响了人们的生活和社会经济的发展以及社会的稳定，利用遥感信息快速、科学、准确的对水资源进行监测和评价，适时掌握水资源信息，及时采取相应对策，对水资源合理规划及保护有着重要的作用。

水体的遥感监测多是通过地表水域面积、深度的变化进行监测，要想知道水域面积的变化必然离不开水体信息的提取，相对于植被、建筑物、城市道路等类型地物的提取而言，从遥感影像上自动提取水体的方法相对比较少，目前单波段阈值法、多波段谱间关系法、归一化水体指数法等水体判别方法以及根据形状信息来进行水体识别和分类的方法相继得到了应用。

其中单波段阈值法是利用某种地物与背景地物在某一波段上的反射率或像元灰度值的差异确定某一数值为区分该地物和背景地物的方法。此方法原理简单操作简便。但其中关键的是阈值的确定，阈值选取的准确性直接决定了最终水体提取的准确性。因此单波段阈值法在地物类型丰富，地物在所选波段上灰度值接近的影像上具有一定的局限性，很容易提取非水体地物，其提取精度低。

多波段谱间关系法的实质是构造波段运算函数对影像进行处理，该方法能够利用多波段的优势综合提取水体信息。此方法综合利用了多个波段的光谱信息，因此提取效果往往要比单波段阈值法要好。但是此方法要根据不同遥感卫星多光谱影像中水体独特的多波段谱间关系特征，构造出谱间关系法水体提取计算模型，如G+R＞NIR+MIR或MIR/G＜a等，G代表绿光波段，R代表红光波段，NIR代表近红外波段，MIR代表短波红外波段，a为阈值。由于不同遥感卫星多光谱数据往往需要不同的计算模型，因此这种方法的普适性不高，很难推广。

归一化水体指数法利用多光谱影像中绿光波段与近红外波段的特点，构成影像中的水体计算公式：NDWI＝(G－NIR)/(G+NIR)，式中：G表示绿光波段，NIR表示近红外波段。此方法采用了归一化比值运算，能够自动消除地形起伏的影响，在一定程度上区分水体和阴影。但事实上在很多情况下用NDWI提取的水体信息中仍然夹杂着许多非水体信息，特别是在提取山区的水体方面不太理想，水体和阴影存在混淆。而在平原地区该方法提取水体效果较好。

总的来说，以上几种方法存在以下几方面缺点：(1)方法较为原始，提取精度低；(2)在提取过程中需要辅助人工干预，甚至还需要实地实测、手工勾绘，费时费力；(3)针对不同遥感卫星多光谱数据需要不同的计算模型和阈值，普适性低；(4)对于有阴影干扰的山区地区提取效果差，工程应用实现难度较大。

发明内容

本发明解决的技术问题为：克服现有技术不足，提供一种遥感卫星多光谱影像的水体信息自动提取方法，该方法采取全局-局部的思想精确地提取水体信息，并能有效地去除山体阴影等干扰。整个方法采用自动处理的方式，流程简单，工程实现容易。

本发明解决的技术方案为：一种遥感卫星多光谱影像的水体信息自动提取方法，包括水体指数提取阶段、水体指数阈值分割阶段和山体阴影剔除阶段，

所述水体指数提取阶段步骤如下：

(1)选取一景无云的国产卫星多光谱遥感影像，根据该遥感影像的绿光波段的灰度值计算表观反射率值ρ_Green，根据近红外波段的灰度值计算表观反射率ρ_NIR，计算过程如下：

先按照如下公式将影像的灰度值转换为表观辐亮度L_a：

L_a＝Gain×DN+Bias，

式中Gain为增益，DN为灰度值，Bias为偏移量，这些参数从影像和影像头文件中获取；

然后，将表观辐射亮度L_a转换为表观反射率ρ_a，公式如下：

${\mathrm{\ρ}}_a=\frac{\mathrm{\π}\×L_a\×d^2}{E_s\×{\mathrm{cos\θ}}_s},$ 即得到ρ_Green和ρ_NIR；

式中d是日-地距离订正因子，E_s是大气外太阳光谱辐照度，θ_s是太阳天顶角，这些参数可从影像头文件和中国资源卫星应用中心网站获得；

(2)根据步骤(1)中计算出的遥感影像的绿光波段的表观反射率ρ_Green和近红外波段的表观反射率ρ_NIR获取归一化水体指数影像NDWI，公式如下：

$NDWI=\frac{{\mathrm{\ρ}}_{Green}-{\mathrm{\ρ}}_{NIR}}{{\mathrm{\ρ}}_{Green}+{\mathrm{\ρ}}_{NIR}}$

该归一化水体指数影像包括水体与非水体目标，水体与非水体目标一般在水体指数影像直方图上呈双峰分布；由于归一化水体指数影像NDWI是范围为[-1，1]的离散浮点数，对归一化水体指数影像NDWI乘以100并四舍五入，得到水体指数影像，该水体指数影像像元值是范围为[-100，100]的整数，至此完成水体指数提取阶段；

所述水体指数阈值分割阶段步骤如下：

(3)将步骤(2)中水体指数影像中像元的值从小到大依次排列，得到水体指数影像中像元的值的集合；

(4)设定水体指数影像的阈值T₀，该T₀为步骤(3)水体指数影像中像元的值的集合中的第i个值，小于阈值T₀的值为非水体目标，大于阈值T₀的值为水体目标；i＝1…P,P为大于1的正整数；水体目标的像元数占水体指数影像的所有像元数的比例为ω₀，所有水体目标的像元的平均值为μ₀，非水体目标的像元数占水体指数影像的所有像元数的比例为ω₁，所有非水体目标的像元的平均值为μ₁，水体指数影像的平均像元值为μ,水体与非水体目标间的方差记为g，i从1到N，依次计算水体指数影像的平均像元值μ，公式如下

μ＝ω₀×μ₀+ω₁×μ₁

然后计算水体与非水体目标间的方差记为g，公式如下：

g＝ω₀×(μ₀-μ)²+ω₁×(μ₁-μ)²

记录每组阈值T₀对应的水体与非水体目标间的方差g，将最大的水体与非水体目标间的方差g对应的阈值T₀，作为全局阈值T；

(5)将水体指数影像中大于等于步骤(4)的全局阈值T的像元设置为水体目标，将水体指数影像中小于步骤(4)的全局阈值T的像元设置为非水体目标；去除孤立的单个水体目标，将所有相邻接的水体目标设置为一个水体单元；

对每个水体单元的边界各向外每次扩展一个像元，直到扩展的像元数之和与水体单元像元数的差最小，此时扩展的水体单元作为疑似水体单元；

(6)将疑似水体单元中像元的值从小到大依次排列，得到疑似水体单元NDWI中像元的值的集合；

(7)设定疑似水体单元的阈值T_j，该T_j为步骤(6)疑似水体单元中像元的值的集合中的第j个值，小于阈值T_j的值为非水体目标，大于阈值T_j的值为水体目标；j＝1…L,L为大于1的正整数；水体目标的像元数占疑似水体单元的所有像元数的比例为ω_j0，所有水体目标的像元的平均值为μ_j0，非水体目标的像元数占疑似水体单元的所有像元数的比例为ω_j1，所有非水体目标的像元的平均值为μ_j1，疑似水体单元的平均像元值为μ_j,水体与非水体目标间的方差记为g_j，j从1到L，依次计算疑似水体单元的平均像元值μ_j，公式如下

μ_j＝ω_j0×μ_j0+ω_j1×μ_j1

然后计算水体与非水体目标间的方差记为g_j，公式如下：

g_j＝ω_j0×(μ_j0-μ_j)²+ω_j1×(μ_j1-μ_j)²

记录每组阈值T_j与水体与非水体目标间的方差g_j，将最大的水体与非水体目标间的方差g_j对应的阈值T_j，作为局部阈值T_h；

(8)将疑似水体单元中大于等于步骤(7)的局部阈值T_h的像元设置为新的疑似水体单元，再对该新的疑似水体单元的边界向外每次扩展一个像元，直到扩展的像元数之和与新的疑似水体单元像元数的差最小，此时扩展的新的水体单元再次作为疑似水体单元；

(9)将步骤(8)得到的疑似水体单元，代替步骤(6)中的疑似水体单元，重复进行步骤(6)至步骤(8)，直至步骤(8)得到的新的疑似水体单元与扩展之前的疑似水体单元像元数之差小于扩展之前的疑似水体单元像元数的1/100为止，此时，设此新的疑似水体单元为最终疑似水体单元，该最终疑似水体单元可能包含山体阴影；

所述山体阴影剔除阶段步骤如下：

(10)将步骤(9)中新的疑似水体单元在国产卫星多光谱遥感影像相同地理范围的DEM数据中标记出来，形成标记的疑似水体信息，DEM为国产卫星多光谱遥感影像范围内的地面高程数据；

(11)从步骤(1)中的一景无云的国产卫星多光谱遥感影像的元数据中获取太阳方位角L，即太阳光线在地平面上的投影与当地子午线的夹角；

(12)计算DEM数据中每一个像元点对给定太阳方位角L方向的坡度β_L；

(13)从步骤(1)中的一景无云的国产卫星多光谱遥感影像的元数据中获取太阳高度角θ，将DEM数据中每一个像元的太阳方位角L方向的坡度β_L与太阳高度角θ进行比较，如果β_L＞θ，则DEM数据中的此像元以及太阳光线经过此像元在地平面上的投影为阴影区域，将DEM数据中所有阴影区域标记为疑似阴影信息；

(14)在DEM数据中，对比步骤(10)所有标记的疑似水体信息与步骤(13)的标记的疑似阴影信息，如果标记的疑似水体信息被完全包含在标记的疑似阴影信息里，则视此标记的疑似水体信息为阴影信息，则从步骤(9)中最终疑似水体单元中去除对应的阴影信息，最终完成水体信息提取。

所述步骤(4)中在国产遥感卫星影像中归一化水体指数影像NDWI的水体与非水体的初始值阈值一般为0.1，因此此处设置T₀的初始值为10。

所述步骤(12)中计算DEM数据中每一个像元点对给定太阳方位角L方向的坡度β_L，公式如下：

${\mathrm{\β}}_L=\mathrm{\β}\×\frac{f_{x}cos\mathrm{\φ}+f_{y}sin\mathrm{\φ}}{\sqrt{{f_x}^2+{f_y}^2}},$

式中φ为L与北方向顺时针方向夹角，范围取[0°，360°]；β是反映单位长度上高程的升降值，即某点梯度方向的坡度，β表达式为：

$\mathrm{\β}=arctan\sqrt{{f_x}^2+{f_y}^2},$

式中，坡度β范围取[-90°，90°]，式中，f_x为X轴方向的地形曲线函数，f_y为Y轴方向的地形曲线函数，一般使用三阶反距离平方权差分的计算表达式表示，公式为：

$\left\{\begin{array}{c}{f}_x=(z_7-z_1+2(z_8-z_2)+z_9-z_3)/8r\\ f_y=(z_3-z_1+2(z_6-z_4)+z_9-z_7)/8r\end{array}\right.$

式中z_i为DEM一个3×3像元窗口中某像元点的高程，i为此3×3像元窗口各像元点的编号，i＝1,2,…9，r为DEM像元的分辨率。

本发明与现有技术相比的优点在于：

(1)本发明使用遥感影像的表观反射率代替影像的灰度值来提取水体指数，因为遥感成像是在大气中对地面物体辐射进行探测和收集，大气是遥感信息传输的必经介质，太阳辐射在大气传输过程中会与大气发生一系列的相互作用，从而影响卫星传感器入瞳处所记录的地表辐射亮度，也就是最后获取的遥感影像会在一定程度上偏离其本来的地表面貌，因此对遥感影像进行大气订正，即使用表观反射率代替原遥感影像中的灰度值可以在一定程度上还原地表真实面貌，能更加精确地提取遥感影像中的水体指数；

(2)本发明采用最大类间方差法对水体指数影像设置分割阈值。在归一化水体指数影像NDWI中，水体信息得到最大亮度的增强，而其他背景地物受到普遍的抑制，这样使得水体目标和非水体目标在归一化水体指数影像直方图中呈双峰分布，这就为最大类间方差法确定水体目标和非水体目标之间的阈值提供了基础。最大类间方差法是基于图像中目标和背景两类别的可分离性提出的，它是建立在假设图像中目标和背景构成的混合密度函数是由两个服从等方差正态分布的子分布构成基础之上的，对于归一化水体指数影像NDWI，它将分割影像分为前景和背景两类，并求出每一类的类内方差及类间方差，选取使类间方差最大、类内方差最小的最佳阈值，较为清晰地区分水体与非水体之间的界限；

(3)本发明采用全局-局部的分割思想精确分割水体。由于初次的类间方差法阈值分割是针对整幅影像而言，而即使在同一幅影像中不同的水体单元其反射波谱特性也存在一些细微差别，这些差别主要体现在水体边缘或含有杂质的水体，所以使用单一阈值分割水体指数影像常常只能突出水体的主要轮廓，而对于湖水的边缘，细长的河流和含有泥沙的江河等往往只显示其中的主体部分，对于水体的细节不能完全提取。因此本发明在全局分割后再针对每一个水体单元进行局部的阈值分割，这样每一个水体单元的细节就能完整的提取出来，大大提高了水体提取的精度；

(4)本发明使用缓冲区扩张法对每一个水体单元进行局部阈值分割。由于对水体指数影像的初次阈值分割往往会漏掉许多水体的细节部分，所以本发明在初次分割后对每一个疑似水体单元反复进行扩张与分割，直至该疑似水体单元稳定为止，这样可以使水体单元的边界更加精确；

(5)本发明利用坡度标记阴影信息排除水体单元中的阴影。由于山体阴影往往是干扰水体信息提取的最大因素，因此本发明根据水体提取的目标影像对应的高程数据标记阴影区域，由于这些阴影区域往往是陡峭的斜坡，因此排除大范围湖水等静态水体存在的可能，而对于流动的水体如瀑布等，其斜坡顶部非阴影区域一定存在河流等水体。因此，在本发明中完全被标记的阴影区域包含在内的水体单元将被视为阴影而排除，这样就大大降低了水体提取过程中错检的概率。

附图说明

图1为本发明的水体信息自动提取方法流程图；

图2为本发明水体光谱曲线图；

图3为本发明归一化水体指数影像的双峰分布直方图；

图4为本发明坡度三维立体示意图；

图5为本发明中太阳投射阴影示意图；

图6为本发明中DEM数据3×3像元窗口示意图。

具体实施方式

下面结合附图对本发明做进一步详细描述，如图1所示，本发明一种遥感卫星多光谱影像的水体信息自动提取方法的流程图，在水体指数提取阶段，首先对一景无云的国产卫星多光谱遥感影像进行大气订正，因为遥感成像是在大气中对地面物体辐射进行探测和收集，太阳辐射在大气传输过程中会与大气发生一系列的相互作用，包括大气吸收、散射等，从而导致太阳辐射传输特性发生改变，影响卫星传感器入瞳处所记录的地表辐射亮度，这对后续的水体指数提取结果会造成一定的偏差，所以在获取精确的水体指数值之前还原影像中地表辐射特性的真实面貌将显得至关重要。在大气订正中，对归一化水体指数、比值植被指数、土壤调整植被指数等暗目标的指数值计算通常先使用在大气顶部等效地物反射率，即表观反射率对原始影像的灰度值进行转换，其转换过程是先按照如下公式将影像的灰度值转换为表观辐亮度L_a：

L_a＝Gain×DN+Bias，

式中Gain为增益，DN为灰度值，Bias为偏移量，这些参数可从影像和影像头文件中获取。

接着将表观辐射亮度转换为表观反射率，公式如下：

${\mathrm{\ρ}}_a=\frac{\mathrm{\π}\×L_a\×d^2}{E_s\×{\mathrm{cos\θ}}_s},$

式中d是日-地距离订正因子，E_s是大气外太阳光谱辐照度，θ_s是太阳天顶角，这些参数可从影像头文件和中国资源卫星应用中心网站获得。本发明根据该遥感影像的绿光波段的灰度值计算表观反射率值ρ_Green，近红外波段的灰度值计算表观反射率ρ_NIR。

然后根据地物反射波谱的特性可知，利用地物反射率与可见光波长之间关系可以有效地识别区分各种不同地物。由于地物反射波谱不仅能反映地物的反射率随波长的变化规律，而且较为直观地反映出某一特定波长不同地物的反射率差异。如图2所示，水体的反射率主要在波长为0.6μm-0.7μm间的蓝绿光波段，其他波段吸收都很强，特别到波长大于0.9μm的近红外波段吸收就更强了，因此在影像上特别是近红外影像上水体呈黑色。水体指数则是根据这一特性，利用波段比值法原理，使水体信息明显的得到最大亮度增强，而其他背景地物受到普遍的抑制。由于水体对入射光具有强吸收性，在大部分遥感传感器的波长范围内，呈现较弱的反射率，且随波长的增加逐渐减小，在近红外和中红外范围吸收性最强，几乎无反射，因此用可见光波段和近红外波段或中红外波段的反差构成的水体指数可以使水体信息得到加强，非水体信息得到抑制，从而达到突出水体信息目的。

因此根据以上步骤中计算出的遥感影像的绿光波段的表观反射率ρ_Green和近红外波段的表观反射率ρ_NIR获取归一化水体指数影像NDWI，公式如下：

$NDWI=\frac{{\mathrm{\ρ}}_{Green}-{\mathrm{\ρ}}_{NIR}}{{\mathrm{\ρ}}_{Green}+{\mathrm{\ρ}}_{NIR}}$

该归一化水体指数影像包括水体与非水体目标，水体与非水体目标一般在水体指数影像直方图上呈双峰分布，如图3所示；由于归一化水体指数影像NDWI是范围为[-1，1]的离散浮点数，为了后续计算方便，对归一化水体指数影像NDWI乘以100并四舍五入，得到水体指数影像，其像元值是范围为[-100，100]的整数。至此完成水体指数影像的获取阶段。

要想从水体指数影像中提取水体单元，需要对水体指数影像进行阈值分割，本发明采用最大类间方差法进行分割，该方法是基于图像中目标和背景两类别的可分离性提出的,它是建立在假设图像中目标和背景构成的混合密度函数是由两个服从等方差正态分布的子分布构成基础之上的。对图像I(x,y)，前景(即目标)和背景的分割阈值记作T，属于前景点数占整幅图像的比例为ω₀，其平均灰度为μ₀；背景像素点数占整幅图像比例为ω₁，平均灰度为μ₁。图像的总平均灰度为μ：类间方差记为g。

μ＝ω₀×μ₀+ω₁×μ₁

g＝ω₀×(μ₀-μ)²+ω₁×(μ₀-μ)²

假设图像的背景较暗，并且图像的大小为M×N，图像中像素的灰度值小于阈值T的像素个数记作N₀，像素灰度大于阈值T的像素个数记作N₁，则有：

${\mathrm{\ω}}_0=\frac{N_0}{M\×N}$

${\mathrm{\ω}}_1=\frac{N_1}{M\×N}$

ω₀+ω₁＝1

N₀+N₁＝M×N

由上式可知类间方差公式为:

g＝ω₀×ω₁×(ω₀+ω₁)×(μ₀+μ₁)²

通过在水体指数影像的像元值范围内采用遍历的方法即可得到使类间方差g最大的阈值T。

具体方法为首先将水体指数影像中像元的值按小到大依次排列，得到水体指数影像中像元的值的集合；设定水体指数影像的阈值T₀，在国产遥感卫星影像中归一化水体指数影像NDWI的水体与非水体的临界点一般为0.1，为了减少计算复杂度，本发明将归一化水体指数影像NDWI的水体与非水体阈值设置为0.1，因此此处设置T₀的初始值为10，该T₀为水体指数影像中像元的值的集合中的第i个值，小于阈值T₀的值为非水体目标，大于阈值T₀的值为水体目标；i＝1…P,P为大于1的正整数；水体目标的像元数占水体指数影像的所有像元数的比例为ω₀，所有水体目标的像元的平均值为μ₀，非水体目标的像元数占水体指数影像的所有像元数的比例为ω₁，所有非水体目标的像元的平均值为μ₁，水体指数影像的平均像元值为μ,水体与非水体目标间的方差记为g，i从1到N，依次计算水体指数影像的平均像元值μ，公式如下：

μ＝ω₀×μ₀+ω₁×μ₁

然后计算水体与非水体目标间的方差记为g，公式如下：

g＝ω₀×(μ₀-μ)²+ω₁×(μ₁-μ)²

记录每组阈值T₀对应的水体与非水体目标间的方差g，将最大的水体与非水体目标间的方差g对应的阈值T₀，作为全局阈值T；

将水体指数影像中大于等于上一步骤的全局阈值T的像元设置为水体目标，将水体指数影像中小于上一步骤的全局阈值T的像元设置为非水体目标；去除孤立的单个水体目标，这里指只有1个像元构成的水体单元，然后利用连通性划分为一个个水体单元；

由于初次的类间方差法阈值分割是针对整幅影像而言，而即使在同一幅影像中不同的水体单元其反射波谱特性也存在一些细微差别，这些差别主要体现在水体边缘或含有杂质的水体，所以使用单一阈值分割水体指数影像常常只能突出水体的主要轮廓，而对于湖水的边缘，细长的河流和含有泥沙的江河等往往只显示其中的主体部分，对于水体的细节不能完全提取。因此本发明在全局分割后再针对每一个水体单元进行局部的阈值分割，其基本思想是对水体单元做局部缓冲区扩张，再对缓冲区结果做直方图分割，得到该水体单元新的像元数N和新的分割阈值T，当-0.01N₀<(N-N₀)<0.01N₀时，可认为水体单元分割稳定，输出水体单元结果，否则继续重复迭代过程直到稳定再输出结果，其中N₀为缓冲区扩张前的水体单元像元数。

具体方法为首先对每个水体单元的边界各向外每次扩展一个像元，直到扩展的像元数之和与水体单元像元数的差最小，此时扩展的水体单元作为疑似水体单元；将疑似水体单元中像元的值从小到大依次排列，得到疑似水体单元NDWI中像元的值的集合；设定疑似水体单元的阈值T_j，该T_j为上一步骤中疑似水体单元中像元的值的集合中的第j个值，小于阈值T_j的值为非水体目标，大于阈值T_j的值为水体目标；j＝1…L,L为大于1的正整数；水体目标的像元数占疑似水体单元的所有像元数的比例为ω_j0，所有水体目标的像元的平均值为μ_j0，非水体目标的像元数占疑似水体单元的所有像元数的比例为ω_j1，所有非水体目标的像元的平均值为μ_j1，疑似水体单元的平均像元值为μ_j,水体与非水体目标间的方差记为g_j，j从1到L，依次计算疑似水体单元的平均像元值μ_j，公式如下

μ_j＝ω_j0×μ_j0+ω_j1×μ_j1

然后计算水体与非水体目标间的方差记为g，公式如下：

g_j＝ω_j0×(μ_j0-μ_j)²+ω_j1×(μ_j1-μ_j)²

记录每组阈值T_j与水体与非水体目标间的方差g_j，将最大的水体与非水体目标间的方差g_j对应的阈值T_j，作为局部阈值T_h；

将疑似水体单元中大于等于上一步骤的局部阈值T_h的像元设置为新的疑似水体单元，再对该新的疑似水体单元的边界向外每次扩展一个像元，直到扩展的像元数之和与新的疑似水体单元像元数的差最小，此时扩展的新的水体单元再次作为疑似水体单元；

然后对该疑似水体单元重复进行上述扩张和分割的操作，直至新的疑似水体单元与扩展之前的疑似水体单元像元数之差小于扩展之前的疑似水体单元像元数的1/100为止。此时，设此新的疑似水体单元为最终疑似水体单元。至此完成了水体指数影像阈值分割阶段，但此时该最终疑似水体单元可能包含山体阴影。

为了去除可能存在的山体阴影，本发明使用国产卫星多光谱遥感影像范围内的地面高程数据，即DEM来获取山体的阴影信息。

如图4所示，山体阴影主要是太阳光经过山体而投射到地面的阴影，跟据太阳高度角、太阳方位角和坡度可以判断该区域是否存在阴影。太阳高度角是指太阳光的入射方向和地平面之间的夹角，太阳方位角是指太阳光线在地平面上的投影与当地子午线的夹角，这两个值可直接从遥感影像的元数据中获得。坡度是指地表单元陡缓的程度，如图5所示用几何表示是指地表P点法线方向与此点铅垂方向之间的夹角。一般使用DEM数据获取某地面点梯度方向的坡度公式如下：

$\mathrm{\&beta;}=arctan\sqrt{{f_x}^2+{f_y}^2},$

坡度β范围取[-90°，90°]，式中f_x为X轴方向的地形曲线函数，f_y为Y轴方向的地形曲线函数，一般使用三阶反距离平方权差分的计算表达式表示，公式为：

$\left\{\begin{array}{c}{f}_x=(z_7-z_1+2(z_8-z_2)+z_9-z_3)/8r\\ f_y=(z_3-z_1+2(z_6-z_4)+z_9-z_7)/8r\end{array}\right.$

式中z_i为DEM一个3×3像元窗口中某像元点的高程，如图6所示，X轴方向为正北方向，i为此3×3像元窗口各像元点的编号，i＝1,2,…9，r为DEM像元的分辨率。

上式中坡度β为梯度方向上的坡度，对于任意给定方向L，设L与北方向顺时针方向夹角，即方位角为φ，由数学分析可知，曲面z＝f(x,y)沿L方向的方向导数为曲面函数z＝f(x,y)沿L方向的坡度β_L公式为 ${\mathrm{\&beta;}}_L=\mathrm{\&beta;}\&times;\frac{f_{x}cos\mathrm{\&phi;}+f_{y}sin\mathrm{\&phi;}}{\sqrt{{f_x}^2+{f_y}^2}}.$

将太阳方位角代入公式可求出DEM每个地面点沿太阳照射方向的坡度，再与太阳高度角相比较可得知该区域是否为阴影区域，当太阳高度角小于坡度时，此点为阴影区域，如图4所示，通过此点高程h与太阳高度角可求出阴影长度。至此可在DEM数据中得知所有的阴影区域。

利用以上DEM数据中标记阴影区域可以排除水体单元中的阴影。因为这些阴影区域往往是陡峭的斜坡，因此排除大范围湖水等静态水体存在的可能，而对于流动的水体如瀑布等，其斜坡顶部非阴影区域一定存在河流等水体。因此，在本发明中完全被标记的阴影区域包含在内的水体单元将被视为阴影而排除，这样就大大降低了水体提取过程中错检的概率。至此完成本发明目标水体信息的提取。

本发明未详细阐述部分属于本领域公知技术。

Claims (3)

1.一种遥感卫星多光谱影像的水体信息自动提取方法，其特征在于：包括水体指数提取阶段、水体指数阈值分割阶段和山体阴影剔除阶段，

所述水体指数提取阶段步骤如下：

(1)选取一景无云的国产卫星多光谱遥感影像，根据该遥感影像的绿光波段的灰度值计算表观反射率值ρ_Green，根据近红外波段的灰度值计算表观反射率ρ_NIR，计算过程如下：

先按照如下公式将影像的灰度值转换为表观辐亮度L_a：

L_a＝Gain×DN+Bias，

式中Gain为增益，DN为灰度值，Bias为偏移量；

然后，将表观辐射亮度L_a转换为表观反射率ρ_a，公式如下：

${\mathrm{\&rho;}}_a=\frac{\mathrm{\&pi;}\&times;L_a\&times;d^2}{E_s\&times;{\mathrm{cos\&theta;}}_s},$ 即得到ρ_Green和ρ_NIR；

式中d是日-地距离订正因子，E_s是大气外太阳光谱辐照度，θ_s是太阳天顶角，这些参数可从影像头文件和中国资源卫星应用中心网站获得；

(2)根据步骤(1)中计算出的遥感影像的绿光波段的表观反射率ρ_Green和近红外波段的表观反射率ρ_NIR获取归一化水体指数影像NDWI，公式如下：

$NDWI=\frac{{\mathrm{\&rho;}}_{Green}-{\mathrm{\&rho;}}_{NIR}}{{\mathrm{\&rho;}}_{Green}+{\mathrm{\&rho;}}_{NIR}}$

该归一化水体指数影像包括水体与非水体目标，水体与非水体目标一般在水体指数影像直方图上呈双峰分布；由于归一化水体指数影像NDWI是范围为[-1，1]的离散浮点数，对归一化水体指数影像NDWI乘以100并四舍五入，得到水体指数影像，该水体指数影像像元值是范围为[-100，100]的整数，至此完成水体指数提取阶段；

所述水体指数阈值分割阶段步骤如下：

(3)将步骤(2)中水体指数影像中像元的值从小到大依次排列，得到水体指数影像中像元的值的集合；

(4)设定水体指数影像的阈值T₀，该T₀为步骤(3)水体指数影像中像元的值的集合中的第i个值，小于阈值T₀的值为非水体目标，大于阈值T₀的值为水体目标；i＝1…P,P为大于1的正整数；水体目标的像元数占水体指数影像的所有像元数的比例为ω₀，所有水体目标的像元的平均值为μ₀，非水体目标的像元数占水体指数影像的所有像元数的比例为ω₁，所有非水体目标的像元的平均值为μ₁，水体指数影像的平均像元值为μ,水体与非水体目标间的方差记为g，i从1到N，依次计算水体指数影像的平均像元值μ，公式如下

μ＝ω₀×μ₀+ω₁×μ₁

然后计算水体与非水体目标间的方差记为g，公式如下：

g＝ω₀×(μ₀-μ)²+ω₁×(μ₁-μ)²

记录每组阈值T₀对应的水体与非水体目标间的方差g，将最大的水体与非水体目标间的方差g对应的阈值T₀，作为全局阈值T；

(5)将水体指数影像中大于等于步骤(4)的全局阈值T的像元设置为水体目标，将水体指数影像中小于步骤(4)的全局阈值T的像元设置为非水体目标；去除孤立的单个水体目标，将所有相邻接的水体目标设置为一个水体单元；

对每个水体单元的边界各向外每次扩展一个像元，直到扩展的像元数之和与水体单元像元数的差最小，此时扩展的水体单元作为疑似水体单元；

(6)将疑似水体单元中像元的值从小到大依次排列，得到疑似水体单元NDWI中像元的值的集合；

(7)设定疑似水体单元的阈值T_j，该T_j为步骤(6)疑似水体单元中像元的值的集合中的第j个值，小于阈值T_j的值为非水体目标，大于阈值T_j的值为水体目标；j＝1…L,L为大于1的正整数；水体目标的像元数占疑似水体单元的所有像元数的比例为ω_j0，所有水体目标的像元的平均值为μ_j0，非水体目标的像元数占疑似水体单元的所有像元数的比例为所有非水体目标的像元的平均值为疑似水体单元的平均像元值为μ_j,水体与非水体目标间的方差记为g_j，j从1到L，依次计算疑似水体单元的平均像元值μ_j，公式如下

μ_j＝ω_j0×μ_j0+ω_j1×μ_j1

然后计算水体与非水体目标间的方差记为g_j，公式如下：

g_j＝ω_j0×(μ_j0-μ_j)²+ω_j1×(μ_j1-μ_j)²

记录每组阈值T_j与水体与非水体目标间的方差g_j，将最大的水体与非水体目标间的方差g_j对应的阈值T_j，作为局部阈值T_h；

(8)将疑似水体单元中大于等于步骤(7)的局部阈值T_h的像元设置为新的疑似水体单元，再对该新的疑似水体单元的边界向外每次扩展一个像元，直到扩展的像元数之和与新的疑似水体单元像元数的差最小，此时扩展的新的水体单元再次作为疑似水体单元；

(9)将步骤(8)得到的疑似水体单元，代替步骤(6)中的疑似水体单元，重复进行步骤(6)至步骤(8)，直至步骤(8)得到的新的疑似水体单元与扩展之前的疑似水体单元像元数之差小于扩展之前的疑似水体单元像元数的1/100为止，此时，设此新的疑似水体单元为最终疑似水体单元，该最终疑似水体单元可能包含山体阴影；

所述山体阴影剔除阶段步骤如下：

(10)将步骤(9)中新的疑似水体单元在国产卫星多光谱遥感影像相同地理范围的DEM数据中标记出来，形成标记的疑似水体信息，DEM为国产卫星多光谱遥感影像范围内的地面高程数据；

(11)从步骤(1)中的一景无云的国产卫星多光谱遥感影像的元数据中获取太阳方位角L，即太阳光线在地平面上的投影与当地子午线的夹角；

(12)计算DEM数据中每一个像元点对给定太阳方位角L方向的坡度β_L；

(13)从步骤(1)中的一景无云的国产卫星多光谱遥感影像的元数据中获取太阳高度角θ，将DEM数据中每一个像元的太阳方位角L方向的坡度β_L与太阳高度角θ进行比较，如果β_L>θ，则DEM数据中的此像元以及太阳光线经过此像元在地平面上的投影为阴影区域，将DEM数据中所有阴影区域标记为疑似阴影信息；

(14)在DEM数据中，对比步骤(10)所有标记的疑似水体信息与步骤(13)的标记的疑似阴影信息，如果标记的疑似水体信息被完全包含在标记的疑似阴影信息里，则视此标记的疑似水体信息为阴影信息，则从步骤(9)中最终疑似水体单元中去除对应的阴影信息，最终完成水体信息提取。

2.根据权利要求1所述的一种遥感卫星多光谱影像的水体信息自动提取方法，其特征在于：所述步骤(4)中在国产遥感卫星影像中归一化水体指数影像NDWI的水体与非水体的初始值阈值一般为0.1，因此此处设置T₀的初始值为10。

3.根据权利要求1所述的一种遥感卫星多光谱影像的水体信息自动提取方法，其特征在于：所述步骤(12)中计算DEM数据中每一个像元点对给定太阳方位角L方向的坡度β_L，公式如下：

${\mathrm{\&beta;}}_L=\mathrm{\&beta;}\&times;\frac{f_{x}cos\mathrm{\&phi;}+f_{y}sin\mathrm{\&phi;}}{\sqrt{{f_x}^2+{f_y}^2}},$

式中φ为L与北方向顺时针方向夹角，范围取[0°，360°]；β是反映单位长度上高程的升降值，即某点梯度方向的坡度，β表达式为：

$\mathrm{\&beta;}=arctan\sqrt{{f_x}^2+{f_y}^2},$

式中，坡度β范围取[-90°，90°]，f_x为X轴方向的地形曲线函数，f_y为Y轴方向的地形曲线函数，使用三阶反距离平方权差分的计算表达式表示，公式为：

$\left\{\begin{array}{c}{f}_x=(z_7-z_1+2(z_8-z_2)+z_9-z_3)/8r\\ f_y=(z_3-z_1+2(z_6-z_4)+z_9-z_7)/8r\end{array}\right.$

式中，z_i为DEM一个3×3像元窗口中某像元点的高程，i为此3×3像元窗口各像元点的编号，i＝1,2,…9，r为DEM像元的分辨率。