CN112166693B

CN112166693B - 基于小卫星的区域地表水资源遥感监测方法

Info

Publication number: CN112166693B
Application number: CN201218002156.5A
Authority: CN
Inventors: 吴双; 周会珍; 王智勇; 冉琼; 于冰洋; 文强
Original assignee: Twenty First Century Aerospace Technology Co ltd
Current assignee: Twenty First Century Aerospace Technology Co ltd
Priority date: 2012-06-29
Filing date: 2012-06-29
Publication date: 2014-10-22
Anticipated expiration: 2032-06-29

Abstract

本发明属于遥感应用技术领域，公开了一种基于小卫星的区域地表水资源遥感监测方法，包括以下步骤：1)确定监测时相和监测区域，采集并甄选小卫星影像和辅助资料；2)小卫星多光谱影像正射校正；3)小卫星多光谱影像辐射归一化；4)地表水资源信息提取。本发明针对小卫星中高分辨率、宽覆盖、高时间分辨率等成像特点，形成了一套适合小卫星的高效的区域地表水资源遥感监测方法，弥补了当前地表水资源遥感监测方法在大区域应用的不足，更准确地反映了区域地表水资源的空间分布格局及其演变特点，实现了连续多年多时相区域中高分辨率遥感影像的地表水资源遥感监测，保证了结果的准确性、连续性和可比性。

Description

基于小卫星的区域地表水资源遥感监测方法

技术领域

本发明属于遥感应用技术领域，涉及一种基于中高分辨率遥感数据高动态监测区域地表水资源的方法。

背景技术

水是维持自然界一切生命和社会经济持续发展所必需的资源。水资源短缺与生态环境恶化是当今世界范围内普遍面临的严重问题。特别是在我国，水资源短缺已严重影响了国家经济发展和社会和谐。

地表水作为水资源的重要组成部分能从特定角度反映水资源状况，遥感技术具有宏观、科学、快速、大区域的特点，成为地表水资源宏观监测的必选技术。在利用遥感技术监测地表水资源时，人们对水资源的监测需求与遥感数据之间存在矛盾，主要表现为以下两点：

1)当前使用的遥感数据主要来源于国外卫星数据，其较高的地面分辨率和较高的时间分辨率难以兼备，从而导致监测结果要么更新慢，要么精度低。

2)当前对于地表水资源的研究多停留在较小尺度的行政区域内，且往往针对某一个或几个离散时间点进行研究，未形成连续性监测，与水文要素的监测原则不一致，监测结果参考价值不够，可比性差。

因此，迫切需要新的遥感数据源和新的水资源遥感监测方法以更好地开展区域的水资源遥感监测研究。国产小卫星具有中高分辨率、宽覆盖、高时间分辨率等特点，为区域水资源遥感监测奠定了数据基础。

发明内容

本发明的目的是提供一种基于小卫星的区域地表水资源遥感监测方法，该方法更准确地反映了区域地表水资源的空间分布格局及其演变特点。

所述基于小卫星的区域地表水资源遥感监测方法包括以下步骤：

步骤一：确定监测时相和监测区域，采集并甄选小卫星影像和辅助资料；

步骤二：对小卫星多光谱影像进行正射校正

对于由单成像仪或两个成水平方向的成像仪获取的小卫星多光谱影像，采用通用的正射校正方法进行正射校正；

对于两个或两个以上呈一定夹角(0，180°)的成像仪获取的小卫星多光谱影像，采用以下方法进行正射校正：

(1)选择每个成像仪分别获取的标准景影像；

(2)分别对每个标准景影像进行列分割；

(3)对分割后的影像块分别进行正射校正；

(4)通过重采样将正射影像输出；

(5)对分割景正射影像进行镶嵌，得到宽刈幅标准景正射校正影像；

步骤三：对小卫星多光谱影像进行辐射归一化处理；

步骤四：提取小卫星影像区域地表水资源信息。

与现有方法相比，本发明的优势在于：本发明针对区域地表水文要素特点，基于小卫星中高分辨率、宽覆盖、高时间分辨率等成像特点，形成了一套适合小卫星的高效的区域地表水资源遥感监测方法，弥补了当前地表水资源遥感监测方法在大区域应用的不足，更准确地反映了区域地表水资源的空间分布格局及其演变特点，实现了连续多年多时相区域中高分辨率遥感影像的地表水资源遥感监测，保证了结果的准确性、连续性和可比性。

优选的，对分割后的影像块采用商业遥感图像处理软件PCI的有理函数模型，通过参考控制点GCP和数字高程模型完成大幅宽遥感影像正射校正，参考影像为1∶10万地形图，参考DEM为ASTER 30m分辨率的数字高程模型。采用该有理函数模型进行正射校正能够得到较高精度的校正结果。

优选的，所述小卫星多光谱影像辐射归一化包括如下步骤：

(1)对经步骤二进行正射校正后的影像进行图像配准；

(2)确定参考影像与待校正影像；

(3)获取参考影像与待校正影像的伪不变特征点集，方法如下：

①初步选择参考影像与待校正影像的伪不变特征点，公式如下：

f≤ratio and(DN_ref＜thres_1 or DN_ref＞thres_2)

thres_1＝DN_min+(DN_max-DN_min)/4

thres_2＝DN_max-(DN_max-DN_min)/4

其中，f为参考影像某像元值与所对应的待校正影像DN值的比值，ratio为参考影像平均DN值与待校正影像平均DN值的比值，DN_ref为参考影像的DN值，DN_min和DN_min分别为参考影像DN值的最小值和最大值；

对参考影像的每个像元点按上述公式进行计算，符合公式要求的像元点作为初步选择的伪不变特征点，与参考影像伪不变特征点相对应位置的待校正影像的点作为初步选择的待校正影像的伪不变特征点；

②伪不变特征点进一步迭代选择

对初步选择的参考影像的伪不变特征点和相应的待校正影像的伪不变特征点做主成分分析得到第一主成分pc1，计算伪不变特征点到pc1的距离l，若l≤4，则保留该点，否则剔除；重复上述过程，直到最后所有的点满足该公式要求，则得到最终的参考影像和待校正影像的伪不变特征点集；

(4)确立参考影像与待归一化影像的线性关系

利用参考影像和待校正影像最终获取的伪不变特征点，采用最小二乘法建立线性关系式y＝a×x+b，其中，参考影像的伪不变特征点为因变量y，待校正影像伪不变特征点为自变量x，得到关系式的截距b和斜率a；

(5)以整幅待归一化影像为自变量，利用已得线性关系式获取归一化处理后影像。

上述辐射归一化方法，一方面完善了区域地表水信息提取流程，另一方面提升了影像制图效果，提高了地表水信息提取精度。

优选的，所述参考影像和待校正影像的确定步骤如下：首先，计算多幅影像的平均梯度G，选择平均梯度最大的三幅影像；然后，对已选取的三幅影像计算影像灰度均值mean，灰度均值最大的影像为参考影像，其余的作为待校正影像。

优选的，所述步骤四中，对于河流，首先采用归一化植被指数阈值法得到河流水系面层，然后将河流面层转换成线层矢量，继而采用设置河流结点计算河流长度均值的提取方法，即：一条河流如果有其他支流的汇入，则设置一个结点，相邻的两个结点间再取中间位置点，以这三个点水面宽度测量值的平均值作为该河段水面宽度，再乘以该河段长度得到河流水面面积。

上述地表水信息提取方法中，河流水系的提取克服了传统方法提取结果不连续以及受河流两侧阴影影响的不足。

优选的，所述归一化植被指数阈值设定为：NDVI≤0.1。

优选的，所述步骤四中，首先采集并甄选监测时相和监测区域内的全色影像，然后利用步骤三处理得到的多光谱影像与全色影像进行融合，采用面向对象分类的方法提取地表水信息，步骤如下：

对枯水期和丰水期获得的遥感影像，按照监测区域的地貌特征分为枯水期山区、枯水期平原、丰水期山区、丰水期平原四种情况，采用不同的分割尺度进行分割，然后，依据三种对象分类特征——近红外图层均值、近红外图层方差、各波段图层纹理特征选取各水体类型的训练样本，采用K近邻监督分类方法进行分类，从而得到各种水体提取结果。

上述地表水信息提取方法可满足区域中的重点地区，或者对较小的水体(如坑塘)等有较高要求的地区的水信息提取。

优选的，所述分割尺度如下：丰水期平原地区：35，丰水期山区：100，枯水期山区：30，枯水期平原：20。

附图说明

图1为本发明小卫星区域地表水资源遥感监测方法流程图；

图2为本发明小卫星多光谱影像辐射归一化流程图；

图3为本发明小卫星区域地表水信息提取流程图；

图4为基于北京一号小卫星的2007年北京及周边地区五大流域地表水资源分布专题图。

具体实施方式

下面以北京一号小卫星对北京及周边五大流域地表水资源进行遥感监测为例，结合附图对本发明作进一步详细的说明。

图1为基于小卫星的区域地表水资源遥感监测方法流程图，包括以下步骤：

步骤一：确定监测时相和监测区域，采集并甄选小卫星影像和辅助资料

地表水资源分布受降水、径流以及蒸发等水文要素的影响而呈现或枯或丰的特点，枯水期内，流域地表水最少，主要依靠地下水补给水源；丰水期内，流域地表水量丰富，江河水流主要依靠降雨或融雪补给，延续时间较长。因此，将地表水资源的遥感监测时相确定为枯水期和丰水期，即每年开展丰水期监测和枯水期监测。

(1)确定监测时相和监测区域

由于流域是水资源管理的最佳自然单元，河流在流域内保持着它的自然完整性。因此，将流域作为地表水资源遥感监测区域。比如：本实施例确定北京及周边五大流域作为地表水资源的遥感监测区域。

北京及周边地区降水多集中在夏季和秋初，为丰水期；冬季和春节降水较少，是枯水期，但冬季由于气温较低，常有结冰现象。基于此，确定影像的采集时限枯水期为3-5月份，丰水期为6-9月份。

(2)采集并甄选采集时限内影像

采集监测时相和监测区域内小卫星的多光谱影像并进行甄选。甄选采集时限内影像的原则如下：a)影像清晰，云量少于20％。b)监测区域内影像时相即影像拍摄时间跨度小于2个月。c)不同年份同一区域内的影像拍摄时间月份跨度小于2个月。d)影像的几何校正精度一致。

(3)获取辅助资料

获取外业检查点、中分辨率数字高程模型和行政地图。

步骤二：小卫星多光谱影像正射校正

北京及周边地区五大流域山区和平原兼有，地貌变化多样，需要对小卫星多光谱影像进行正射校正。

小卫星搭载的遥感器可分为以下几种形式：单成像仪，两个成水平(180°)方向的成像仪，两个或两个以上呈一定夹角(0，180°)的成像仪。对于由单成像仪或两个成水平(180°)方向的成像仪获取的小卫星多光谱影像，采用通用的正射校正方法——直接利用商业遥感图像处理软件ERDAS，采用通用推扫式模型对影像进行正射校正，参见张兆霞，朱龙文的“基于遥感卫星的系列正射影像图制作”(《遥感与航空摄影》，2004年第3期)一文。

北京一号小卫星多光谱遥感器为宽刈幅双成像仪器的CCD线阵，相对于地平面是交叉倾斜投影，使得地形高差对左右影像造成的投影差畸变是反向的，通用的正射校正方法难以满足校正精度，需要修正成像仪自身的侧视角度。因此，采用如下方式对多光谱影像进行正射校正：

(1)选择左右成像仪分别获取的标准景影像

对每个成像仪获取的影像单独进行校正处理，选择每个成像仪获取的影像数据，数据大小为标准景(9984*9984像元)。

(2)对标准景影像进行列分割

为了保证正射校正模型对地形起伏、地球曲率等畸变的拟合精度，需要对单成像仪的影像在列方向进行分割，使其分成左右两个影像块，每个影像块大小为9984*5500，中间重叠1000个像元。

(3)对分割后的影像块分别进行正射校正

对分割后的影像块分别采用商业遥感图像处理软件PCI的有理函数模型(Retional Function Model，RFM)或商业遥感图像处理软件ERDAS的通用推扫式模型，通过参考控制点GCP和数字高程模型(DEM)完成大幅宽遥感影像正射校正。采用有理函数模型时，参考影像为1∶10万地形图，参考DEM为ASTER 30m分辨率的数字高程模型。

(4)正射影像采样输出

通过重采样将正射影像输出。

(5)分割景正射影像镶嵌

由于两个分割景之间保留了1000列的重叠部分，只需选取10个左右的同名控制点，就能达到镶嵌精度的要求，因此，采用PCI的镶嵌功能，通过选取少量的同名控制点，将经过正射校正的2个分割景影像进行镶嵌处理，即得到了宽刈幅标准景正射校正影像。

上述小卫星多光谱影像按成像仪单独校正、分割标准景单独校正的正射校正方法进行正射校正，提高了影像几何精度。

步骤三：小卫星多光谱影像辐射归一化

一方面，北京及周边五大流域地区范围广阔，即使宽刈幅的北京一号多光谱影像也无法用一幅覆盖整个区域。另一方面，在地物的连续动态监测中，不同时相拍摄的影像会受传感器本身、太阳高度角、大气、地形等因素对地面光谱反射信号的影响，导致多时相影像同一目标的光谱特征存在差异，一方面影响了影像制图，另一方面影响了地表水信息的提取。因此，为了提高影像制图效果和地表水提取精度，本发明提出了小卫星多光谱影像的辐射归一化方法，图2为北京一号小卫星多光谱影像辐射归一化流程。步骤如下：

(1)图像配准

对选取的监测区域内多幅多时相影像进行图像配准，配准精度控制在1个像元以内。

(2)参考影像与待校正影像的确定

从多幅多时相影像中选取图像清晰、整体亮度高的影像作为参考影像，采用平均梯度G以及影像灰度均值mean做为评判指标，其中，平均梯度G反映了图像的清晰度水平，平均梯度越大，图像的锐化效果越强，局部对比度越高；影像灰度均值反映了图像的整体亮度。平均梯度G的计算公式如下：

式中，G表示图像的平均梯度值，Z(x_i，y_j)为图像在第i行，第j列的灰度值，M为图像的总行数，N为图像总列数，i，j分别表示依次取的像素行列号，δ表示对函数求偏微分。

参考影像和待校正影像的确定步骤如下：首先，计算多幅影像的平均梯度G，选择平均梯度最大的三幅影像；然后，对已选取的三幅影像计算影像灰度均值mean，灰度均值最大的影像为参考影像，其余的作为待校正影像。

(3)伪不变特征点自动获取

本发明建立了主成分分析耦合迭代控制的伪不变特征点自动获取模型，利用该模型得到参考影像的伪不变特征点集以及待校正影像相对应的伪不变特征点集，进行辐射归一化。模型如下：

①初步选择参考影像与待校正影像的伪不变特征点

公式如下：

f≤ratio and(DN_ref＜thres_1 or DN_ref＞thres_2)

thres_1＝DN_min+(DN_max-DN_min)/4

thres_2＝DN_max-(DN_max-DN_min)/4

其中，f为参考影像某像元值与所对应的待校正影像DN值的比值，ratio为参考影像平均DN值与待校正影像平均DN值的比值，DN_ref为参考影像的DN值，DN_min和DN_min分别为参考影像DN值的最小值和最大值。

②伪不变特征点进一步迭代选择

(4)确立参考影像与待归一化影像的线性关系

利用参考影像和待校正影像最终获取的伪不变特征点，采用最小二乘法建立线性关系式y＝a×x+b，其中，参考影像的伪不变特征点为因变量y，待校正影像的伪不变特征点为自变量x，得到关系式的截距b和斜率a。

(5)获取辐射归一化处理后影像

以整幅待归一化影像为自变量，利用已得线性关系式获取归一化处理后影像。

步骤四：小卫星影像区域地表水资源信息提取

对于监测区域内面积较大的水体——水库、湖泊、河流和成片的坑塘，可以直接利用小卫星中分辨率的多光谱影像进行信息提取；而对于区域中的重点地区，或者对较小的水体(如坑塘)等有较高要求的地区，则需要借助4米全色影像，通过全色与多光谱影像融合进行地表水信息提取。

本实施例中，北京及周边五大流域地区以河流、水库、湖泊、坑塘四类地表水体进行分类，采用的地表水信息提取方法如下：

北京境外的五大流域地区，直接利用步骤三处理得到的多光谱影像进行信息提取，具体流程见图3(a)。其中，对于水库、湖泊和成片的坑塘，直接采用归一化植被指数(NDVI)阈值法，归一化植被指数阈值设定为：NDVI≤0.1；对于河流，首先采用归一化植被指数阈值法得到河流水系面层，为了减少传统方法提取结果不连续以及河流两侧阴影的影响，将河流面层转换成线层矢量，继而采用设置河流结点计算河流长度均值的提取方法，即：一条河流如果有其他支流的汇入，则设置一个结点，相邻的两个结点间再取中间位置点，以这三个点水面宽度测量值的平均值作为该河段水面宽度，再乘以该河段长度得到河流水面面积。

北京境内地区为监测区域中的重点地区，利用步骤三处理得到的多光谱影像与采集并处理后的4米全色影像进行融合，采用面向对象分类的方法提取地表水信息，具体流程见图3(b)。对枯水期和丰水期获得的遥感影像，按照北京的地貌特征(具体做法是：利用等高线图，以100米等高线作为山区和平原的划分线)分为枯水期山区、枯水期平原、丰水期山区、丰水期平原四种情况，采用不同的分割尺度进行分割，优选的分割尺度如下：丰水期平原地区：35，丰水期山区：100，枯水期山区：30，枯水期平原：20。然后，依据三种对象分类特征——近红外图层均值、近红外图层方差、各波段图层纹理特征选取各水体类型的训练样本，采用K近邻监督分类方法进行分类，从而得到各种水体提取结果。

利用已得各种水体的提取结果(即各种水体类型的面层矢量)，就可计算各种水体面积，制作水体分布专题图。

图4为基于北京一号小卫星的2007年北京及周边地区五大流域地表水资源分布专题图，其中图4(a)为枯水期，(b)为丰水期。从图上可以清楚地看出河流、水库、湖泊流域分布特征——五大流域地表水资源空间分布不均匀，西部、北部的山区主要以河流和水库分布为主，湖泊、坑塘几乎没有分布；东南部平原地区以湖泊、坑塘分布为主。北运河流域河网密集，坑塘众多。潮白河流域下游地区河流几乎断流，主要以坑塘为主。永定河流域河流比较稀疏，甚至出现“无流区”，坑塘分布很少。大清河流域东南边缘、蓟运河流域的西部和南部也有坑塘分布。

说明书中描述的只是该发明的具体实施方式。虽然结合附图描述了本发明的实施方式，但是本领域内熟练的技术人员可以在所附权利要求的范围内作出各种变形或修改。

Claims

1.一种基于小卫星的区域地表水资源监测方法，其特征在于包括以下步骤：

步骤一：确定监测时相和监测区域，采集并甄选小卫星多光谱影像和辅助资料；

步骤二：对小卫星多光谱影像进行正射校正；

对于两个或两个以上呈一定夹角的成像仪获取的小卫星多光谱影像，所述夹角在0～180°之间，采用以下方法进行正射校正：

(1)选择每个成像仪分别获取的标准景影像；

(2)分别对每个标准景影像进行列分割；

(3)对分割后的影像块分别进行正射校正；

(4)通过重采样将分割景正射影像输出；

步骤三：对小卫星多光谱影像进行辐射归一化处理，包括如下步骤：

(1)对经步骤二进行正射校正后的影像进行图像配准；

(2)确定参考影像与待校正影像；

f≤ratio and(DN_ref＜thres_1or DN_ref＞thres_2)

thres_1＝DN_min+(DN_max-DN_min)/4

thres_2＝DN_max-(DN_max-DN_min)/4

其中，f为参考影像某像元点的DN值与所对应的待校正影像像元点的DN值的比值，ratio为参考影像平均DN值与待校正影像平均DN值的比值，DN_ref为参考影像的DN值，DN_min和DN_max分别为参考影像DN值的最小值和最大值；

②伪不变特征点进一步迭代选择

(4)确立参考影像与待归一化影像的线性关系

(5)以整幅待归一化影像为自变量，利用已得线性关系式获取归一化处理后影像；步骤四：提取小卫星多光谱影像区域地表水资源信息。

2.如权利要求1所述的基于小卫星的区域地表水资源监测方法，其特征在于对分割后的影像块采用遥感图像处理软件PCI的有理函数模型，通过参考控制点GCP和数字高程模型完成宽刈幅遥感影像正射校正，参考影像为1∶10万地形图，参考数字高程模型为ASTER 30m分辨率的数字高程模型。

3.如权利要求1所述的基于小卫星的区域地表水资源监测方法，其特征在于所述参考影像和待校正影像的确定步骤如下：首先，计算多幅影像的平均梯度G，选择平均梯度最大的三幅影像；然后，对已选取的三幅影像计算影像灰度均值mean，灰度均值最大的影像为参考影像，其余的作为待校正影像。

4.如权利要求1所述的基于小卫星的区域地表水资源监测方法，其特征在于所述步骤四中，对于河流，首先采用归一化植被指数阈值法得到河流水系面层，然后将河流面层转换成线层矢量，继而采用设置河流结点计算河流长度均值的提取方法，即：一条河流如果有其他支流的汇入，则设置一个结点，相邻的两个结点间的河流为一个河段，

河段取中间位置点，以这三个点水面宽度测量值的平均值作为该河段水面宽度，再乘以该河段长度得到河流水面面积。

5.如权利要求4所述的基于小卫星的区域地表水资源监测方法，其特征在于所述归一化植被指数阈值设定为：NDVI≤0.1。

6.如权利要求1所述的基于小卫星的区域地表水资源监测方法，其特征在于所述步骤四中，首先采集并甄选监测时相和监测区域内的全色影像，然后利用步骤三处理得到的多光谱影像与全色影像进行融合，采用面向对象分类的方法提取地表水信息，步骤如下：

7.如权利要求6所述的基于小卫星的区域地表水资源监测方法，其特征在于所述面向对象分类的分割尺度如下：丰水期平原地区：35，丰水期山区：100，枯水期山区：30，枯水期平原：20。