CN109300133B - 一种城市河网区水体提取方法 - Google Patents

Abstract

本发明公开了一种城市河网区水体提取方法。城市区域存在大量的建筑物阴影，有的阴影位于水上，有的位于陆地上。由于阴影光谱和水体相似，在进行水体提取时，经常会误提许多阴影区。本发明选择位于阴影区的水体作为暗像元，由于被阴影遮挡，阴影区水体没有太阳直射光照射，可以认为其出水辐射为0，利用其计算成像时的气溶胶参数，由此对图像上的所有水体像元进行大气校正；利用大气校正后的归一化水体指数配合单波段阈值法，可以非常有效地去除城市区域存在的大量建筑物阴影，快速准确地提取城市区域的水体。本发明解决了城市河网区提取水体中难以去除大量城市阴影的问题。

Description

一种城市河网区水体提取方法

技术领域

本发明属于遥感信息提取技术领域，涉及一种基于高分辨率遥感图像，进行城市河网区水体提取的方法。

背景技术

城市河网区承担着蓄积雨洪、调节行洪和滨水景观美化等方面的功能，是社会经济发展和生态环境建设的基础与支撑。准确地提取地表水体位置、面积及相应的变化信息，是开展城市水资源调查、水利规划、环境监测与保护等工作的先决条件。

由于卫星遥感数据监测具有范围广、实时动态性强的特点，遥感技术能够大大降低获取地表水体信息的时间成本与人力成本，在水域监测的相关领域得到了越来越多的关注。但是，受限于遥感数据源的庞杂以及水体提取算法的局限性，基于卫星遥感影像的水体提取精度尚不能比肩传统实地调查手段得到的结果。这是遥感技术在该领域业务化应用过程中亟需改进的关键所在。

目前，遥感水体提取主要基于水体的光谱特征和空间关系进行分析，应用较为广泛的有阈值法、多波段运算法、谱间关系法、区域生长法等。其中，阈值法根据水体在某个特征波段(一般为诸如近红外波段的长波波段)的反射率与其余地物之间的差异，通过图像阈值分割算法实现水体范围的提取；多波段运算法基于绿光、近红外波段反射率构建归一化差异水体指数，抑制植被信息而突出水体信息；谱间关系法综合考虑了水体与其他地物的波谱曲线差异，以决策树分类等方法进行水体识别；区域生长法则是通过人工选取纯水体像元作为种子，通过相邻像元与种子像元的关系来确定关系。

将上述方法应用于大范围的开阔水域提取时一般均可取得较好的结果，但是，应用场景变为城市河网区时，则普遍出现水体提取精度下降的情况。这是因为城市水域形态参差不齐且河涌宽度较小，对遥感影像空间分辨率的要求较高，进而凸显出以下三方面的问题：1)高分辨率遥感影像可以清晰地观察到城市高层建筑的阴影，而阴影与水体光谱特征容易混淆，使得大量的阴影像元被误判为水体像元；2)城市河网区的水体光学类型较多，诸如清洁水体、富营养化水体、浑浊水体、黑臭水体等的光谱形态各不相同，容易出现水体像元的漏提；3)通过影像辐射定标、大气校正等预处理过程“还原”地物实际光谱形态，是所有遥感水体识别方法的根本，但是目前宽谱段、高分辨率卫星传感器的设计初衷多是针对陆地观测，势必降低了水体像元预处理结果的精度。

为此，研究一种针对城市河网区，具有较高精度的水体识别方法具有重要的价值。

发明内容

针对以上问题，本发明基于高分辨率遥感影像，综合了水体像元的大气校正算法和光谱特征，提出一种城市河网区水体提取方法，有效地解决了阴影像元的干扰，提高了城市河网区水体识别的精度。

本发明的目的通过以下的技术方案实现：一种城市河网区水体提取方法，包括针对遥感图像，选择城市河网区中因临水建筑形成的阴影区域作为暗像元，认为其出水辐射为0；利用暗像元推算各个波段气溶胶参数，进而进行所有波段的大气校正；根据上述校正后图像数据获取水体离水反射率，利用水体离水反射率计算归一化水体指数，通过归一化水体指数判定是否是水体像元。本发明可以非常有效地减少城市区域大量建筑物阴影的干扰，快速准确地提取城市区域的水体。

优选的，在判断是否是水体像元时，判断归一化水体指数是否大于第一阈值，且单波段水体离水反射率是否低于第二阈值，同时满足时判定当前像元为水体像元。这样归一化水体指数配合单波段阈值法进行综合判断，识别更为准确。

具体的，包括以下步骤：

A、对遥感图像数据进行预处理，获得辐射亮度图像；

B、在辐射亮度图像中选择若干个城市河网区中被涉水建筑物遮挡的水体阴影像元作为暗像元，统计所有选中的暗像元在各个波段辐亮度的最小值，认为最小值等于天空光在水表形成的漫反射和大气散射之和；

C、利用步骤B中所得各个波段辐亮度的最小值，计算成像时的气溶胶浑浊度系数及波长指数，由此推算各个波段的气溶胶光学厚度；

D、用步骤A所述辐射亮度图像各个波段减去步骤B中所得各个波段辐亮度的最小值，再除以各个波段大气上行透过率，得到水体的离水辐亮度，完成大气校正；

E、将离水辐亮度换算成离水反射率，利用近红外波段及绿光波段离水反射率计算归一化水体指数NDWI；

F、像元的归一化水体指数NDWI大于第一阈值且近红外波段水体离水反射率低于第二阈值，判别其为水体像元。

优选地，所述步骤A具体包括：

A1：获取卫星数据的辐射定标参数：国产卫星数据每年数据定标参数都有变化，一般在中国资源卫星中心官网(URL：http://www.cresda.com/CN/)公布；国外数据定标参数一般从头文件中读取。

A2：对遥感图像数据进行预处理：根据卫星数据公布的辐射定标方法及定标参数，对原始数据进行辐射定标，获得辐射亮度图像L(λ)。对于水体像元而言：

L(λ)＝L_w(λ)T_↑(λ)+rL_sky(λ)T_↑(λ)+L_g(λ)+L_p(λ) (1)

式中，λ为卫星传感器各通道的中心波长，单位μm；L_w(λ)是进入水体的光被水体散射回来的离水辐亮度，包含了水质信息；T_↑为沿着观测方向的大气上行透过率；rL_sky(λ)T_↑(λ)为天空光被水面反射后被卫星观测到的辐亮度，没有任何水体信息，其中L_sky(λ)为天空光辐亮度，r为气水界面反射率，一般取值0.025；L_g(λ)为水面对太阳直射光的镜面反射，传感器成像时的观测天顶角和太阳天顶角一般不等，认为可避开水面对太阳光的镜面反射，该项为0；L_p(λ)为大气分子、气溶胶等颗粒散射形成的大气程辐射。因此，在忽略水体镜面反射时，卫星传感器接收水体总的辐亮度L(λ)为：

L(λ)＝L_w(λ)T_↑(λ)+rL_sky(λ)T_↑(λ)+L_p(λ) (2)

优选地，所述步骤B具体包括：

B1：尽可能多的选取辐射亮度图像上位于阴影区的水体阴影像元，包括建筑物、桥梁等在水中形成的阴影。

B2：统计选中的暗像元在各个波段的最小辐亮度值L_shd(λ)(一般而言，在大气情况相同的情况下，阴影像元的反射率应当近似相等。但是由于周围环境光、混合像元、传感器辐射误差等原因，选中的阴影像元反射率存在些许的差异。此处取最小值是为了避免出现过校正的情况)，认为该值等于天空光在水表形成的漫反射和大气散射之和，如图3所示，即：

L_shd(λ)＝rL_sky(λ)T_↑(λ)+L_p(λ) (3)

优选地，所述步骤C具体包括：

C1：从遥感数据的头文件中读取成像时的太阳天顶角θ、卫星观测天顶角

太阳方位角φ_s、卫星方位角φ_v及大气层位太阳辐照度E₀(λ)。

C2：利用水体阴影的最小辐亮度值L_shd(λ)计算气溶胶光学厚度τ_a(λ)，具体求解过程如下：

首先，大气散射形成的程辐射辐亮度值可表示为：

其中，P(Θ)为散射相函数，Θ为散射角。成像时天气晴朗，气溶胶以细粒为主，气体分子以及气溶胶的散射类型均为瑞利散射，则：

为整层大气散射系数，它等于大气分子散射系数ω_r(λ)和气溶胶散射系数ω_a(λ)之和，即

而大气整层光学厚度

α(λ)为整层大气吸收系数。由于多波段卫星传感器的工作波段均位于大气窗口，大气吸收作用很小，忽略大气吸收有：

则式(4)可改写为：

其次，天空光漫射到地表面的光谱辐照度

c为散射光下行分量的比例，成像时天气晴朗，气溶胶以细粒为主，气体分子以及气溶胶的散射类型均为瑞利散射，则c＝1/2。假设天空光辐射各项同向性，则：

τ(λ)＝τ_a(λ)+τ_r(λ) (10)

τ_r(λ)＝0.008524λ^-4(1+0.0113λ^-2+0.00013λ^-4) (11)

气溶胶光学厚度根据

公式可表示为：

τ_a(λ)＝βλ^-α (12)

式中，α为气溶胶波长指数，与气溶胶平均半径有关，平均半径越小，其值越小；β为混浊度系数，表征了整层大气中气溶胶的数量。

然后，大气上行透过率T_↑(λ)可根据比尔定律表示为：

最后，将式8、9、10、12、13代入式3可知，L_shd(λ)是气溶胶波长指数α、混浊度系数β的函数，即：

常用高分辨率卫星遥感影像一般包括蓝光、绿光、红光、近红外四个波段，对应上述2个未知参数，可构成一个非线性最小二乘问题。同时，为降低非线性系统求解过程中带来的异常解，这里参考相关文献(Sasithorn，1983；Cornette and Shanks，1992；Gilabertet al，1993；王跃思等，2006；陈启东等，2011；于兴娜等，2012；徐舟等，2015)对未知参数α、β设置上下限，分别为[0.2,4]、[0.008,1.5]。其中，气溶胶波长指数α的上限4对应气溶胶以细粒子为主的情况，下限0.2对应较为灰霾的天气状况；气溶胶浑浊度系数β的上限、下限则分别对应气溶胶浓度较高和天气晴朗的状况。故此，以α、β为未知数x，构造目标函数f(x)(式15)，形成带边界的非线性优化问题(式16)进行求解。

f(x)＝F(α，β)-L_shd(λ) (15)

优选地，所述步骤D具体包括：

D1：成像时天气晴朗，认为待校正水域入射天空光辐亮度L_sky(λ)一致，用原始辐射亮度图像各个波段L(λ)减去步骤B中统计的L_shd(λ)，根据式(2)、式(3)可得：

L(λ)-L_shd(λ)＝L_w(λ)T_↑(λ) (17)

D2：根据步骤C中所得的气溶胶光学厚度，带入式(13)求得大气上行透过率T_↑(λ)，进而得到离水辐亮度L_w(λ)：

优选地，所述步骤E具体包括：

E1：获取水体离水反射率，则根据下式计算：

E2：利用水体离水反射率计算归一化水体指数：

其中，R_w(NIR)、R_w(GREEN)表示近红外波段、绿光波段的离水反射率。

本发明与现有技术相比，具有如下优点和有益效果：

本发明的城市河网区水体提取方法，不需复杂的图像处理及阈值分析，仅仅基于图像上的水体阴影进行大气校正，利用大气校正后的图像计算水体指数配合单波段阈值法就可以实现城市水体提取。而且水体阴影像元的选取简单易行，没有针对于不同图像的阈值判定。该方法应用于高楼密布的城市河网区提取时，能减少临水高层建筑带来的阴影干扰，具有快速有效的优点。

附图说明

图1是本发明方法的流程图。

图2是传感器接收的水体信号组成示意图。

图3是传感器接收到的位于阴影区的水体信号组成示意图。

图4是实例1中原始影像假彩色合成图。

图5是实例1中采用现有技术辐射定标后NDWI指数提取水体的效果图。

图6是实例1中采用现有技术FLAASH大气校正后NDWI指数提取水体的效果图。

图7是实例1中采用本发明方法提取水体的效果图。

图8是实例2中原始影像假彩色合成图。

图9是实例2中采用现有技术辐射定标后NDWI指数提取水体的效果图。

图10是实例2中采用现有技术FLAASH大气校正后NDWI指数提取水体的效果图。

图11是实例2中采用本发明方法提取水体的效果图。

图12是实例3中原始影像假彩色合成图。

图13是实例3中采用现有技术辐射定标后NDWI指数提取水体的效果图。

图14是实例3中采用现有技术FLAASH大气校正后NDWI指数提取水体的效果图。

图15是实例3中采用本发明方法提取水体的效果图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

城市区域存在大量的建筑物阴影，有的阴影位于水上，有的阴影位于陆地上，由于其光谱和水体相似，在进行水体提取时，经常会误提许多阴影区。本发明针对遥感图像，选择城市河网区中因临水建筑形成的阴影区域作为暗像元，认为其出水辐射为0；利用暗像元推算各个波段气溶胶参数，进而进行所有波段的大气校正；根据上述校正后图像数据获取水体离水反射率，利用水体离水反射率计算归一化水体指数，通过归一化水体指数，配合单波段阈值法，判定是否是水体像元。本发明可以非常有效地去除城市区域存在的大量建筑物阴影，快速准确的提取城市区域的水体。

本实施例以广州、澳门作为城市河网区的示例，以高分2号数据、SPOT7数据、worldview2数据作为示例，详细说明利用本发明方法进行水体提取的过程。根据图1所示的方法流程图，方法包括步骤：

1：购得研究区2016年12月18日的高分2号遥感影像，并从中国资源卫星应用中心官网上查询对应辐射定标参数。

2：对原始数据进行辐射定标，获得辐射亮度图像，记辐亮度L(λ)。

如图2所示，对于水体像元而言：

L(λ)＝L_w(λ)T_↑(λ)+rL_sky(λ)T_↑(λ)+L_g(λ)+L_p(λ) (1)

式中，λ为卫星传感器各通道的中心波长，单位μm；L_w(λ)是进入水体的光被水体散射回来后的离水辐亮度，包含了水质信息；T_↑为沿着观测方向的大气上行透过率；rL_sky(λ)T_↑(λ)为天空光被水面反射后被卫星观测到的辐亮度，没有任何水体信息，其中L_sky(λ)为天空光辐亮度，r为气水界面反射率，一般取值0.025；L_g(λ)为水面对太阳直射光的镜面反射，传感器成像时的观测天顶角和太阳天顶角一般不等，认为可避开水面对太阳光的镜面反射，该项为0；L_p(λ)为大气分子、气溶胶等颗粒散射形成的大气程辐射。因此，在忽略水体镜面反射时，卫星传感器接收水体总的辐亮度L(λ)为：

L(λ)＝L_w(λ)T_↑(λ)+rL_sky(λ)T_↑(λ)+L_p(λ) (2)。

3：尽可能多的选取图像上由高层建筑、桥梁等造成的水体阴影区像元。该步骤可以通过ArcGIS软件，用面状矢量把图像中的水体阴影区勾画出来。

4：一般而言，在大气情况相同的情况下，阴影像元的反射率应当近似相等。但是由于周围环境光、混合像元、传感器辐射误差等原因，选中的阴影像元反射率存在些许的差异。为了避免出现过校正的情况，统计所有选中的阴影像元辐亮度在各个波段的最小值L_shd(λ)，认为该值等于天空光在水表形成的漫反射和大气散射之和，如图3所示，即：

L_shd(λ)＝rL_sky(λ)T_↑(λ)+L_p(λ) (3)

该步骤可以通过ArcGIS软件的统计功能，利用勾画出来的阴影区矢量图层对图像辐射亮度L(λ)进行最小值统计。

5：从遥感数据头文件中读取成像时的太阳天顶角θ、卫星观测天顶角

太阳方位角φ_s、卫星方位角φ_v及大气层外太阳辐照度E₀(λ)。

6：根据上述参数，按以下步骤构建有边界的非线性优化方程组求解气溶胶光学厚度τ_a(λ)。

首先，大气散射形成的程辐射辐亮度值可表示为：

其中，P(Θ)为散射相函数，Θ为散射角。成像时天气晴朗，气溶胶以细粒为主，气体分子以及气溶胶的散射类型可近似为瑞利散射，则：

为整层大气散射系数，它等于大气分子散射系数ω_r(λ)和气溶胶散射系数ω_a(λ)之和，即

而大气整层光学厚度

α(λ)为整层大气吸收系数。由于多波段卫星传感器的工作波段均位于大气窗口，大气吸收作用很小，忽略大气吸收有：

则式(4)可改写为：

其次，天空光漫射到地表面的光谱辐照度

c为散射光下行分量的比例，成像时天气晴朗，气溶胶以细粒为主，气体分子以及气溶胶的散射类型均为瑞利散射，则c＝1/2。假设天空光辐射各项同向性，则：

τ(λ)＝τ_a(λ)+τ_r(λ) (10)

τ_r(λ)＝0.008524λ^-4(1+0.0113λ^-2+0.00013λ^-4) (11)

气溶胶光学厚度根据

公式可表示为：

τ_a(λ)＝βλ^-α (12)

式中，α为气溶胶波长指数，与气溶胶平均半径有关，平均半径越小，其值越小；β为混浊度系数，表征了整层大气中气溶胶的数量。

然后，大气上行透过率T_↑(λ)可根据比尔定律表示为：

最后，综合上述过程可知L_shd(λ)是气溶胶波长指数α、混浊度系数β的函数，即：

常用高分辨率卫星遥感影像一般包括蓝光、绿光、红光、近红外四个波段，对应上述2个未知参数，可构成一个非线性最小二乘问题。同时，为降低非线性系统求解过程中带来的异常解，这里参考现有技术对未知参数α、β设置上下限，分别为[0.2,4]、[0.008,1.5]。其中，气溶胶波长指数α的上限4对应气溶胶以细粒子为主的情况，下限0.2对应较为灰霾的天气状况；气溶胶浑浊度系数β的上限、下限则分别对应气溶胶浓度较高和天气晴朗的状况。故此，以α、β为未知数x，构造目标函数f(x)，形成带边界的非线性优化问题进行求解。

7：成像时天气晴朗，认为待校正水域入射天空光辐亮度L_sky(λ)一致，用原始辐亮度影像各个波段L(λ)减去上述统计的L_shd(λ)，可得：

L(λ)-L_shd(λ)＝L_w(λ)T_↑(λ) (17)

8：根据所得的气溶胶光学厚度，带入式(10)求得大气上行透过率T_↑(λ)，进而得到离水辐亮度：

10：计算水体离水反射率：

11：利用水体离水反射率计算归一化水体指数：

12：在PCI、ENVI、ArcGIS中使用波段运算提取水体，运算条件为：NDWI>0且NIR<0.05。

本发明具体给出3个实例，来显示本发明的城市河网区水体提取方法相较于现有方法的优越性。

参见图4-7，该实例1针对的对象是2016年12月18日广州主城区的高分2号影像，该影像原始影像假彩色合成图如图4所示。第一种比对方法是辐射定标后直接利用水体指数提取，结果如图5所示；第二种比对方法利用ENVIFLAASH模块大气校正后，再用水体指数提取，结果如图6所示；本发明方法提取水体结果如图7所示。结果表明，相对于其他两种方法，本发明方法可以有效地去除大量城市阴影区域。

参见图8-11，该实例2针对的对象是2017年12月29日广州主城区的SPOT7影像，该影像原始影像假彩色合成图如图8所示。第一种比对方法是辐射定标后直接利用水体指数提取，结果如图9所示；第二种比对方法利用ENVI FLAASH模块大气校正后，再用水体指数提取，结果如图10所示；本发明方法提取水体结果如图11所示。结果表明，相对于其他两种方法，本发明方法可以有效地去除大量城市阴影区域。

参见图12-15，该实例3针对的对象是2015年10月21日澳门半岛的worldview2遥感影像，该影像原始影像假彩色合成图如图12所示。第一种比对方法是辐射定标后直接利用水体指数提取，结果如图13所示；第二种比对方法利用ENVI FLAASH模块大气校正后，再用水体指数提取，结果如图14所示；本发明方法提取水体结果如图15所示。结果表明，相对于其他两种方法，本发明方法可以有效地去除大量城市阴影区域。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (7)

1.一种城市河网区水体提取方法，其特征在于，包括以下步骤：

A、对遥感图像数据进行预处理，获得辐射亮度图像；

B、在辐射亮度图像中选择若干个城市河网区中被涉水建筑物遮挡的水体阴影像元作为暗像元，统计所有选中的暗像元在各个波段辐亮度的最小值，认为最小值等于天空光在水表形成的漫反射和大气散射之和；

C、利用步骤B中所得各个波段辐亮度的最小值，计算成像时的气溶胶浑浊度系数及波长指数，由此推算各个波段的气溶胶光学厚度；

D、用步骤A所述辐射亮度图像各个波段减去步骤B中所得各个波段辐亮度的最小值，再除以各个波段大气上行透过率，得到水体的离水辐亮度，完成大气校正；

E、将离水辐亮度换算成离水反射率，利用近红外波段及绿光波段离水反射率计算归一化水体指数NDWI；

F、像元的归一化水体指数NDWI大于第一阈值且近红外波段水体离水反射率低于第二阈值，判别其为水体像元。

2.根据权利要求1所述的城市河网区水体提取方法，其特征在于，所述步骤A具体包括：

A1：获取卫星数据的辐射定标参数；

A2：对遥感图像数据进行预处理：根据卫星数据公布的辐射定标方法及定标参数，对原始数据进行辐射定标，获得辐射亮度图像L(λ)；

L(λ)＝L_w(λ)T_↑(λ)+rL_sky(λ)T_↑(λ)+L_p(λ)

式中，λ为卫星传感器各通道的中心波长；L_w(λ)是进入水体的光被水体散射回来的离水辐亮度，包含了水质信息；T_↑为沿着观测方向的大气上行透过率；rL_sky(λ)T_↑(λ)为天空光被水面反射后被卫星观测到的辐亮度，没有任何水体信息，其中L_sky(λ)为天空光辐亮度，r为气水界面反射率；L_p(λ)为大气分子、气溶胶颗粒散射形成的大气程辐射。

3.根据权利要求2所述的城市河网区水体提取方法，其特征在于，所述步骤B具体包括：

B1：选取多个辐射亮度图像上位于阴影区的水体阴影像元；

B2：统计选中的暗像元在各个波段的最小辐亮度值L_shd(λ)，认为该值等于天空光在水表形成的漫反射和大气散射之和，即：

L_shd(λ)＝rL_sky(λ)T_↑(λ)+L_p(λ)。

4.根据权利要求3所述的城市河网区水体提取方法，其特征在于，所述步骤C具体包括：

C1：从遥感数据的头文件中读取成像时的太阳天顶角θ、卫星观测天顶角

太阳方位角φ_s、卫星方位角φ_v及大气层外太阳辐照度E₀(λ)；

C2：利用水体阴影的最小辐亮度值L_shd(λ)计算气溶胶光学厚度τ_a(λ)，具体求解过程如下：

首先，大气散射形成的程辐射辐亮度值表示为：

其中，τ(λ)表示大气整层光学厚度，令其等于整层大气散射系数

等于大气分子散射系数ω_r(λ)和气溶胶散射系数ω_a(λ)之和，即

P(Θ)为散射相函数，Θ为散射角；

其次，天空光漫射到地表面的光谱辐照度

c为散射光下行分量的比例，假设天空光辐射各项同向性，则：

τ(λ)＝τ_a(λ)+τ_r(λ)

τ_r(λ)＝0.008524λ^-4(1+0.0113λ^-2+0.00013λ^-4)

气溶胶光学厚度根据

公式表示为：

τ_a(λ)＝βλ^-α

式中，α为气溶胶波长指数；β为混浊度系数；

然后，大气上行透过率T_↑(λ)根据比尔定律表示为：

最后，将上述各个参数的表示形式代入L_shd(λ)＝rL_sky(λ)T_↑(λ)+L_p(λ)中，得到：

以α、β为未知数x，构造目标函数f(x)，形成带边界的非线性优化问题进行求解：

f(x)＝F(α，β)-L_shd(λ)

通过上述目标函数得到最优的气溶胶光学厚度τ_a(λ)，lowerBound、upperBound为设置的未知数的上限和下限值。

5.根据权利要求4所述的城市河网区水体提取方法，其特征在于，在进行优化问题求解时，对未知数α、β设置上下限，分别为[O.2，4]、[O.008，1.5]。

6.根据权利要求4或5所述的城市河网区水体提取方法，其特征在于，所述步骤D具体包括：

D1：用原始辐射亮度图像各个波段L(λ)减去统计的L_shd(λ)，得到：

L(λ)-L_shd(λ)＝L_w(λ)T_↑(λ)

D2：根据步骤C中所得的气溶胶光学厚度，带入式(13)求得大气上行透过率T_↑(λ)，进而得到离水辐亮度L_w(λ)：

7.根据权利要求6所述的城市河网区水体提取方法，其特征在于，所述步骤E具体包括：

E1：获取水体离水反射率，则根据下式计算：

E2：利用水体离水反射率计算归一化水体指数：

其中，R_w(NIR)、R_w(GREEN)表示近红外波段、绿光波段的离水反射率。

Images