一种遥感卫星多光谱影像的水体信息自动提取方法
【技术领域】
[0001] 本发明涉及一种遥感卫星多光谱影像的水体信息自动提取方法，适用于无云覆盖 的国产遥感卫星多光谱影像的水体信息提取，属于遥感信息处理领域
【背景技术】
[0002] 随着经济的发展，全球环境破坏也日趋严重，干旱洪涝、水环境污染等问题频繁发 生严重影响了人们的生活和社会经济的发展以及社会的稳定，利用遥感信息快速、科学、准 确的对水资源进行监测和评价，适时掌握水资源信息，及时采取相应对策，对水资源合理规 划及保护有着重要的作用。
[0003] 水体的遥感监测多是通过地表水域面积、深度的变化进行监测，要想知道水域面 积的变化必然离不开水体信息的提取，相对于植被、建筑物、城市道路等类型地物的提取 而言，从遥感影像上自动提取水体的方法相对比较少，目前单波段阈值法、多波段谱间关系 法、归一化水体指数法等水体判别方法以及根据形状信息来进行水体识别和分类的方法相 继得到了应用。
[0004] 其中单波段阈值法是利用某种地物与背景地物在某一波段上的反射率或像元灰 度值的差异确定某一数值为区分该地物和背景地物的方法。此方法原理简单操作简便。但 其中关键的是阈值的确定，阈值选取的准确性直接决定了最终水体提取的准确性。因此单 波段阈值法在地物类型丰富，地物在所选波段上灰度值接近的影像上具有一定的局限性， 很容易提取非水体地物，其提取精度低。
[0005] 多波段谱间关系法的实质是构造波段运算函数对影像进行处理，该方法能够利用 多波段的优势综合提取水体信息。此方法综合利用了多个波段的光谱信息，因此提取效果 往往要比单波段阈值法要好。但是此方法要根据不同遥感卫星多光谱影像中水体独特的多 波段谱间关系特征，构造出谱间关系法水体提取计算模型，如G+R > NIR+MIR或MIR/G < a 等，G代表绿光波段，R代表红光波段，NIR代表近红外波段，MIR代表短波红外波段，a为阈 值。由于不同遥感卫星多光谱数据往往需要不同的计算模型，因此这种方法的普适性不高， 很难推广。
[0006] 归一化水体指数法利用多光谱影像中绿光波段与近红外波段的特点，构成影像中 的水体计算公式：NDWI = (G - NIRV(G+NIR)，式中：G表示绿光波段，NIR表示近红外波 段。此方法采用了归一化比值运算，能够自动消除地形起伏的影响，在一定程度上区分水体 和阴影。但事实上在很多情况下用NDWI提取的水体信息中仍然夹杂着许多非水体信息，特 别是在提取山区的水体方面不太理想，水体和阴影存在混淆。而在平原地区该方法提取水 体效果较好。
[0007] 总的来说，以上几种方法存在以下几方面缺点：（1)方法较为原始，提取精度低； (2)在提取过程中需要辅助人工干预，甚至还需要实地实测、手工勾绘，费时费力；(3)针对 不同遥感卫星多光谱数据需要不同的计算模型和阈值，普适性低；(4)对于有阴影干扰的 山区地区提取效果差，工程应用实现难度较大。
【发明内容】
[0008] 本发明解决的技术问题为：克服现有技术不足，提供一种遥感卫星多光谱影像的 水体?目息自动提取方法，该方法米取全局 -局郃的思想精确地提取水体?目息，并能有效地 去除山体阴影等干扰。整个方法采用自动处理的方式，流程简单，工程实现容易。
[0009] 本发明解决的技术方案为：一种遥感卫星多光谱影像的水体信息自动提取方法， 包括水体指数提取阶段、水体指数阈值分割阶段和山体阴影剔除阶段，
[0010] 所述水体指数提取阶段步骤如下：
[0011] (1)选取一景无云的国产卫星多光谱遥感影像，根据该遥感影像的绿光波段的灰 度值计算表观反射率值P ，根据近红外波段的灰度值计算表观反射率Pnir，计算过程如 下：
[0012] 先按照如下公式将影像的灰度值转换为表观辐亮度La:
[0013] La= GainXDN+Bias，
[0014] 式中Gain为增益，DN为灰度值，Bias为偏移量，这些参数从影像和影像头文件中 获取；
[0015] 然后，将表观辐射亮度La转换为表观反射率P a，公式如下：
[0017] 式中d是日-地距离订正因子，Es是大气外太阳光谱辐照度，Θ 5是太阳天顶角， 这些参数可从影像头文件和中国资源卫星应用中心网站获得；
[0018] (2)根据步骤⑴中计算出的遥感影像的绿光波段的表观反射率P 和近红外 波段的表观反射率Pnir获取归一化水体指数影像NDWI，公式如下：
[0020] 该归一化水体指数影像包括水体与非水体目标，水体与非水体目标一般在水体指 数影像直方图上呈双峰分布；由于归一化水体指数影像NDWI是范围为[_1，1]的离散浮点 数，对归一化水体指数影像NDWI乘以100并四舍五入，得到水体指数影像，该水体指数影像 像元值是范围为[-100,100]的整数，至此完成水体指数提取阶段；
[0021] 所述水体指数阈值分割阶段步骤如下：
[0022] (3)将步骤⑵中水体指数影像中像元的值从小到大依次排列，得到水体指数影 像中像元的值的集合；
[0023] (4)设定水体指数影像的阈值T。，该T。为步骤（3)水体指数影像中像元的值的集 合中的第i个值，小于阈值T。的值为非水体目标，大于阈值T。的值为水体目标；i = 1···Ρ，P 为大于1的正整数；水体目标的像元数占水体指数影像的所有像元数的比例为ω。，所有水 体目标的像元的平均值为μ c，非水体目标的像元数占水体指数影像的所有像元数的比例 为CO1，所有非水体目标的像元的平均值为μ i，水体指数影像的平均像元值为μ，水体与非 水体目标间的方差记为g，i从1到Ν，依次计算水体指数影像的平均像元值μ，公式如下
[0024] μ = CO0Xyt^co1Xy1
[0025] 然后计算水体与非水体目标间的方差记为g，公式如下：
[0026] g =
[0027] 记录每组阈值T。对应的水体与非水体目标间的方差g，将最大的水体与非水体目 标间的方差g对应的阈值T。，作为全局阈值T ;
[0028] (5)将水体指数影像中大于等于步骤⑷的全局阈值T的像元设置为水体目标，将 水体指数影像中小于步骤（4)的全局阈值T的像元设置为非水体目标；去除孤立的单个水 体目标，将所有相邻接的水体目标设置为一个水体单元；
[0029] 对每个水体单元的边界各向外每次扩展一个像元，直到扩展的像元数之和与水体 单元像元数的差最小，此时扩展的水体单元作为疑似水体单元；
[0030] (6)将疑似水体单元中像元的值从小到大依次排列，得到疑似水体单元NDWI中像 元的值的集合；
[0031] (7)设定疑似水体单元的阈值Tj，该Tj为步骤（6)疑似水体单元中像元的值的集 合中的第j个值，小于阈值T j的值为非水体目标，大于阈值T j的值为水体目标；j = 1-L，L 为大于1的正整数；水体目标的像元数占疑似水体单元的所有像元数的比例为CO j。，所有水 体目标的像元的平均值为μ j(]，非水体目标的像元数占疑似水体单元的所有像元数的比例 为ωη，所有非水体目标的像元的平均值为μ U，疑似水体单元的平均像元值为P j,水体与 非水体目标间的方差记为gj，j从1到L，依次计算疑似水体单元的平均像元值μ j，公式如 下
[0032] μ J= ω j0X μ j0+ω.ηΧ μ .η
[0033] 然后计算水体与非水体目标间的方差记为gj，公式如下：
[0034] gj = ω j〇 X ( μ j〇- μ j)2+ ω η X ( μ η- μ j)2
[0035] 记录每组阈值Tj与水体与非水体目标间的方差g j，将最大的水体与非水体目标间 的方差gi对应的阈值T j，作为局部阈值Th;
[0036] (8)将疑似水体单元中大于等于步骤（7)的局部阈值Th的像元设置为新的疑似水 体单元，再对该新的疑似水体单元的边界向外每次扩展一个像元，直到扩展的像元数之和 与新的疑似水体单元像元数的差最小，此时扩展的新的水体单元再次作为疑似水体单元；
[0037] (9)将步骤