CN104881659A - 一种不透水层的提取方法及装置 - Google Patents

Abstract

本发明提供了一种不透水层的提取方法及装置，包括：获取卫星影像数据；计算影像数据在λ波段的大气层辐亮度ρ_λ，将影像数据预处理为辐亮度影像数据；根据辐亮度影像数据中的绿光波段数据和近红外波段数据计算调整归一化水体指数MNDWI；当MNDWI不大于0时，在辐亮度影像数据中提取纯净像元数据；利用建筑线方程、土壤线方程及土壤阴影线方程确定PBI的系数n；利用公式PBI＝nBlue-NIR计算PBI指数；根据均值标准差法确定PBI指数的不透水层的提取阈值T；其中，当PBI指数大于阈值T时，确定影像为不透水层区域；所述λ取值为自然数；如此，不透水层的提取方法构建的PBI指数除了能够区分不透水层与水体、植被地物，还能够有效区分各类不透水层与裸土及阴影，降低不透水层的提取误差。

Description

一种不透水层的提取方法及装置

技术领域

本发明属于遥感影像处理技术领域，尤其涉及一种不透水层的提取方法及装置。

背景技术

近年来，城市作为人口最为集聚和人类活动最为强烈的地区，其环境和气候等日益受到关注；城市的扩展直接改变了土地覆盖、土地利用的状况，也直接影响局部的气候、生态和环境。不透水层作为城市化进程和城市化率评估的重要指标，对分析评价城市地区的气候、环境、水循环等具有重要的意义。

现有技术中，利用遥感影像进行的不透水层提取方法主要包括：监督分类法、混合像元分解法、谱间关系法、不透水指数法等；但这些方法在进行大数据的处理过程中主要依靠人工提取，不仅耗时耗力，而且人为主观影响严重。实际应用中，普遍采用的方法是指数提取法，但指数提取法利用已有的不透水层指数在提取不透水层时无法将裸土部分剔除，导致结果不精确，在后续的提取过程中产生较大误差；特别是在较干旱、裸土地表较多的西北地区，以及东部地区旱季影像上的提取中产生更大的误差。

基于此，目前亟需一种新型的提取方法，重新构建能够有效区分裸土的不透水层提取指数，降低不透水层的提取误差。

发明内容

针对现有技术存在的问题，本发明实施例提供了一种不透水层的提取方法及装置。

本发明提供一种不透水层的提取方法，所述方法包括：

获取卫星影像数据；

计算所述影像数据在λ波段的大气层辐亮度ρ_λ，将所述影像数据预处理为辐亮度影像数据；

根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI；

当所述MNDWI不大于0时，在所述辐亮度影像数据中选取各类地物的纯净像元数据；

利用建筑线方程、土壤线方程及土壤阴影线方程确定垂直不透水层指数PBI的系数n；

利用公式PBI＝nBlue-NIR计算PBI指数；

根据均值标准差法确定所述PBI指数的不透水层的提取阈值T；其中，

当所述PBI指数大于所述阈值T时，确定影像像元为不透水层区域；所述λ取值为自然数。

上述方案中，所述计算所述影像数据在λ波段的大气层辐亮度ρ_λ包括：

根据公式ρ′_λ＝M_ρQ_cal-A_ρ计算未校正的λ波段大气层辐亮度ρ′_λ；

根据公式对所述未校正的λ波段大气层辐亮度ρ′_λ进行校正，获取校正后的λ波段的大气层辐亮度ρ_λ；其中，所述M_ρ为增益系数；所述Q_cal为λ波段的所述影像像元亮度DN值；所述A_ρ为偏置系数；所述θ_SE为当地的太阳高度角。

上述方案中，所述根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI包括：

利用公式计算所述调整归一化水体指数MNDWI；其中，所述Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中近红外波段像元数据。

上述方案中，所述利用建筑线方程确定垂直不透水层指数PBI的系数n包括：

对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y＝a_bx+b；

当确定所述不透水层所在区域存在山体阴影时，构建土壤和山体阴影样本数据，对所述土壤和山体阴影样本数据进行最小二乘拟合，确定土壤阴影线方程y＝a_ssx+b；

利用公式计算所述垂直不透水层指数PBI的系数n；其中，n与a相等。

上述方案中，当确定所述不透水层所在区域不存在山体阴影时，所述方法还包括：

构建土壤样本数据，对所述土壤样本数据进行最小二乘拟合，确定土壤线方程y＝a_sx+b；

利用公式计算所述垂直不透水层指数PBI的系数n；其中，n与a′相等。

上述方案中，所述根据均值标准差法确定所述PBI指数的不透水层提取阈值T包括：

对辐亮度影像数据进行统计获取PBI影像像元的均值μ和标准差ε；

利用公式T＝μ+ε确定所述PBI指数的不透水层提取阈值T。

本发明同时还提供了一种不透水层的提取装置，所述装置包括：

获取单元，所述获取单元用于获取卫星影像数据；

预处理单元，所述预处理单元用于计算所述影像数据在λ波段的大气层光谱反射率ρ_λ，将所述影像数据预处理为辐亮度影像数据；

计算单元，所述计算单元用于根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI；

选取单元，所述选取单元用于确定当所述调整归一化水体指数不大于0时，在所述辐亮度影像数据中选取各类地物的纯净像元数据；

第一确定单元，所述第一确定单元用于利用建筑线方程、土壤线方程及土壤阴影线方程确定垂直不透水层指数PBI的系数n；

第二确定单元，所述第二确定单元用于利用公式PBI＝nBlue-NIR计算PBI指数；

第三确定单元，所述第三确定单元用于根据均值标准差法确定所述PBI指数的不透水层的提取阈值T；

第四确定单元，所述第四确定单元用于确定当所述PBI指数大于所述阈值T时，确定影像像元为不透水层区域；其中，所述λ取值为自然数。

上述方案中，所述预处理单元具体用于：

根据公式ρ′_λ＝M_ρQ_cal-A_ρ计算未校正的λ波段大气层辐亮度ρ′_λ；

根据公式对所述未校正的λ波段大气层辐亮度ρ′_λ进行校正，获取校正后的λ波段的大气层辐亮度ρ_λ；其中，所述M_ρ为增益系数；所述Q_cal为λ波段的DN值；所述A_ρ为偏置系数；所述θ_SE为当地的太阳高度角。

上述方案中，所述计算单元具体用于：

利用公式计算所述调整归一化水体指数MNDWI；其中，所述Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中的近红外波段像元数据。

上述方案中，所述第一确定单元具体用于：

对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y＝a_bx+b；

当确定所述不透水层所在区域存在山体阴影时，构建土壤和山体阴影样本数据，对所述土壤和山体阴影样本数据进行最小二乘拟合，确定土壤阴影线方程y＝a_ssx+b；

利用公式计算所述垂直不透水层指数PBI的系数n；其中，n与a相等；

当确定所述不透水层所在区域不存在山体阴影时，构建土壤样本数据，对所述土壤样本数据进行最小二乘拟合，确定土壤线方程y＝a_sx+b；

利用公式计算所述垂直不透水层指数PBI的系数n；其中，n与a′相等。

本发明提供了一种不透水层的提取方法及装置，所述方法包括：获取卫星影像数据；计算所述影像数据在λ波段的大气层辐亮度ρ_λ，将所述影像数据预处理为辐亮度影像数据；根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI，当所述MNDWI不大于0时，在所述辐亮度影像数据中选取各类地物的纯净像元数据；利用建筑线方程、土壤线方程及土壤阴影线方程确定垂直不透水层指数PBI的系数n；利用公式PBI＝nBlue-NIR计算PBI指数；根据均值标准差法确定所述PBI指数的不透水层的提取阈值T；其中，当所述PBI指数大于所阈值T时，确定影像像元为不透水层区域；所述λ取值为自然数；如此，所述不透水层的提取方法构建的PBI指数能够有效区分各类不透水层、水体、植被及裸土地表，降低了不透水层的提取误差；并且该方法对不同数据源的卫星影像数据均具有适用性。

附图说明

图1为本发明实施例一提供的不透水层的提取方法流程示意图；

图2为本发明实施例二提供的不透水层的提取装置结构示意图；

图3为本发明实施例三提供的有山体阴影区域的原始影像图；

图4为本发明实施例三提供的有山体阴影的纯净像元样本散点示意图；

图5为本发明实施例三提供的有山体阴影的PBI指数的图像直方图；

图6为本发明实施例三提供的有山体阴影区域的不透水层提取结果示意图；

图7为本发明实施例四提供的无山体阴影区域的原始影像图；

图8为本发明实施例四提供的无山体阴影的纯净像元样本散点示意图；

图9为本发明实施例四提供的无山体阴影的PBI指数的图像直方图；

图10为本发明实施例四提供的无山体阴影的不透水层提取结果示意图。

具体实施方式

为了能够在不透水层的提取中有效区分各类不透水层、水体、植被及裸土地表，降低了不透水层的提取误差，本发明提供了一种不透水层的提取方法及装置，所述方法包括：获取卫星影像数据；计算所述影像数据在λ波段的大气层辐亮度ρ_λ，将所述影像数据预处理为辐亮度影像数据；根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI；当所述MNDWI不大于0时，在所述辐亮度影像数据中选取各类地物纯净像元数据；利用建筑线方程、土壤线方程及土壤阴影线方程确定垂直不透水层指数PBI的系数n；利用公式PBI＝nBlue-NIR计算PBI指数；根据均值标准差法确定所述PBI指数的不透水层的提取阈值T；其中，当所述PBI指数大于所述阈值T时，确定影像像元为不透水层区域；其中，所述λ取值为自然数。

下面通过附图及具体实施例对本发明的技术方案做进一步的详细说明。

实施例一

本实施例提供一种不透水层的提取方法，如图1所示，所述方法主要包括以下步骤：

步骤110，获取卫星影像数据。

本步骤中，可以直接从云服务器或者本地服务器上获取卫星影像数据，所述卫星影像数据可以包括：Landsat卫星影像数据、Aster卫星影像数据等。

步骤111，计算所述影像数据在λ波段的大气层辐亮度ρ_λ，将所述影像数据预处理为辐亮度影像数据。

本步骤中，根据公式(1)计算未校正λ波段的大气层辐亮度ρ′_λ；

ρ′_λ＝M_ρQ_cal-A_ρ  (1)

其中，在公式(1)中，所述M_ρ为增益系数；所述Q_cal为λ波段的影像像元亮度DN值；所述A_ρ为偏置系数，所述λ取值为自然数。

根据公式(2)对所述未校正的λ波段大气层光辐亮度ρ′_λ进行校正，获取校正后的λ波段大气层辐亮度ρ_λ；

${\mathrm{\ρ}}_{\mathrm{\λ}}=\frac{{\mathrm{\ρ}}_{\mathrm{\λ}}^{\′}}{\sin \left({\mathrm{\θ}}_{\mathrm{SE}}\right)}---\left(2\right)$

其中，在公式(2)中，所述θ_SE为当地的太阳高度角。

这里，当λ波段大气层辐亮度ρ_λ计算出以后，所述影像数据即被预处理为辐亮度影像数据。

步骤112，根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI。

本步骤中，所述调整归一化水体指数(MNDWI，Modified NormalizedDifference Water Index)是在归一化水体指数(NDWI，Normalized DifferenceWater Index)的基础上进行的计算，用于移除影像上的水体，减小后续的提取误差。

具体地，根据公式(3)计算所述调整归一化水体指数MNDWI；

$\mathrm{MNDWI}=\frac{\mathrm{Green}-\mathrm{NIR}}{\mathrm{Green}+\mathrm{NIR}}---\left(3\right)$

其中，在公式(3)中，所述Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中的近红外波段像元数据。

步骤113，当所述调整归一化水体指数不大于0时，在所述辐亮度影像数据中提取纯净像元数据。

本步骤中，当所述MNDWI不大于0时，表明此影像区域为非水体区域，在影像中保留此区域，在此影像区域中通过目视解译选取各类典型地物，包括三种不透水层、土壤、植被、水体和阴影的纯净像元作为不透水层纯净像元数据样本；其中，所述纯净像元是指该像元中只包含一种地物，没有其他地物混合。所述纯净像元应尽量选取位于地物中心的像元作为纯净像元。

当所述MNDWI大于0时，表明此影像区域为水体区域，对所述区域进行剔除。

步骤114，利用建筑线方程、土壤线方程和土壤阴影线方程确定垂直不透水层指数PBI的系数n。

本步骤中，因城市中不透水层类型可按其光谱特征的不同分为三种：亮色建筑、暗色建筑和蓝色建筑。其中，亮色和蓝色不透水层在蓝光波段的反射率较高，其余几种地物在蓝光波段反射率低，因此蓝光波段能够有效区分亮、蓝色建筑与裸土。而暗色建筑的光谱值在近红外波段无较大变化，裸土则显著升高，因此可利用蓝波段与近红外波段值的差异来区分暗色建筑与裸土。利用以上几种地物的光谱特性，可以确定垂直不透水层指数PBI。

这里，利用建筑线方程确定垂直不透水层指数PBI的系数n具体包括：在蓝-近红外特征空间内作不透水层的样本散点图，对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y＝a_bx+b；其中，所述建筑线斜率为a_b。

其中，在蓝-近红外特征空间内作不透水层的样本散点图具体包括：在提取了三种不透水层的纯净像元数据样本的基础上，对纯净像元数据样本的蓝、近红外波段值进行导出，以蓝波段像元值为x轴横坐标，以近红外波段像元值为y轴纵坐标，将纯净像元数据样本散点进行作图，得到蓝-近红外特征空间内的三种不透水层的样本散点图。

进一步地，判断该影像区域所在的地区是否有大型山体及山体阴影，当确定所述不透水层所在区域存在山体阴影时，在影像上选取纯净的土壤像元数据及山体影像数据构建土壤和山体阴影样本数据，对所述土壤和山体阴影样本数据进行最小二乘拟合，确定土壤阴影线方程y＝a_ssx+b；其中，所述土壤阴影线斜率为a_ss。

利用公式(4)计算所述垂直不透水层指数PBI的系数n；其中，当n与a相等时，不透水层的提取效果最好。

$a=\frac{a_b+a_{\mathrm{ss}}}{2}---\left(4\right)$

进一步地，当确定所述不透水层所在区域不存在山体阴影时，在影像上选取纯净的土壤像元数据据构建土壤样本数据，对所述土壤样本数据进行最小二乘拟合，确定土壤线方程y＝a_sx+b；其中，所述土壤线斜率为a_s。

利用公式(5)计算所述垂直不透水层指数PBI的系数n；其中，当n与a′相等时，不透水层的提取效果最好。

$a^{\′}=\frac{a_b+a_s}{2}---\left(5\right)$

步骤115，利用公式PBI＝nBlue-NIR计算PBI指数。

本步骤中，当确定出垂直不透水层指数PBI的系数n后，可以利用公式(6)计算PBI指数。

PBI＝nBlue-NIR  (6)

其中，所述Blue为辐亮度影像数据中的蓝光波段像元数据，所述NIR为辐亮度影像数据中的近红外波段像元数据。

步骤116，根据均值标准差法确定所述PBI指数的不透水层提取阈值T。

本步骤中，可以根据均值标准差法确定所述PBI指数的不透水层提取阈值T；具体地，对辐亮度影像数据进行快速统计获取所述PBI影像像元的均值μ和标准差ε；

利用公式(7)确定所述PBI指数的不透水层提取阈值T。

T＝μ+ε  (7)

这里，当阈值T为均值μ与一倍标准差ε之和时，不透水层的提取效果最好。

步骤117，根据所述不透水层提取阈值T对不透水层进行提取。

本步骤中，当所述PBI指数大于所阈值T时，确定所述影像像元为不透水层区域，所述不透水层包括：亮色不透水层、暗色不透水层和蓝色不透水层。

当所述PBI指数不大于所阈值T时，确定所述影像像元为土壤、植被等透水层区域。

本实施例提供的不透水层提取方法，重新构建能够有效区分裸土的不透水层提取指数，除了能够区分不透水层与水体、植被地物，还能够有效区分各类不透水层与裸土及阴影，降低不透水层的提取误差，相比已有的提取方法，在较干旱、裸土地表较多的西北地区，以及东部地区旱季影像上的提取效果都有很大的提高。且对于不同的影像数据源均具有适用性，因此，本实施例提供的不透水层的提取方法有较广泛的应用前景，较强的实用性和推广性。

实施例二

相对应实施例一，本实施例还提供了一种不透水层的提取装置，如图2所示，所述装置包括：获取单元21、预处理单元22、计算单元23、选取单元24、第一确定单元25、第二确定单元26，第三确定单元27、第四确定单元28；其中，

所述获取单元21用于获取卫星影像数据；具体地，所述获取单元21可以直接从云服务器或者本地服务器上获取卫星影像数据，所述卫星影像数据可以包括：Landsat卫星影像数据、Aster卫星影像数据等。

当所述获取单元21获取到卫星影像数据后，所述预处理单元22用于计算所述影像数据在λ波段的大气层光谱反射率ρ_λ，将所述影像数据进行预处理为辐亮度影像数据。具体地，所述预处理单元22根据公式(1)计算未校正λ波段的大气层辐亮度ρ′_λ；

ρ′_λ＝M_ρQ_cal-A_ρ  (1)

其中，在公式(1)中，所述M_ρ为增益系数；所述Q_cal为λ波段的影像像元亮度DN值；所述A_ρ为偏置系数，所述λ取值为自然数。。

根据公式(2)对所述未校正的λ波段大气层光辐亮度ρ′_λ进行校正，获取校正后的λ波段大气层辐亮度ρ_λ；

${\mathrm{\ρ}}_{\mathrm{\λ}}=\frac{{\mathrm{\ρ}}_{\mathrm{\λ}}^{\′}}{\sin \left({\mathrm{\θ}}_{\mathrm{SE}}\right)}---\left(2\right)$

其中，在公式(2)中，所述θ_SE为当地的太阳高度角。

这里，当λ波段的大气层辐亮度ρ_λ计算出以后，所述影像数据即被预处理为辐亮度影像数据。

当所述预处理单元22将影像数据处理为辐亮度影像数据后，所述计算单元23用于根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI。

这里，所述调整归一化水体指数MNDWI是在归一化水体指数NDWI的基础上进行的计算，用于移除影像上的水体，减小后续的提取误差。

具体地，所述计算单元23根据公式(3)计算所述调整归一化水体指数MNDWI；

$\mathrm{MNDWI}=\frac{\mathrm{Green}-\mathrm{NIR}}{\mathrm{Green}+\mathrm{NIR}}---\left(3\right)$

其中，在公式(3)中所述Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中的近红外波段像元数据。

当所述计算单元23计算出所述调整归一化水体指数MNDWI后，所述选取单元24用于判断所述调整归一化水体指数MNDWI是否大于0，当确定所述调整归一化水体指数MNDWI大于0时，表明此影像区域为非水体区域，在影像中保留此区域，通过目视解译选取此影像区域中各类典型地物，包括三种不透水层、土壤、植被、水体和阴影的纯净像元作为不透水层纯净像元数据样本；其中，所述纯净像元是指该像元中只包含一种地物，没有其他地物混合。所述纯净像元应尽量选取位于地物中心的像元作为纯净像元。

当所述选取单元24确定MNDWI大于0时，表明此影像区域为水体区域，对所述区域进行剔除。

当所述选取单元24确定述调整归一化水体指数MNDWI大于0时，所述第一确定单元25利用建筑线方程确定垂直不透水层指数PBI的系数n。

具体地，所述第一确定单元25在此影像区域中确定三种不透水层地物的纯净像元作为不透水层纯净像元数据；在蓝-近红外特征空间内作不透水层的样本散点图，对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y＝a_bx+b；其中，所述建筑线斜率为a_b。

当所述第一确定单元25确定出建筑线方程后，判断该影像区域所在的地区是否有大型山体及山体阴影，当确定所述不透水层所在区域存在山体阴影时，在影像上选取纯净的土壤像元数据及山体影像数据构建土壤和山体阴影样本数据，对所述土壤和山体阴影样本数据进行最小二乘拟合，确定土壤阴影线方程y＝a_ssx+b；其中，所述土壤阴影线斜率为a_ss。

当所述第一确定单元25确定出土壤阴影线方程后，利用公式(4)计算所述垂直不透水层指数PBI的系数n；其中，当n与a相等时，不透水层的提取效果最好。

$a=\frac{a_b+a_{\mathrm{ss}}}{2}---\left(4\right)$

进一步地，当所述第一确定单元25确定所述不透水层所在区域不存在山体阴影时，在影像上选取纯净的土壤像元数据据构建土壤样本数据，对所述土壤样本数据进行最小二乘拟合，确定土壤线方程y＝a_sx+b；其中，所述土壤线斜率为a_s。

所述第一确定单元25确定出土壤线方程后，利用公式(5)计算所述垂直不透水层指数PBI的系数n；其中，当n与a′相等时，不透水层的提取效果最好。

$a^{\′}=\frac{a_b+a_s}{2}---\left(5\right)$

当所述第一确定单元25确定出垂直不透水层指数PBI的系数n后，所述第二确定单元26用于根据公式(6)计算PBI指数。

PBI＝nBlue-NIR  (6)

其中，所述Blue为辐亮度影像数据中的蓝光波段像元数据，所述NIR为辐亮度影像数据中的近红外波段像元数据。

这里，因城市中不透水层类型可按其光谱特征的不同分为三种：亮色建筑、暗色建筑和蓝色建筑。其中，亮色和蓝色不透水层在蓝光波段的反射率较高，其余几种地物在蓝光波段反射率低，因此蓝光波段能够有效区分亮、蓝色建筑与裸土。而暗色建筑的光谱值在近红外波段无较大变化，裸土则显著升高，因此，所述第二确定单元26可利用蓝波段与近红外波段值的差异来区分暗色建筑与裸土；利用以上几种地物的光谱特性，可以确定PBI指数。

当所述第二确定单元26确定出PBI指数后，所述第三确定单元27用于：根据均值标准差法确定所述PBI指数的不透水层提取阈值T。具体地，所述第三确定单元27对辐亮度影像数据进行快速统计获取所述PBI影像像元的均值μ和标准差ε；

利用公式(7)确定所述PBI指数的不透水层提取阈值T。

T＝μ+ε  (7)

这里，当阈值T为均值μ与一倍标准差ε之和时，不透水层的提取效果最好。

当所述第三确定单元27确定出所述阈值T后，所述第四确定单元28具体用于：判断所述PBI指数是否大于所阈值T，当所述第四确定单元28所述PBI指数大于所阈值T时，确定所述影像像元为不透水层区域，所述不透水层包括：亮色不透水层、暗色不透水层和蓝色不透水层。

当所述第四确定单元28确定所述PBI指数不大于所阈值T时，确定所述影像像元为土壤、植被等透水层区域。

实际应用中，所述获取单元21、所述预处理单元22、所述计算单元23、所述选取单元24、所述第一确定单元25、所述第二确定单元26，所述第三确定单元27、所述第四确定单元28可由该装置中的中央处理器(CPU，CentralProcessing Unit)、数字信号处理器(DSP，Digtal Signal Processor)、可编程逻辑阵列(FPGA，Field Programmable Gate Array)、微控制单元(MCU，MicroController Unit)实现。

实施例三

实际应用中，本实施例是以北京市作为有山体阴影的实验区进行不透水层的提取，所述有山体阴影区域的原始影像图如图3所示。

当对北京有山体阴影的实验区进行不透水层的提取时，按照实施例一提供的方法，具体实施如下：

直接从云服务器或者本地服务器上获取北京市的卫星影像数据，所述卫星影像数据可以包括：Landsat卫星影像数据。

根据公式(1)计算未校正的λ波段大气层辐亮度ρ′_λ；

ρ′_λ＝M_ρQ_cal-A_ρ  (1)

其中，在公式(1)中，所述M_ρ为增益系数；所述Q_cal为λ波段的影像像元亮度DN值；所述A_ρ为偏置系数；所述λ取值为自然数。

根据公式(2)对所述未校正的λ波段大气层光辐亮度ρ′_λ进行校正，获取校正后的λ波段大气层辐亮度ρ_λ；

${\mathrm{\ρ}}_{\mathrm{\λ}}=\frac{{\mathrm{\ρ}}_{\mathrm{\λ}}^{\′}}{\sin \left({\mathrm{\θ}}_{\mathrm{SE}}\right)}---\left(2\right)$

其中，在公式(2)中，所述θ_SE为当地的太阳高度角。

这里，当λ波段大气层辐亮度ρ_λ计算出以后，所述影像数据即被预处理为辐亮度影像数据。

根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI。

其中，所述调整归一化水体指数MNDWI是在归一化水体指数NDWI的基础上进行的计算，用于移除影像上的水体，减小后续的提取误差。

具体地，根据公式(3)计算所述调整归一化水体指数MNDWI；

$\mathrm{MNDWI}=\frac{\mathrm{Green}-\mathrm{NIR}}{\mathrm{Green}+\mathrm{NIR}}---\left(3\right)$

在公式(3)中，所述Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中的近红外波段像元数据。

进一步地，当所述调整归一化水体指数不大于0时，在所述辐亮度影像数据中提取纯净像元数据。

具体地，当所述MNDWI不大于0时，表明此影像区域为非水体区域，在影像中保留此区域。

当所述MNDWI大于0时，表明此影像区域为水体区域，对所述区域进行剔除。

进一步地，利用建筑线方程、土壤阴影线方程确定垂直不透水层指数PBI的系数n；

具体地，在上一步保留的区域中通过目视解译选取各类典型地物，包括三种不透水层、土壤、植被、水体和阴影的纯净像元作为不透水层纯净像元数据样本；其中，所述纯净像元是指该像元中只包含一种地物，没有其他地物混合。所述纯净像元应尽量选取位于地物中心的像元作为纯净像元。

这里，因城市中不透水层类型可按其光谱特征的不同分为三种：亮色建筑、暗色建筑和蓝色建筑。其中，亮色和蓝色不透水层在蓝光波段的反射率较高，其余几种地物在蓝光波段反射率低，因此蓝光波段能够有效区分亮、蓝色建筑与裸土。而暗色建筑的光谱值在近红外波段无较大变化，裸土则显著升高，因此可利用蓝波段与近红外波段值的差异来区分暗色建筑与裸土。利用以上几种地物的光谱特性，可以确定垂直不透水层指数PBI。

这里，利用建筑线方程确定垂直不透水层指数PBI的系数n具体包括：

在蓝-近红外特征空间内作不透水层的样本散点图，对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y＝a_bx+b；其中，所述建筑线斜率为a_b。这里，所述a_b的值为1.2。

其中，在蓝-近红外特征空间内作不透水层的样本散点图具体包括：在提取了三种地物的纯净像元数据样本的基础上，对纯净像元数据样本的蓝、近红外波段值进行导出，以蓝波段像元值为x轴横坐标，以近红外波段像元值为y轴纵坐标，将纯净像元数据样本散点进行作图，得到蓝-近红外特征空间内的三种地物的样本散点图。

具体地，纯净像元样本散点图如图4所示，在图4中，所述“●”代表亮色不透水层，所述代表蓝色不透水层，所述“○”代表暗色不透水层，所述“□”代表土壤，所述“△”代表阴影，所述“+”代表植被，所述“×”代表水体，所述“-----”代表建筑线，所述“——”代表土壤阴影线。

进一步地，在影像上选取纯净的土壤像元数据及山体影像数据构建土壤和山体阴影样本数据，对所述土壤和山体阴影样本数据进行最小二乘拟合，确定土壤阴影线方程y＝a_ssx+b；其中，所述土壤阴影线斜率为a_ss，所述a_ss的值为3.8。

利用公式(4)计算所述垂直不透水层指数PBI的系数n；其中，当n与a相等时，不透水层的提取效果最好。

$a=\frac{a_b+a_{\mathrm{ss}}}{2}---\left(4\right)$

其中，所述n值为2.5。

当确定出垂直不透水层指数PBI的系数n后，可以利用公式(6)计算PBI指数。

PBI＝nBlue-NIR  (6)

其中，所述Blue为辐亮度影像数据中的蓝光波段像元数据，所述NIR为辐亮度影像数据中的近红外波段像元数据。

根据均值标准差法确定所述PBI指数的不透水层提取阈值T；具体地，对辐亮度影像数据进行快速统计获取所述PBI影像像元的均值μ和标准差ε；

利用公式(7)确定所述PBI指数的不透水层提取阈值T。

T＝μ+ε  (7)

这里，如图5所示，所述μ值为0.074，所述ε值为0.03，当阈值T为均值μ与一倍标准差ε之和时，不透水层的提取效果最好。

当所述PBI指数大于所阈值T时，确定所述影像像元为不透水层区域，所述不透水层包括：亮色不透水层、暗色不透水层和蓝色不透水层。其中，有山体阴影的不透水层提取结果如图6所示。

当所述PBI指数不大于所阈值T时，确定所述影像像元为土壤、植被等透水层区域。

实施例四

实际应用中，本实施例是以武汉市作为无山体阴影的实验区进行不透水层的提取。其中，无山体阴影区域的原始影像图如图7所示。

当对无山体阴影的实验区进行不透水层的提取时，按照实施例一提供的方法，具体实施如下：

直接从云服务器或者本地服务器上获取武汉市的卫星影像数据，所述卫星影像数据可以包括：Landsat卫星影像数据。

根据公式(1)计算未校正的λ波段大气层辐亮度ρ′_λ；

ρ′_λ＝M_ρQ_cal-A_ρ  (1)

其中，在公式(1)中，所述M_ρ为增益系数；所述Q_cal为λ波段的影像像元亮度DN值；所述A_ρ为偏置系数；所述λ取值为自然数。

根据公式(2)对所述未校正的λ波段大气层光辐亮度ρ′_λ进行校正，获取校正后的λ波段大气层辐亮度ρ_λ；

${\mathrm{\ρ}}_{\mathrm{\λ}}=\frac{{\mathrm{\ρ}}_{\mathrm{\λ}}^{\′}}{\sin \left({\mathrm{\θ}}_{\mathrm{SE}}\right)}---\left(2\right)$

其中，在公式(2)中，所述θ_SE为当地的太阳高度角。

这里，当λ波段大气层辐亮度ρ_λ计算出以后，所述影像数据即被预处理为辐亮度影像数据。

根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI。

其中，所述调整归一化水体指数MNDWI是在归一化水体指数NDWI的基础上进行的计算，用于移除影像上的水体，减小后续的提取误差。

具体地，根据公式(3)计算所述调整归一化水体指数MNDWI；

$\mathrm{MNDWI}=\frac{\mathrm{Green}-\mathrm{NIR}}{\mathrm{Green}+\mathrm{NIR}}---\left(3\right)$

在公式(3)中，所述Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中的近红外波段像元数据。

进一步地，当所述调整归一化水体指数不大于0时，在所述辐亮度影像数据中提取纯净像元数据。

具体地，当所述MNDWI不大于0时，表明此影像区域为非水体区域，在影像中保留此区域。

当所述MNDWI大于0时，表明此影像区域为水体区域，对所述区域进行剔除。

进一步地，利用建筑线方程、土壤线方程确定垂直不透水层指数PBI的系数n；

具体地，在上一步保留的区域中通过目视解译选取各类典型地物，包括三种不透水层、土壤、植被、水体和阴影的纯净像元作为不透水层纯净像元数据样本；其中，所述纯净像元是指该像元中只包含一种地物，没有其他地物混合。所述纯净像元应尽量选取位于地物中心的像元作为纯净像元。

这里，因城市中不透水层类型可按其光谱特征的不同分为三种：亮色建筑、暗色建筑和蓝色建筑。其中，亮色和蓝色不透水层在蓝光波段的反射率较高，其余几种地物在蓝光波段反射率低，因此蓝光波段能够有效区分亮、蓝色建筑与裸土。而暗色建筑的光谱值在近红外波段无较大变化，裸土则显著升高，因此可利用蓝波段与近红外波段值的差异来区分暗色建筑与裸土。利用以上几种地物的光谱特性，可以确定垂直不透水层指数PBI。

这里，利用建筑线方程确定垂直不透水层指数PBI的系数n具体包括：

在蓝-近红外特征空间内作不透水层的样本散点图，对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y＝a_bx+b；其中，所述建筑线斜率为a_b。这里，所述a_b的值为1.4。

其中，在蓝-近红外特征空间内作不透水层的样本散点图具体包括：在提取了三种地物的纯净像元数据样本的基础上，对纯净像元数据样本的蓝、近红外波段值进行导出，以蓝波段像元值为x轴横坐标，以近红外波段像元值为y轴纵坐标，将纯净像元数据样本散点进行作图，得到蓝-近红外特征空间内的三种地物的样本散点图。

具体地，纯净像元样本散点图如图8所示，在图7中，所述“●”代表亮色不透水层，所述代表蓝色不透水层，所述“○”代表暗色不透水层，所述“□”代表土壤，所述“△”代表阴影，所述“+”代表植被，所述“×”代表水体，所述“-----”代表建筑线，所述“——”代表土壤线。

进一步地，在影像上选取纯净的土壤像元数据据构建土壤样本数据，对所述土壤样本数据进行最小二乘拟合，确定土壤线方程y＝a_sx+b；其中，所述土壤线斜率为a_s。其中，a_s的值为3.6。

利用公式(5)计算所述垂直不透水层指数PBI的系数n；其中，当n与a′相等时，不透水层的提取效果最好。

$a^{\′}=\frac{a_b+a_s}{2}---\left(5\right)$

其中，所述n值为2.5。

当确定出垂直不透水层指数PBI的系数n后，可以利用公式(6)计算PBI指数。

PBI＝nBlue-NIR  (6)

其中，所述Blue辐亮度影像数据中的蓝光波段像元数据，所述NIR为辐亮度影像数据中的近红外波段像元数据。

根据均值标准差法确定所述PBI指数的不透水层提取阈值T；具体地，对辐亮度影像数据进行快速统计获取所述PBI影像像元的均值μ和标准差ε；

利用公式(7)确定所述PBI指数的不透水层提取阈值T。

T＝μ+ε  (7)

这里，如图9所示，所述μ值为0.029，所述ε值为0.097，当阈值T为均值μ与一倍标准差ε之和时，不透水层的提取效果最好。

当所述PBI指数大于所阈值T时，确定所述影像像元为不透水层区域，所述不透水层包括：亮色不透水层、暗色不透水层和蓝色不透水层。其中，无山体阴影的不透水层提取结果如图10所示。

当所述PBI指数不大于所阈值T时，确定所述影像像元为土壤、植被等透水层区域。

以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种不透水层的提取方法，其特征在于，所述方法包括：

获取卫星影像数据；

计算所述影像数据在λ波段的大气层辐亮度ρ_λ，将所述影像数据预处理为辐亮度影像数据；

根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI；

当所述MNDWI不大于0时，在所述辐亮度影像数据中选取各类地物的纯净像元数据；

利用建筑线方程、土壤线方程及土壤阴影线方程确定垂直不透水层指数PBI的系数n；

利用公式PBI＝nBlue-NIR计算PBI指数；

根据均值标准差法确定所述PBI指数的不透水层的提取阈值T；其中，

当所述PBI指数大于所述阈值T时，确定影像像元为不透水层区域；所述λ取值为自然数。

2.如权利要求1所述的方法，其特征在于，所述计算所述影像数据在λ波段的大气层辐亮度ρ_λ包括：

根据公式ρ′_λ＝M_ρQ_cal-A_ρ计算未校正的λ波段大气层辐亮度ρ′_λ；

根据公式对所述未校正的λ波段大气层辐亮度ρ′_λ进行校正，获取校正后的λ波段的大气层辐亮度ρ_λ；其中，所述M_ρ为增益系数；所述Q_cal为λ波段的所述影像像元亮度DN值；所述A_ρ为偏置系数；所述θ_SE为当地的太阳高度角。

3.如权利要求1所述的方法，其特征在于，所述根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI包括：

利用公式计算所述调整归一化水体指数MNDWI；其中，所述Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中近红外波段像元数据。

4.如权利要求1所述的方法，其特征在于，所述利用建筑线方程确定垂直不透水层指数PBI的系数n包括：

对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y＝a_bx+b；

当确定所述不透水层所在区域存在山体阴影时，构建土壤和山体阴影样本数据，对所述土壤和山体阴影样本数据进行最小二乘拟合，确定土壤阴影线方程y＝a_ssx+b；

利用公式计算所述垂直不透水层指数PBI的系数n；其中，n与a相等。

5.如权利要求4所述的方法，其特征在于，当确定所述不透水层所在区域不存在山体阴影时，所述方法还包括：

构建土壤样本数据，对所述土壤样本数据进行最小二乘拟合，确定土壤线方程y＝a_sx+b；

利用公式计算所述垂直不透水层指数PBI的系数n；其中，n与a′相等。

6.如权利要求1所述的方法，其特征在于，所述根据均值标准差法确定所述PBI指数的不透水层提取阈值T包括：

对辐亮度影像数据进行统计获取PBI影像像元的均值μ和标准差ε；

利用公式T＝μ+ε确定所述PBI指数的不透水层提取阈值T。

7.一种不透水层的提取装置，其特征在于，所述装置包括：

获取单元，所述获取单元用于获取卫星影像数据；

预处理单元，所述预处理单元用于计算所述影像数据在λ波段的大气层光谱反射率ρ_λ，将所述影像数据预处理为辐亮度影像数据；

计算单元，所述计算单元用于根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI；

选取单元，所述选取单元用于确定当所述调整归一化水体指数不大于0时，在所述辐亮度影像数据中选取各类地物的纯净像元数据；

第一确定单元，所述第一确定单元用于利用建筑线方程、土壤线方程及土壤阴影线方程确定垂直不透水层指数PBI的系数n；

第二确定单元，所述第二确定单元用于利用公式PBI＝nBlue-NIR计算PBI指数；

第三确定单元，所述第三确定单元用于根据均值标准差法确定所述PBI指数的不透水层的提取阈值T；

第四确定单元，所述第四确定单元用于确定当所述PBI指数大于所述阈值T时，确定影像像元为不透水层区域；其中，所述λ取值为自然数。

8.如权利要求7所述的装置，其特征在于，所述预处理单元具体用于：

根据公式ρ′_λ＝M_ρQ_cal-A_ρ计算未校正的λ波段大气层辐亮度ρ′_λ；

根据公式对所述未校正的λ波段大气层辐亮度ρ′_λ进行校正，获取校正后的λ波段的大气层辐亮度ρ_λ；其中，所述M_ρ为增益系数；所述Q_cal为λ波段的DN值；所述A_ρ为偏置系数；所述θ_SE为当地的太阳高度角。

9.如权利要求7所述的装置，其特征在于，所述计算单元具体用于：

利用公式计算所述调整归一化水体指数MNDWI；其中，所述Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中的近红外波段像元数据。

10.如权利要求7所述的装置，其特征在于，所述第一确定单元具体用于：

对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y＝a_bx+b；

当确定所述不透水层所在区域存在山体阴影时，构建土壤和山体阴影样本数据，对所述土壤和山体阴影样本数据进行最小二乘拟合，确定土壤阴影线方程y＝a_ssx+b；

利用公式计算所述垂直不透水层指数PBI的系数n；其中，n与a相等；

当确定所述不透水层所在区域不存在山体阴影时，构建土壤样本数据，对所述土壤样本数据进行最小二乘拟合，确定土壤线方程y＝a_sx+b；

利用公式计算所述垂直不透水层指数PBI的系数n；其中，n与a′相等。