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

摘要

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

主权项

一种遥感卫星多光谱影像的水体信息自动提取方法，其特征在于：包括水体指数提取阶段、水体指数阈值分割阶段和山体阴影剔除阶段，所述水体指数提取阶段步骤如下：(1)选取一景无云的国产卫星多光谱遥感影像，根据该遥感影像的绿光波段的灰度值计算表观反射率值ρ_Green，根据近红外波段的灰度值计算表观反射率ρ_NIR，计算过程如下：先按照如下公式将影像的灰度值转换为表观辐亮度L_a：L_a＝Gain×DN+Bias，式中Gain为增益，DN为灰度值，Bias为偏移量；然后，将表观辐射亮度L_a转换为表观反射率ρ_a，公式如下：${\rho }_a=\frac{\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{{\rho }_{Green}-{\rho }_{NIR}}{{\rho }_{Green}+{\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)中最终疑似水体单元中去除对应的阴影信息，最终完成水体信息提取。