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

一种不透水层的提取方法及装置的制造方法
【技术领域】
[0001] 本发明属于遥感影像处理技术领域，尤其设及一种不透水层的提取方法及装置。
【背景技术】
[0002] 近年来，城市作为人口最为集聚和人类活动最为强烈的地区，其环境和气候等日 益受到关注；城市的扩展直接改变了±地覆盖、±地利用的状况，也直接影响局部的气候、 生态和环境。不透水层作为城市化进程和城市化率评估的重要指标，对分析评价城市地区 的气候、环境、水循环等具有重要的意义。
[0003] 现有技术中，利用遥感影像进行的不透水层提取方法主要包括；监督分类法、混合 像元分解法、谱间关系法、不透水指数法等；但该些方法在进行大数据的处理过程中主要依 靠人工提取，不仅耗时耗力，而且人为主观影响严重。实际应用中，普遍采用的方法是指数 提取法，但指数提取法利用已有的不透水层指数在提取不透水层时无法将裸±部分剔除， 导致结果不精确，在后续的提取过程中产生较大误差；特别是在较干旱、裸±地表较多的西 北地区，W及东部地区旱季影像上的提取中产生更大的误差。
[0004] 基于此，目前亟需一种新型的提取方法，重新构建能够有效区分裸±的不透水层 提取指数，降低不透水层的提取误差。

【发明内容】

[0005] 针对现有技术存在的问题，本发明实施例提供了一种不透水层的提取方法及装 置。
[0006] 本发明提供一种不透水层的提取方法，所述方法包括：
[0007] 获取卫星影像数据；
[000引计算所述影像数据在A波段的大气层福亮度P 将所述影像数据预处理为福亮 度影像数据；
[0009] 根据所述福亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调 整归一化水体指数MNDWI;
[0010] 当所述MNDWI不大于0时，在所述福亮度影像数据中选取各类地物的纯净像元数 据；
[0011] 利用建筑线方程、±壤线方程及±壤阴影线方程确定垂直不透水层指数PBI的系 数n;
[001引利用公式PBI=nBlue-NIR计算PBI指数；
[0013] 根据均值标准差法确定所述PBI指数的不透水层的提取阔值T;其中，
[0014] 当所述PBI指数大于所述阔值T时，确定影像像元为不透水层区域；所述A取值 为自然数。
[0015] 上述方案中，所述计算所述影像数据在A波段的大气层福亮度P,包括；
[0016] 根据公式P' .=MpQt,i-Ap计算未校正的A波段大气层福亮度P'
[0017] 根据公式
对所述未校正的A波段大气层福亮度P 进行校正， 获取校正后的A波段的大气层福亮度P其中，所述MP为增益系数；所述为A波段 的所述影像像元亮度DN值；所述Ap为偏置系数；所述0SC为当地的太阳高度角。
[0018] 上述方案中，所述根据所述福亮度影像数据中的绿光波段像元数据和近红外波段 像元数据计算调整归一化水体指数MNDWI包括：
[0019] 利用公式
计算所述调整归一化水体指数MNDWI;其中，所 述Green为所述福亮度影像数据中的绿光波段像元数据，所述NIR为所述福亮度影像数据 中近红外波段像元数据。
[0020] 上述方案中，所述利用建筑线方程确定垂直不透水层指数PBI的系数n包括：
[0021] 对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y= BfaX+b；
[0022] 当确定所述不透水层所在区域存在山体阴影时，构建±壤和山体阴影样本数据， 对所述±壤和山体阴影样本数据进行最小二乘拟合，确定±壤阴影线方程y=a,,x+b;
[0023] 利用公式
计算所述垂直不透水层指数PBI的系数n;其中，n与a相 等。
[0024] 上述方案中，当确定所述不透水层所在区域不存在山体阴影时，所述方法还包 括：
[0025] 构建±壤样本数据，对所述±壤样本数据进行最小二乘拟合，确定±壤线方程y =a曰x+b;
[0026] 利用公式
计算所述垂直不透水层指数PBI的系数n;其中，n与a' 相等。
[0027] 上述方案中，所述根据均值标准差法确定所述PBI指数的不透水层提取阔值T包 括：
[002引对福亮度影像数据进行统计获取PBI影像像元的均值y和标准差e;
[0029] 利用公式T=y+e确定所述PBI指数的不透水层提取阔值T。
[0030] 本发明同时还提供了一种不透水层的提取装置，所述装置包括：
[0031] 获取单元，所述获取单元用于获取卫星影像数据；
[0032] 预处理单元，所述预处理单元用于计算所述影像数据在A波段的大气层光谱反 射率P将所述影像数据预处理为福亮度影像数据；
[0033] 计算单元，所述计算单元用于根据所述福亮度影像数据中的绿光波段像元数据和 近红外波段像元数据计算调整归一化水体指数MNDWI;
[0034] 选取单元，所述选取单元用于确定当所述调整归一化水体指数不大于0时，在所 述福亮度影像数据中选取各类地物的纯净像元数据；
[0035] 第一确定单元，所述第一确定单元用于利用建筑线方程、±壤线方程及±壤阴影 线方程确定垂直不透水层指数PBI的系数n;
[0036] 第二确定单元，所述第二确定单元用于利用公式PBI=nBlue-NIR计算PBI指数；
[0037] 第=确定单元，所述第=确定单元用于根据均值标准差法确定所述PBI指数的不 透水层的提取阔值T;
[003引第四确定单元，所述第四确定单元用于确定当所述PBI指数大于所述阔值T时，确 定影像像元为不透水层区域；其中，所述A取值为自然数。
[0039] 上述方案中，所述预处理单元具体用于：
[0040] 根据公式P' .=MpQt,i-Ap计算未校正的A波段大气层福亮度P'
[0041] 根据公式，
对所述未校正的A波段大气层福亮度P进行校正， 获取校正后的A波段的大气层福亮度P其中，所述MP为增益系数；所述为A波段 的DN值；所述Ap为偏置系数；所述0SE为当地的太阳高度角。
[0042] 上述方案中，所述计算单元具体用于：
[0043] 利用公式
计算所述调整归一化水体指数MNDWI;其中，所 述Green为所述福亮度影像数据中的绿光波段像元数据，所述NIR为所述福亮度影像数据 中的近红外波段像元数据。
[0044] 上述方案中，所述第一确定单元具体用于：
[0045] 对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y= BfaX+b；
[0046] 当确定所述不透水层所在区域存在山体阴影时，构建±壤和山体阴影样本数据， 对所述±壤和山体阴影样本数据进行最小二乘拟合，确定±壤阴影线方程y=a,,x+b;
[0047] 利用公式
计算所述垂直不透水层指数PBI的系数n;其中，n与a相 等；
[0048] 当确定所述不透水层所在区域不存在山体阴影时，构建±壤样本数据，对所述± 壤样本数据进行最小二乘拟合，确定±壤线方程y=a,x+b;
[0049] 利用公式
计算所述垂直不透水层指数PBI的系数n;其中，n与a' 相等。
[0050] 本发明提供了一种不透水层的提取方法及装置，所述方法包括：获取卫星影像数 据；计算所述影像数据在A波段的大气层福亮度P将所述影像数据预处理为福亮度影 像数据；根据所述福亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整 归一化水体指数MNDWI，当所述MNDWI不大于0时，在所述福亮度影像数据中选取各类地物 的纯净像元数据；利用建筑线方程、±壤线方程及±壤阴影线方程确定垂直不透水层指数 PBI的系数n;利用公式PBI=nBlue-NIR计算PBI指数；根据均值标准差法确定所述PBI 指数的不透水层的提取阔值T;其中，当所述PBI指数大于所阔值T时，确定影像像元为不 透水层区域；所述A取值为自然数；如此，所述不透水层的提取方法构建的PBI指数能够 有效区分各类不透水层、水体、植被及裸±地表，降低了不透水层的提取误差；并且该方法 对不同数据源的卫星影像数据均具有适用性。
【附图说明】
[0051] 图1为本发明实施例一提供的不透水层的提取方法流程示意图；
[0化2] 图2为本发明实施例二提供的不透水层的提取装置结构示意图；
[0053]图3为本发明实施例S提供的有山体阴影区域的原始影像图； [0化4] 图4为本发明实施例=提供的有山体阴影的纯净像元样本散点示意图；
[0055] 图5为本发明实施例S提供的有山体阴影的PBI指数的图像直方图；
[0056] 图6为本发明实施例S提供的有山体阴影区域的不透水层提取结果示意图；
[0057] 图7为本发明实施例四提供的无山体阴影区域的原始影像图；
[005引图8为本发明实施例四提供的无山体阴影的纯净像元样本散点示意图；
[0化9] 图9为本发明实施例四提供的无山体阴影的PBI指数的图像直方图；
[0060] 图10为本发明实施例四提供的无山体阴影的不透水层提取结果示意图。
【具体实施方式】
[0061] 为了能够在不透水层的提取中有效区分各类不透水层、水体、植被及裸±地表，降 低了不透水层的提取误差，本发明提供了一种不透水层的提取方法及装置，所述方法包括： 获取卫星影像数据；计算所述影像数据在A波段的大气层福亮度P 将所述影像数据预 处理为福亮度影像数据；根据所述福亮度影像数据中的绿光波段像元数据和近红外波段像 元数据计算调整归一化水体指数MNDWI;当所述MNDWI不大于0时，在所述福亮度影像数据 中选取各类地物纯净像元数据；利用建筑线方程、上壤线方程及上壤阴影线方程确定垂直 不透水层指数PBI的系数n;利用公式PBI=nBlue-NIR计算PBI指数；根据均值标准差法 确定所述PBI指数的不透水层的提取阔值T;其中，当所述PBI指数大于所述阔值T时，确 定影像像元为不透水层区域；其中，所述A取值为自然数。
[0062] 下面通过附图及具体实施例对本发明的技术方案做进一步的详细说明。
[0063] 实施例一
[0064] 本实施例提供一种不透水层的提取方法，如图1所示，所述方法主要包括W下步 骤：
[00化]步骤110,获取卫星影像数据。
[0066] 本步骤中，可W直接从云服务器或者本地服务器上获取卫星影像数据，所述卫星 影像数据可W包括；Landsat卫星影像数据、Aster卫星影像数据等。
[0067] 步骤111，计算所述影像数据在A波段的大气层福亮度P将所述影像数据预处 理为福亮度影像数据。
[0068] 本步骤中，根据公式（1)计算未校正A波段的大气层福亮度P'
[00 例 P' ,=MpQeal-Ap(1)
[0070] 其中，在公式（1)中，所述Mp为增益系数；所述Qcal为A波段的影像像元亮度DN 值；所述Ap为偏置系数，所述A取值为自然数。
[0071] 根据公式（2)对所述未校正的A波段大气层光福亮度P进行校正，获取校 正后的A波段大气层福亮度P\ ;
[00 巧 （2)
[0073] 其中，在公式（2)中，所述0SE为当地的太阳高度角。
[0074] 该里，当A波段大气层福亮度P,计算出W后，所述影像数据即被预处理为福亮 度影像数据。
[0075] 步骤112,根据所述福亮度影像数据中的绿光波段像元数据和近红外波段像元数 据计算调整归一化水体指数MNDWI。
[0076]本步骤中，所述调整归一化水体指数（MNDWI，ModifiedNormalizedDifference WaterIndex)是在归一化水体指数（NDWI，NormalizedDifferenceWaterIndex)的基础 上进行的计算，用于移除影像上的水体，减小后续的提取误差。
[0077] 具体地，根据公式做计算所述调整归一化水体指数MNDWI ;
[0078]
(3)
[0079] 其中，在公式（3)中，所述Green为所述福亮度影像数据中的绿光波段像元数据， 所述NIR为所述福亮度影像数据中的近红外波段像元数据。
[0080] 步骤113,当所述调整归一化水体指数不大于0时，在所述福亮度影像数据中提取 纯净像元数据。
[0081] 本步骤中，当所述MNDWI不大于0时，表明此影像区域为非水体区域，在影像中保 留此区域，在此影像区域中通过目视解译选取各类典型地物，包括=种不透水层、±壤、植 被、水体和阴影的纯净像元作为不透水层纯净像元数据样本；其中，所述纯净像元是指该像 元中只包含一种地物，没有其他地物混合。所述纯净像元应尽量选取位于地物中屯、的像元 作为纯净像元。
[0082] 当所述MNDWI大于0时，表明此影像区域为水体区域，对所述区域进行剔除。
[0083] 步骤114,利用建筑线方程、±壤线方程和±壤阴影线方程确定垂直不透水层指数 PBI的系数n。
[0084] 本步骤中，因城市中不透水层类型可按其光谱特征的不同分为=种：亮色建筑、暗 色建筑和藍色建筑。其中，亮色和藍色不透水层在藍光波段的反射率较高，其余几种地物在 藍光波段反射率低，因此藍光波段能够有效区分亮、藍色建筑与裸±。而暗色建筑的光谱值 在近红外波段无较大变化，裸±则显著升高，因此可利用藍波段与近红外波段值的差异来 区分暗色建筑与裸±。利用W上几种地物的光谱特性，可W确定垂直不透水层指数PBI。
[0085] 该里，利用建筑线方程确定垂直不透水层指数PBI的系数n具体包括；在藍-近红 外特征空间内作不透水层的样本散点图，对所述不透水层的样本散点进行线性最小二乘拟 合，确定所述建筑线方程y= 3bX+b;其中，所述建筑线斜率为3b。
[0086] 其中，在藍-近红外特征空间内作不透水层的样本散点图具体包括；在提取了 = 种不透水层的纯净像元数据样本的基础上，对纯净像元数据样本的藍、近红外波段值进行 导出，W藍波段像元值为X轴横坐标，W近红外波段像元值为y轴纵坐标，将纯净像元数据 样本散点进行作图，得到藍-近红外特征空间内的=种不透水层的样本散点图。
[0087] 进一步地，判断该影像区域所在的地区是否有大型山体及山体阴影，当确定所述 不透水层所在区域存在山体阴影时，在影像上选取纯净的±壤像元数据及山体影像数据构 建±壤和山体阴影样本数据，对所述±壤和山体阴影样本数据进行最小二乘拟合，确定± 壤阴影线方程y=屯,x+b;其中，所述±壤阴影线斜率为a,,。
[008引利用公式（4)计算所述垂直不透水层指数PBI的系数n;其中，当n与a相等时， 不透水层的提取效果最好。
[0089]
(4)
[0090] 进一步地，当确定所述不透水层所在区域不存在山体阴影时，在影像上选取纯净 的±壤像元数据据构建±壤样本数据，对所述±壤样本数据进行最小二乘拟合，确定±壤 线方程y=a,x+b;其中，所述±壤线斜率为a,。
[OOW] 利用公式（5)计算所述垂直不透水层指数PBI的系数n;其中，当n与a'相等时， 不透水层的提取效果最好。
[009引
（5)
[009引步骤115,利用公式PBI=nBlue-NIR计算PBI指数。
[0094] 本步骤中，当确定出垂直不透水层指数PBI的系数n后，可W利用公式（6)计算 roi指数。
[0095] PBI = nBlue-NIR (6)
[0096] 其中，所述Blue为福亮度影像数据中的藍光波段像元数据，所述NIR为福亮度影 像数据中的近红外波段像元数据。
[0097] 步骤116,根据均值标准差法确定所述PBI指数的不透水层提取阔值T。
[009引本步骤中，可W根据均值标准差法确定所述PBI指数的不透水层提取阔值T ;具体 地，对福亮度影像数据进行快速统计获取所述PBI影像像元的均值y和标准差e ;
[0099] 利用公式（7)确定所述PBI指数的不透水层提取阔值T。
[0100] T=y+e(7)
[0101] 该里，当阔值T为均值y与一倍标准差e之和时，不透水层的提取效果最好。
[0102] 步骤117,根据所述不透水层提取阔值T对不透水层进行提取。
[0103] 本步骤中，当所述PBI指数大于所阔值T时，确定所述影像像元为不透水层区域， 所述不透水层包括；亮色不透水层、暗色不透水层和藍色不透水层。
[0104] 当所述PBI指数不大于所阔值T时，确定所述影像像元为±壤、植被等透水层区 域。
[01化]本实施例提供的不透水层提取方法，重新构建能够有效区分裸±的不透水层提取 指数，除了能够区分不透水层与水体、植被地物，还能够有效区分各类不透水层与裸±及阴 影，降低不透水层的提取误差，相比已有的提取方法，在较干旱、裸±地表较多的西北地区， W及东部地区旱季影像上的提取效果都有很大的提高。且对于不同的影像数据源均具有适 用性，因此，本实 施例提供的不透水层的提取方法有较广泛的应用前景，较强的实用性和推 广性。
[0106] 实施例二
[0107] 相对应实施例一，本实施例还提供了一种不透水层的提取装置，如图2所示，所述 装置包括：获取单元21、预处理单元22、计算单元23、选取单元24、第一确定单元25、第二 确定单元26,第=确定单元27、第四确定单元28 ;其中，
[0108] 所述获取单元21用于获取卫星影像数据；具体地，所述获取单元21可W直接从云 服务器或者本地服务器上获取卫星影像数据，所述卫星影像数据可W包括；Landsat卫星 影像数据、Aster卫星影像数据等。
[0109] 当所述获取单元21获取到卫星影像数据后，所述预处理单元22用于计算所述影 像数据在A波段的大气层光谱反射率P将所述影像数据进行预处理为福亮度影像数 据。具体地，所述预处理单元22根据公式（1)计算未校正A波段的大气层福亮度P'
[0110] P' ,=MpQcai_Ap(1)
[011U其中，在公式（1)中，所述Mp为增益系数；所述Qtai为A波段的影像像元亮度DN值；所述Ap为偏置系数，所述A取值为自然数。。
[0112] 根据公式（2)对所述未校正的A波段大气层光福亮度P进行校正，获取校 正后的A波段大气层福亮度P\ ;
[0113]
(2)
[0114] 其中，在公式（2)中，所述0SE为当地的太阳高度角。
[0115] 该里，当A波段的大气层福亮度P,计算出W后，所述影像数据即被预处理为福 亮度影像数据。
[0116] 当所述预处理单元22将影像数据处理为福亮度影像数据后，所述计算单元23用 于根据所述福亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归一 化水体指数MNDWI。
[0117] 该里，所述调整归一化水体指数MNDWI是在归一化水体指数NDWI的基础上进行的 计算，用于移除影像上的水体，减小后续的提取误差。
[0118] 具体地，所述计算单元23根据公式（3)计算所述调整归一化水体指数MNDWI;
[0119]
(3)
[0120] 其中，在公式（3)中所述Green为所述福亮度影像数据中的绿光波段像元数据，所 述NIR为所述福亮度影像数据中的近红外波段像元数据。
[0121] 当所述计算单元23计算出所述调整归一化水体指数MNDWI后，所述选取单元24 用于判断所述调整归一化水体指数MNDWI是否大于0,当确定所述调整归一化水体指数 MNDWI大于0时，表明此影像区域为非水体区域，在影像中保留此区域，通过目视解译选取 此影像区域中各类典型地物，包括=种不透水层、±壤、植被、水体和阴影的纯净像元作为 不透水层纯净像元数据样本；其中，所述纯净像元是指该像元中只包含一种地物，没有其他 地物混合。所述纯净像元应尽量选取位于地物中屯、的像元作为纯净像元。
[0122] 当所述选取单元24确定MNDWI大于0时，表明此影像区域为水体区域，对所述区 域进行剔除。
[0123] 当所述选取单元24确定述调整归一化水体指数MNDWI大于0时，所述第一确定单 元25利用建筑线方程确定垂直不透水层指数PBI的系数n。
[0124] 具体地，所述第一确定单元25在此影像区域中确定S种不透水层地物的纯净像 元作为不透水层纯净像元数据；在藍-近红外特征空间内作不透水层的样本散点图，对所 述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y=abX+b;其中，所 述建筑线斜率为3b。
[0125] 当所述第一确定单元25确定出建筑线方程后，判断该影像区域所在的地区是否 有大型山体及山体阴影，当确定所述不透水层所在区域存在山体阴影时，在影像上选取纯 净的±壤像元数据及山体影像数据构建±壤和山体阴影样本数据，对所述±壤和山体阴影 样本数据进行最小二乘拟合，确定±壤阴影线方程y=a,,x+b;其中，所述±壤阴影线斜率 为日55。
[01%] 当所述第一确定单元25确定出±壤阴影线方程后，利用公式（4)计算所述垂直不 透水层指数PBI的系数n;其中，当n与a相等时，不透水层的提取效果最好。
[0127]
(4)
[0128] 进一步地，当所述第一确定单元25确定所述不透水层所在区域不存在山体阴影 时，在影像上选取纯净的±壤像元数据据构建±壤样本数据，对所述±壤样本数据进行最 小二乘拟合，确定上壤线方程y=a,x+b;其中，所述±壤线斜率为a,。
[0129] 所述第一确定单元25确定出±壤线方程后，利用公式（5)计算所述垂直不透水层 指数PBI的系数n;其中，当n与a'相等时，不透水层的提取效果最好。
[0130]
(5)
[0131] 当所述第一确定单元25确定出垂直不透水层指数PBI的系数n后，所述第二确定 单元26用于根据公式做计算PBI指数。
[0132] PBI = nBlue-NIR (6)
[0133] 其中，所述Blue为福亮度影像数据中的藍光波段像元数据，所述NIR为福亮度影 像数据中的近红外波段像元数据。
[0134] 该里，因城市中不透水层类型可按其光谱特征的不同分为=种；亮色建筑、暗色建 筑和藍色建筑。其中，亮色和藍色不透水层在藍光波段的反射率较高，其余几种地物在藍光 波段反射率低，因此藍光波段能够有效区分亮、藍色建筑与裸±。而暗色建筑的光谱值在近 红外波段无较大变化，裸上则显著升高，因此，所述第二确定单元26可利用藍波段与近红 外波段值的差异来区分暗色建筑与裸± ;利用W上几种地物的光谱特性，可W确定PBI指 数。
[01巧]当所述第二确定单元26确定出PBI指数后，所述第=确定单元27用于；根据均值 标准差法确定所述PBI指数的不透水层提取阔值T。具体地，所述第=确定单元27对福亮 度影像数据进行快速统计获取所述PBI影像像元的均值y和标准差e;
[0136] 利用公式（7)确定所述PBI指数的不透水层提取阔值T。
[0137]T=y+e(7)
[0138] 该里，当阔值T为均值y与一倍标准差e之和时，不透水层的提取效果最好。
[0139] 当所述第=确定单元27确定出所述阔值T后，所述第四确定单元28具体用于：判 断所述PBI指数是否大于所阔值T，当所述第四确定单元28所述PBI指数大于所阔值T时， 确定所述影像像元为不透水层区域，所述不透水层包括；亮色不透水层、暗色不透水层和藍 色不透水层。
[0140] 当所述第四确定单元28确定所述PBI指数不大于所阔值T时，确定所述影像像元 为±壤、植被等透水层区域。
[0141] 实际应用中，所述获取单元21、所述预处理单元22、所述计算单元23、所述选取 单元24、所述第一确定单元25、所述第二确定单元26,所述第=确定单元27、所述第四确 定单元28可由该装置中的中央处理器（CPU,CentralProcessing化it)、数字信号处理 器值SP，DigtalSi即alProcessor)、可编程逻辑阵列（FPGA，FieldProgramm油leGate Array)、微控制单元（MCU，MicroControllerUnit)实现。
[0142] 实施例S
[0143] 实际应用中，本实施例是W北京市作为有山体阴影的实验区进行不透水层的提 取，所述有山体阴影区域的原始影像图如图3所示。
[0144] 当对北京有山体阴影的实验区进行不透水层的提取时，按照实施例一提供的方 法，具体实施如下；
[0145] 直接从云服务器或者本地服务器上获取北京市的卫星影像数据，所述卫星影像数 据可W包括；Landsat卫星影像数据。
[0146] 根据公式（1)计算未校正的A波段大气层福亮度P'
[0147] P' ,=MpQcai-Ap(1)
[0148] 其中，在公式（1)中，所述Mp为增益系数；所述Qtai为入波段的影像像元亮度DN 值；所述Ap为偏置系数；所述A取值为自然数。
[0149] 根据公式（2)对所述未校正的A波段大气层光福亮度P进行校正，获取校 正后的A波段大气层福亮度P\ ;
[0150]
(2；
[0151] 其中，在公式（2)中，所述0SE为当地的太阳高度角。
[0152] 该里，当A波段大气层福亮度P,计算出W后，所述影像数据即被预处理为福亮 度影像数据。
[0153] 根据所述福亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调 整归一化水体指数 MNDWI。
[0154] 其中，所述调整归一化水体指数MNDWI是在归一化水体指数NDWI的基础上进行的 计算，用于移除影像上的水体，减小后续的提取误差。
[0巧5] 具体地，根据公式（3)计算所述调整归一化水体指数MNDWI;
[0156]
(3)
[0157] 在公式（3)中，所述Green为所述福亮度影像数据中的绿光波段像元数据，所述 NIR为所述福亮度影像数据中的近红外波段像元数据。
[0158] 进一步地，当所述调整归一化水体指数不大于0时，在所述福亮度影像数据中提 取纯净像元数据。
[0159] 具体地，当所述MNDWI不大于0时，表明此影像区域为非水体区域，在影像中保留 此区域。
[0160]当所述MNDWI大于0时，表明此影像区域为水体区域，对所述区域进行剔除。
[01W]进一步地，利用建筑线方程、±壤阴影线方程确定垂直不透水层指数PBI的系数n；
[0162] 具体地，在上一步保留的区域中通过目视解译选取各类典型地物，包括S种不透 水层、±壤、植被、水体和阴影的纯净像元作为不透水层纯净像元数据样本；其中，所述纯净 像元是指该像元中只包含一种地物，没有其他地物混合。所述纯净像元应尽量选取位于地 物中屯、的像元作为纯净像元。
[0163] 该里，因城市中不透水层类型可按其光谱特征的不同分为=种；亮色建筑、暗色建 筑和藍色建筑。其中，亮色和藍色不透水层在藍光波段的反射率较高，其余几种地物在藍光 波段反射率低，因此藍光波段能够有效区分亮、藍色建筑与裸±。而暗色建筑的光谱值在近 红外波段无较大变化，裸±则显著升高，因此可利用藍波段与近红外波段值的差异来区分 暗色建筑与裸±。利用W上几种地物的光谱特性，可W确定垂直不透水层指数PBI。
[0164] 该里，利用建筑线方程确定垂直不透水层指数PBI的系数n具体包括：
[0165] 在藍-近红外特征空间内作不透水层的样本散点图，对所述不透水层的样本散点 进行线性最小二乘拟合，确定所述建筑线方程y=abX+b;其中，所述建筑线斜率为ab。该 里，所述却的值为1.2。
[0166] 其中，在藍-近红外特征空间内作不透水层的样本散点图具体包括；在提取了S 种地物的纯净像元数据样本的基础上，对纯净像元数据样本的藍、近红外波段值进行导出， W藍波段像元值为X轴横坐标，W近红外波段像元值为y轴纵坐标，将纯净像元数据样本散 点进行作图，得到藍-近红外特征空间内的S种地物的样本散点图。
[0167] 具体地，纯净像元样本散点图如图4所示，在图4中，所述"?"代表亮色不透水层， 所述"?"代表藍色不透水层，所述"0"代表暗色不透水层，所述"□"代表±壤，所述"A" 代表阴影，所述代表植被，所述"X"代表水体，所述"--"代表建筑线，所述"一一'' 代表±壤阴影线。
[0168] 进一步地，在影像上选取纯净的±壤像元数据及山体影像数据构建±壤和山体阴 影样本数据，对所述±壤和山体阴影样本数据进行最小二乘拟合，确定±壤阴影线方程y =assX+b;其中，所述±壤阴影线斜率为心所述aj勺值为3. 8。
[0169] 利用公式（4)计算所述垂直不透水层指数PBI的系数n;其中，当n与a相等时， 不透水层的提取效果最好。
[0170] a =ri', :d、、(4) 厶
[0171] 其中，所述n值为2. 5。
[0172] 当确定出垂直不透水层指数PBI的系数n后，可W利用公式（6)计算PBI指数。 [017引PBI=nBlue-NIR(6)
[0174]其中，所述Blue为福亮度影像数据中的藍光波段像元数据，所述NIR为福亮度影 像数据中的近红外波段像元数据。
[01巧]根据均值标准差法确定所述PBI指数的不透水层提取阔值T;具体地，对福亮度影 像数据进行快速统计获取所述PBI影像像元的均值y和标准差e;
[0176]利用公式（7)确定所述PBI指数的不透水层提取阔值T。
[0177]T=y+e(7)
[0178] 该里，如图5所示，所述y值为0. 074,所述e值为0. 03,当阔值T为均值y与 一倍标准差e之和时，不透水层的提取效果最好。
[0179] 当所述PBI指数大于所阔值T时，确定所述影像像元为不透水层区域，所述不透水 层包括；亮色不透水层、暗色不透水层和藍色不透水层。其中，有山体阴影的不透水层提取 结果如图6所示。
[0180] 当所述PBI指数不大于所阔值T时，确定所述影像像元为±壤、植被等透水层区 域。
[0181] 实施例四
[0182] 实际应用中，本实施例是W武汉市作为无山体阴影的实验区进行不透水层的提 取。其中，无山体阴影区域的原始影像图如图7所示。
[0183] 当对无山体阴影的实验区进行不透水层的提取时，按照实施例一提供的方法，具 体实施如下：
[0184] 直接从云服务器或者本地服务器上获取武汉市的卫星影像数据，所述卫星影像数 据可W包括；Landsat卫星影像数据。
[0化5] 根据公式（1)计算未校正的A波段大气层福亮度P'
[0186] P' ,=MpQcai-Ap(1)
[0187] 其中，在公式（1)中，所述Mp为增益系数；所述Qcai为A波段的影像像元亮度DN 值；所述Ap为偏置系数；所述A取值为自然数。
[0188] 根据公式（2)对所述未校正的A波段大气层光福亮度P 进行校正，获取校 正后的A波段大气层福亮度P\ ;
[0189]
C2)
[0190] 其中，在公式（2)中，所述0SE为当地的太阳高度角。
[0191] 该里，当A波段大气层福亮度P,计算出W后，所述影像数据即被预处理为福亮 度影像数据。
[0192] 根据所述福亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调 整归一化水体指数MNDWI。
[0193] 其中，所述调整归一化水体指数MNDWI是在归一化水体指数NDWI的基础上进行的 计算，用于移除影像上的水体，减小后续的提取误差。
[0194] 具体地，根据公式（3)计算所述调整归一化水体指数MNDWI;
[0195]
(3)
[0196] 在公式（3)中，所述Green为所述福亮度影像数据中的绿光波段像元数据，所述 NIR为所述福亮度影像数据中的近红外波段像元数据。
[0197] 进一步地，当所述调整归一化水体指数不大于0时，在所述福亮度影像数据中提 取纯净像元数据。
[019引具体地，当所述MNDWI不大于0时，表明此影像区域为非水体区域，在影像中保留 此区域。
[0199] 当所述MNDWI大于0时，表明此影像区域为水体区域，对所述区域进行剔除。
[0200] 进一步地，利用建筑线方程、±壤线方程确定垂直不透水层指数PBI的系数n;
[0201] 具体地，在上一步保留的区域中通过目视解译选取各类典型地物，包括=种不透 水层、±壤、植被、水体和阴影的纯净像元作为不透水层纯净像元数据样本；其中，所述纯净 像元是指该像元中只包含一种地物，没有其他地物混合。所述纯净像元应尽量选取位于地 物中屯、的像元作为纯净像元。
[0202] 该里，因城市中不透水层类型可按其光谱特征的不同分为=种；亮色建筑、暗色建 筑和藍色建筑。其中，亮色和藍色不透水层在藍光波段的反射率较高，其余几种地物在藍光 波段反射率低，因此藍光波段能够有效区分亮、藍色建筑与裸±。而暗色建筑的光谱值在近 红外波段无较大变化，裸±则显著升高，因此可利用藍波段与近红外波段值的差异来区分 暗色建筑与裸±。利用W上几种地物的光谱特性，可W确定垂直不透水层指数PBI。
[0203] 该里，利用建筑线方程确定垂直不透水层指数PBI的系数n具体包括：
[0204] 在藍-近红外特征空间内作不透水层的样本散点图，对所述不透水层的样本散点 进行线性最小二乘拟合，确定所述建筑线方程y=abX+b;其中，所述建筑线斜率为ab。该 里，所述却的值为1.4。
[0205] 其中，在藍-近红外特征空间内作不透水层的样本散点图具体包括；在提取了 = 种地物的纯净像元数据样本的基础上，对纯净像元数据样本的藍、近红外波段值进行导出， W藍波段像元值为X轴横坐标，W近红外波段像元值为y轴纵坐标，将纯净像元数据样本散 点进行作图，得到藍-近红外特征空间内的S种地物的样本散点图。
[0206] 具体地，纯净像元样本散点图如图8所示，在图7中，所述"?"代表亮色不透水层， 所述"?"代表藍色不透水层，所述"0"代表暗色不透水层，所述"□"代表±壤，所述"A" 代表阴影，所述"+"代表植被，所述"X"代表水体，所述"--"代表建筑线，所述"一一，' 代表±壤线。
[0207] 进一步地，在影像上选取纯净的±壤像元数据据构建±壤样本数据，对所述±壤 样本数据进行最小二乘拟合，确定±壤线方程y=a,x+b;其中，所述±壤线斜率为a,。其 中，3s的值为3. 6。
[0208] 利用公式（5)计算所述垂直不透水层指数PBI的系数n;其中，当n与a'相等时， 不透水层的提取效果最好。
[0209]
(5)
[0210] 其中，所述n值为2. 5。
[0211] 当确定出垂直不透水层指数PBI的系数n后，可W利用公式做计算PBI指数。 [021 引PBI=nBlue-NIR(6)
[0213] 其中，所述Blue福亮度影像数据中的藍光波段像元数据，所述NIR为福亮度影像 数据中的近红外波段像元数据。
[0214] 根据均值标准差法确定所述PBI指数的不透水层提取阔值T;具体地，对福亮度影 像数据进行快速统计获取所述PBI影像像元的均值y和标准差e;
[0215] 利用公式（7)确定所述PBI指数的不透水层提取阔值T。
[0216]T=y+e (7)
[0217] 该里，如图9所示，所述y值为0. 029,所述e值为0. 097,当阔值T为均值y与 一倍标准差e之和时，不透水层的提取效果最好。
[021引当所述PBI指数大于所阔值T时，确定所述影像像元为不透水层区域，所述不透水 层包括；亮色不透水层、暗色不透水层和藍色不透水层。其中，无山体阴影的不透水层提取 结果如图10所示。
[0219] 当所述PBI指数不大于所阔值T时，确定所述影像像元为±壤、植被等透水层区 域。
[0220]W上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围，凡在 本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护 范围之内。
【主权项】
1. 一种不透水层的提取方法，其特征在于，所述方法包括： 获取卫星影像数据； 计算所述影像数据在λ波段的大气层辐亮度P λ，将所述影像数据预处理为辐亮度影 像数据； 根据所述辐亮度影像数据中的绿光波段像元数据和近红外波段像元数据计算调整归 一化水体指数MNDWI ; 当所述MNDWI不大于O时，在所述辐亮度影像数据中选取各类地物的纯净像元数据； 利用建筑线方程、土壤线方程及土壤阴影线方程确定垂直不透水层指数PBI的系数η ; 利用公式PBI = nBlue-NIR计算PBI指数； 根据均值标准差法确定所述PBI指数的不透水层的提取阈值T ;其中， 当所述PBI指数大于所述阈值T时，确定影像像元为不透水层区域；所述λ取值为自 然数。2. 如权利要求1所述的方法，其特征在于，所述计算所述影像数据在λ波段的大气层 辐亮度P λ包括： 根据公式P ' λ = M pQMl-Ap计算未校正的λ波段大气层辐亮度P ' λ; 根据公式Ai = 对所述未校正的λ波段大气层辐亮度p ',进行校正，获取 校正后的λ波段的大气层辐亮度ρ λ;其中，所述M ρ为增益系数；所述〇。31为λ波段的所 述影像像元亮度DN值；所述Ap为偏置系数；所述Θ SE为当地的太阳高度角。3. 如权利要求1所述的方法，其特征在于，所述根据所述辐亮度影像数据中的绿光波 段像元数据和近红外波段像元数据计算调整归一化水体指数MNDWI包括： - MTP 利用公式皿《?7 = -^计算所述调整归一化水体指数MNDWI ;其中，所述 Green + NIR Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中 近红外波段像元数据。4. 如权利要求1所述的方法，其特征在于，所述利用建筑线方程确定垂直不透水层指 数PBI的系数η包括： 对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y = abx+b ; 当确定所述不透水层所在区域存在山体阴影时，构建土壤和山体阴影样本数据，对所 述土壤和山体阴影样本数据进行最小二乘拟合，确定土壤阴影线方程y = assx+b ; 利用公式a = 计算所述垂直不透水层指数PBI的系数η;其中，η与a相等。5. 如权利要求4所述的方法，其特征在于，当确定所述不透水层所在区域不存在山体 阴影时，所述方法还包括： 构建土壤样本数据，对所述土壤样本数据进行最小二乘拟合，确定土壤线方程y = asx+b ； 利用公式f = $ &计算所述垂直不透水层指数PBI的系数η ;其中，η与a'相等。6. 如权利要求1所述的方法，其特征在于，所述根据均值标准差法确定所述PBI指数的 不透水层提取阈值T包括： 对辐亮度影像数据进行统计获取PBI影像像元的均值μ和标准差ε ; 利用公式T= μ+ε确定所述PBI指数的不透水层提取阈值Τ。7. -种不透水层的提取装置，其特征在于，所述装置包括： 获取单元，所述获取单元用于获取卫星影像数据； 预处理单元，所述预处理单元用于计算所述影像数据在λ波段的大气层光谱反射率 P λ，将所述影像数据预处理为辐亮度影像数据； 计算单元，所述计算单元用于根据所述辐亮度影像数据中的绿光波段像元数据和近红 外波段像元数据计算调整归一化水体指数MNDWI ; 选取单元，所述选取单元用于确定当所述调整归一化水体指数不大于O时，在所述辐 亮度影像数据中选取各类地物的纯净像元数据； 第一确定单元，所述第一确定单元用于利用建筑线方程、土壤线方程及土壤阴影线方 程确定垂直不透水层指数PBI的系数η ; 第二确定单元，所述第二确定单元用于利用公式PBI = nBlue-NIR计算PBI指数； 第三确定单元，所述第三确定单元用于根据均值标准差法确定所述PBI指数的不透水 层的提取阈值T ; 第四确定单元，所述第四确定单元用于确定当所述PBI指数大于所述阈值T时，确定影 像像元为不透水层区域；其中，所述λ取值为自然数。8. 如权利要求7所述的装置，其特征在于，所述预处理单元具体用于： 根据公式P ' λ = M pQMl-Ap计算未校正的λ波段大气层辐亮度ρ ' λ; 根据公5对所述未校正的λ波段大气层辐亮度ρ ',进行校正，获取 校正后的λ波段的大气层辐亮度ρ λ;其中，所述M ρ为增益系数；所述〇。31为λ波段的DN 值；所述Ap为偏置系数；所述Θ SE为当地的太阳高度角。9. 如权利要求7所述的装置，其特征在于，所述计算单元具体用于： 利用公式M十算所述调整归一化水体指数MNDWI ;其中，所述 Green为所述辐亮度影像数据中的绿光波段像元数据，所述NIR为所述辐亮度影像数据中 的近红外波段像元数据。10. 如权利要求7所述的装置，其特征在于，所述第一确定单元具体用于： 对所述不透水层的样本散点进行线性最小二乘拟合，确定所述建筑线方程y = abx+b ; 当确定所述不透水层所在区域存在山体阴影时，构建土壤和山体阴影样本数据，对所 述土壤和山体阴影样本数据进行最小二乘拟合，确定土壤阴影线方程y = assx+b ; 利用公式a = 计算所述垂直不透水层指数PBI的系数η;其中，η与a相等； 当确定所述不透水层所在区域不存在山体阴影时，构建土壤样本数据，对所述土壤样 本数据进行最小二乘拟合，确定土壤线方程y = asx+b ; 利用公式Y = a& ^as计算所述垂直不透水层指数PBI的系数η;其中，η与a'相等。
【专利摘要】本发明提供了一种不透水层的提取方法及装置，包括：获取卫星影像数据；计算影像数据在λ波段的大气层辐亮度ρλ，将影像数据预处理为辐亮度影像数据；根据辐亮度影像数据中的绿光波段数据和近红外波段数据计算调整归一化水体指数MNDWI；当MNDWI不大于0时，在辐亮度影像数据中提取纯净像元数据；利用建筑线方程、土壤线方程及土壤阴影线方程确定PBI的系数n；利用公式PBI＝nBlue-NIR计算PBI指数；根据均值标准差法确定PBI指数的不透水层的提取阈值T；其中，当PBI指数大于阈值T时，确定影像为不透水层区域；所述λ取值为自然数；如此，不透水层的提取方法构建的PBI指数除了能够区分不透水层与水体、植被地物，还能够有效区分各类不透水层与裸土及阴影，降低不透水层的提取误差。
【IPC分类】G06K9/00
【公开号】CN104881659
【申请号】CN201510324293
【发明人】田玉刚, 徐韵, 杨晓楠, 邓雪彬
【申请人】中国地质大学（武汉）
【公开日】2015年9月2日
【申请日】2015年6月12日