CN107507200A

CN107507200A - 一种基于连通检测与噪声抑制的sar影像高精度大范围水域提取方法

Info

Publication number: CN107507200A
Application number: CN201710766993.7A
Authority: CN
Inventors: 周纪; 张继荣; 曹入尹
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2017-08-31
Filing date: 2017-08-31
Publication date: 2017-12-22
Anticipated expiration: 2037-08-31
Also published as: CN107507200B

Abstract

本发明公开了一种基于连通检测与噪声抑制的SAR影像高精度大范围水域提取方法。首先分别对待提取对象的SAR影像进行图像预处理，再并将其分割为数据层图像和地理信息模块两部分；对待提取对象的SRTM地形数据进行灰度转换、几何校正与配准以及拉伸处理，得到得到高程灰度图；基于双阈值提取得到初始水域图，并对其进行进行斑点去燥、膨胀操作、斑块去燥、腐蚀操作处理后，再结合初始水域图完成边缘细节恢复；最后，将其与地理信息模块进行像元级整合，得到最终提取结果。本发明在单景大范围的合成孔径雷达数据中提取河流、湖泊、水库等水域，且提取出来的水域精度大幅提高并且保留水域的精细边缘形态信息。

Description

一种基于连通检测与噪声抑制的SAR影像高精度大范围水域提取方法

技术领域

本发明属于合成孔径雷达(SAR)图像处理，具体涉及一种基于连通检测与噪声抑制的SAR影像高精度大范围水域提取。

背景技术

合成孔径雷达(SAR)是主动式微波成像系统，相对于遥感体系中其他的被动遥感系统，如搭载在Terra和Aqua卫星上的中分辨率成像光谱仪(MODIS)、搭载在美国陆地卫星4-5号的专题制图仪(TM)等，具有显著的优势。合成孔径雷达不受云雨天气影响，不受昼夜时间影响，可以全天时全天候获取数据，它是地表探测以及地表变化监测的有力工具。同时，水域的探测与变化监测在生产规划及灾害评估等方面作用显著。合成孔径雷达作为遥感应用中的一个重要传感器，在空间信息获取方面发挥着重要作用。

然而，SAR影像的成像过程，易受成像机理与地形自身特性影响，产生大量噪声。因SAR成像系统是基于相干原理的，而这一理论基础存在着原理性缺陷，该缺陷表现为：在雷达回波信号中，相邻像素点的灰度值会由于相干性而产生一些随机的变化，并且这种随机变化是围绕着某一均值而进行的，这样就在图像中产生了斑点噪声。而斑点噪声的产生是由于SAR成像所基于的相干原理所造成的缺陷，因此是不可避免的。从产生机理上讲斑点噪声是由于雷达目标回波信号的衰落现象所引起的。信号的衰落过程是这样产生的：同时被照射的有多个散射体，当雷达目标和雷达站之间具有相对运动时，这多个散射体与雷达之间具有不同的路程长和不同的径向速度，这使得雷达接收机接收到的信号产生一定的随机起伏，从而使SAR对目标散射系数的测量产生很大的偏差。最终表现在图像上，就产生了不可避免的斑点噪声现象。目前在SAR图像的应用前通常会进行斑点噪声抑制，其方法包括中值滤波、均值滤波、Lee滤波、Gamma MAP滤波、基于小波分解的阈值滤波等。以上的各种斑点噪声抑制滤波，都不会完全抑制斑点噪声，通常会引入影像融合的方法进一步抑制斑点噪声，但这对于单景大范围的研究区来说，很难实现。

目前，对SAR影像的水域提取方法有很多种，例如基于纹理信息的提取、结合地形信息的提取、独立成分分析、阈值法分割等，其中阈值法应用最为广泛。阈值法的依据是水体自身的散射特性，即水体的后向散射系数会小于非水体等其它地物。阈值法通过设定合适的阈值，将水体与非水体分开。例如对包含多个波段信息的多波段扫描影像(如TM、MODIS等)，可以通过波段间的组合运算得到一些特殊的描述因子来描述水体，例如用描述因子其中green为可见光绿波段，swir为短波红外并用green和swir波段来区分噪声中的阴影和道路。

用阈值法进行水域提取时，阈值的确定是最为关键的一步，现有的技术假设水体的散射回波强度与该强度下像元出现的概率呈Gamma分布，假设背景的散射回波强度与该强度下像元出现的概率呈高斯分布，又由于水体的散射回波强度相对非水体整体偏低，这样，研究区的回波强度与像元概率直方图会呈现一前一后双峰分布，阈值法的最优阈值确定往往在两峰相交的谷底处。当被提取目标面积占研究区面积30％以上时，会呈现很好的双峰效果，阈值的确定简单方便，提取效果较好。目前关于合成孔径雷达数据水域提取的研究区范围选取大都很小，如湖泊局部、河流交汇处局部等，其水体面积占研究区总面积比例较高，符合双峰阈值分割的假设。但是，当研究区范围扩大时，水域面积占研究区总面积远小于30％，并且很可能包含多条干流、支流、多个湖泊及水库等，阈值法的双峰完全消失；加之地形复杂的特性和合成孔径雷达自身的缺点，水域的提取要么遗漏，要么产生大量误判，产生大量干扰噪声，提取出来的水体无法用于大范围的水域监测。

发明内容

本发明的发明目的在于：针对上述存在的问题，提供一种在单景大范围的合成孔径雷达数据中提取河流、湖泊、水库等水域，并且使用一种有别于现有技术的去噪方法去除噪声和山体阴影等误判，使得提取出来的水域精度大幅提高并且保证水域的边缘形态信息。

由雷达方程可以得出后向散射系数的表达式为：其中σ⁰为后向散射系数，p_r为接收功率，P_t为发射功率，R为雷达到目标的单程距离，G为天线增益，λ为波长，ΔA为分辨单元。

雷达图像的灰度是地面目标后向散射回波强度在图像上的反映，它与接收信号的振幅成正相关。水体的后向散射系数会远小于裸土、植被等非水体，根据雷达参数与极化方式对雷达图像进行辐射校正，计算得到后向散射系数图像。由于本发明针对单景大范围雷达图像的水域提取，水体面积占研究区面积非常小，又因为单景图像既同时包含平原、丘陵、山脉等，地形情况复杂多样，阈值确定直方图已远不具备双峰的特点，兴趣区早已淹没在雷达图像的多类噪声中，即便找到了最合适的分割阈值，分割结果依然会包含大量噪声及误判，产生原因下文说明。

因此，本发明的思想是，既然无法利用双峰分割，则可先扩大常用回波阈值(避免水体漏判)的选取，保证初步提取出的对象包含河流、支流、湖泊及水库，并且保证足够多的轮廓信息量；之后，再对因此产生的大量噪声及误判进行分析去除。大范围初步提取出的水域效果不好，平原地区表现在图像冗乱，斑点斑块噪声大量存在，包含大量误判水体，山地地区表现在基本无法区分山体阴影与水体，产生大面积误判。

为了有效的抑制本发明的大范围水域提取处理时的斑点噪声，本发明先对雷达图像进行初步水域提取，兴趣区(可能为水体的像元)编辑为高亮其余编辑为背景，即对初步水域提取结果进行二值化处理，此时斑点噪声混杂在提取结果中，由于待提取对象为河流、支流、湖泊及水库等，虽然其面积相对研究区很小，但是，这些兴趣区的面积远大于斑点噪声的面积。支流及湖泊的面积大概为几万至上百万个像元，而斑点噪声的像元面积大多不会超过一百个像元。因此，本发明计算四邻域连通区(两个像素点是上下或左右相连的，则称其为四连通的)的面积，为了计算效率，将图像转换为逻辑二值图，将像元小于斑点噪声阈值的对象全部剔除，其中噪声阈值通常设置为遥感影像总面积的两千万分之一左右，即取值范围为T₁±δ，其中T₁表示遥感影像总面积的两千万分之一，δ表示阈值偏差(经验值)。

斑块噪声也是导致提取结果冗乱的一个重要因素，其部分是斑点噪声的扩大，更多的是由于平原地区丘陵等大面积地形起伏和近山脉地区地形变化造成的叠掩。这些斑块噪声大量存在，形状各异，其像元面积大于斑点噪声，斑块噪声的像元面积相对于主干河流、湖泊和水库等大面积的水域是非常小的。但是有的支流河道窄、水深浅，再加上枯水期水位下降，其后向散射系数显著增大，在加上桥梁的阻断，其表现在图像上就是支流或者多级支流之间的间断、不连续。斑块噪声的像元面积相对于这些间断的水域面积没有显著的量级差。如果继续像斑点噪声那样直接剔除小于某面积的对象，那么大量支流会被遗漏。况且即使非旱季，细小支流也会被后向散射系数较大的桥梁截断，仅从影像上看是非连通的，其连通面积也是小于某些斑块噪声。如果支流没有因为枯水期，或者桥梁的架设这些原因变得间断，那么它本身连通，加上与干流连通，其连通区面积会非常的大，远远大于斑块噪声。本发明用一种创新性的方法消除斑块噪声，流程可以概括为：连接兴趣区→剔除斑块→初步恢复细节→精细恢复细节。

连接兴趣区：对兴趣区(前景)进行形态学图像处理中的膨胀操作。在膨胀操作时，优选圆盘状半径为R_discal的结构元对二值图进行膨胀连接，其中R_discal选择为遥感影像的长或宽像素中较大者的1/2000左右，例如14个像元。

连接兴趣区后，支流断续处得到有效的连接，这样，斑块噪声的像素面积远小于膨胀连接后的兴趣区。再次计算各四邻域连通区像素面积，对小于斑块噪声阈值T₂(例如14000像素面积)的对象进行剔除。膨胀连接操作的同时，河流会丢失边缘细节，河流上桥梁丢失，对于分岔丰富的湖泊水库，分岔间完全黏连，丢失轮廓信息。

为初步恢复影像细节引入形态学图像处理中的腐蚀操作。为尽可能恢复河流的轮廓信息，或为下一步精细复原做准备，选用与膨胀连接相同参数的结构元进行腐蚀操作。

在膨胀与腐蚀过程之间进行了斑块噪声的剔除，但由于先膨胀再腐蚀的操作会填充图像内细小空洞、连接邻近物体、平滑边界，所带来的后果就是河流的桥梁消失，还有分岔密集的湖泊水库分岔部分会黏连在一起。为了恢复以上细节，进行下述运算：

其中，A_dilate_erode_(i,j)为经过初始分割、斑点噪声抑制、膨胀连接、斑块噪声抑制、腐蚀操作的矩阵；B_initial_seg_(i,j)为仅经过初始分割的矩阵。D_(i,j)为细节恢复后的输出矩阵，(i,j)表示像元坐标。即基于经过初始分割、斑点噪声抑制、膨胀连接、斑块噪声抑制、腐蚀操作后的二值图和初始分割后的二值图的各像素值得到恢复细节的二值图。

经过上述运算斑块噪声得到有效抑制，同时，保留了最精细的水域轮廓信息，桥梁细节完全得到恢复。

除去斑点与斑块噪声，山体阴影是影响SAR数据水域提取质量的关键噪声。因雷达采用侧视斜距投影成像方式，当雷达波束不能覆盖整个地表时，会出现雷达阴影。阴影出现在朝向远距的方向，在垂直地物背后或陡峭的坡地背后尤为明显，山体阴影处并无散射回波，后向散射系数近似为0，在影像上呈现为黑色。依据阈值提取水体的方法无法分辨山体阴影，山体阴影为块状，面积大，已经大于某些湖泊，远大于支流，成块出现非常聚集。目前去除SAR影像山体阴影多采用异源影像融合的方法，融合数字高程模型(DEM)数据去除山体阴影是一个常用的方法。本发明依然采用融合DEM去除山体阴影的方法，不同之处在于对DEM数据进行图像方面的处理，使其针对性更强，更适用于水域提取。DEM反映了研究区域地面高程的分布情况，从DEM数据中提取数字河网结构及相关流域信息是分布式水文模型开发与应用的基础。SRTM(ShuttleRadarTopographyMission)地形数据按精度可以分为SRTM1和SRTM3，分别对应的分辨率精度为30米和90米数据(目前公开数据为90米分辨率的数据)。本发明对SRTM3高程数据进行处理，为了从山区高程值较低的地区提取洼地和模拟河网，扩大山地地区高程之间的差异，先将DEM数据转换为相对高程下的灰度图，并对其进行黑白色阶的分离，进行对比度拉伸，具体拉伸处理方式如下：

其中，elevation_(i,j)为原始DEM灰度矩阵，elevationlog_(i,j)为对数拉伸后的灰度矩阵，NaN为遥感影像经过几何校正后边缘特有的空值。

经过拉伸处理后，将高亮处的坐标标记为非兴趣区，与水域初提取结果下的兴趣区进行融合判定，进行山体阴影的抑制。对于整个噪声去除流程，在水域最初阈值提取后即进行山体阴影的去除，大面积山区和山脉地区的山体阴影或得到直接剔除，或被小流域分割降级为斑块噪声，然后再进行前述斑点噪声和斑块噪声的抑制，完成整个去噪流程。

因此，本发明基于连通检测与噪声抑制的SAR影像高精度大范围水域提取方法包括下列步骤：

步骤1：输入待提取对象的SAR影像并进行图像预处理，得到SAR影像灰度图，其中图像预处理包括辐射校正、初步斑点噪声抑制、几何校正、对数刻度转换等；

并将预处理后的SAR影像分割为数据层图像和地理信息模块两部分，并记数据层图像为图像I_S；

步骤2：输入待提取对象的SRTM地形数据并转换为灰度图，然后行几何校正与配准处理后，再进行拉伸处理，得到高程灰度图I_D；

其中拉伸处理为：若当前像元的灰度值p_e(i,j)为空值，则将当前像元的灰度值置为0；否则将当前像元的灰度值设置为10*log|p_e(i,j)|；其中(i,j)表示像元坐标；

步骤3：结合图像I_S与高程灰度图I_D进行基于阈值的水域提取处理，得到初始水域图I_{W_0}，其中水域提取处理为：对于同一地理位置的像元坐标(i,j)，若I_S(i,j)≤T_S且I_D(i,j)≥T_D，则将当前像元标记为表示水体的灰度值；否则标记为非水体的灰度值；I_S(i,j)、I_D(i,j)分别表示图像I_S、I_D中位于像元坐标(i,j)处的灰度值，T_S、T_D为预设的图像阈值，设置T_S时，在尽量降低漏选率的前提下进行设置。

步骤4：对初始水域图I_{W_0}进行去燥及边缘细节恢复处理：

401：计算关于水体的四领域连通区域的像元面积，并剔除像元面积小于斑点噪声阈值的对象剔除；

402：采用圆盘状的结构元对步骤401处理后的图像的水体进行膨胀操作；

403：计算关于水体的四领域连通区域的像元面积，剔除像元面积小于斑块噪声阈值的对象剔除；

404：采用圆盘状的结构元对步骤403处理后的图像的水体进行腐蚀操作，得到水域图像I_{W_1}；

405：恢复边缘细节：对于同一地理位置的像元坐标(i,j)，若灰度值I_{W_0}(i,j)和I_{W_1}(i,j)均为水体，则将前像元标记为水体的灰度值；否则标记为非水体的灰度值，从而得到水域图I_W；

其中，步骤401、404的结构元的半径相同。

步骤5：将水域图I_W与地理信息模块进行像元级整合，得到最终提取结果。

进一步的，图像阈值T_S的优选设置方式为：对图像I_S的不同后向散射值进行像元数统计，记占比最多的后向散射值所对应的像元数为A，获取占比为0.03的像元数所对应的后向散射值，取其中的最小后向散射值作为图像阈值T_S。例如：计算I_S的直方图，直方图横坐标为后向散射系数(也就是回波强度，用灰度图的强度表征)，纵坐标为对应的像元数量，表示I_S中某种灰度强度下的像元个数，则T_S为纵坐标最大值的1/30所对应的横坐标(灰度值)。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

现有的合成孔径雷达水域提取办法，大都是围绕局部影像数据展开的相关研究，如湖泊局部，河流交汇处局部，水库等小范围的水域提取，关于大范围(200km×150km及以上)、全局影像的水域高精度提取研究还比较匮乏，本发明在单景大范围的合成孔径雷达数据中提取河流、湖泊、水库等水域，并且使用一种有别于现有技术的去噪方法去除噪声和山体阴影等误判，使得提取出来的水域精度大幅提高并且保留水域的精细边缘形态信息。

附图说明

图1是本发明的处理流程图；

图2嘉陵江、培江、渠江交汇处水域提取结果示意图，其中2-a为SAR原始图，2-b为水域提取结果图；

图3川西高原近横断山脉处提取结果示意图，其中3-a为SAR原始图，3-b为水域提取结果图；

图4噪声抑制效果示意图。

图5为局部光学影像(三岔湖)与水域提取细节示意图，其中5-a为水域提取结果图示意图(图中白色框所示)，5-b是提取结果的高分辨率光学影像图，5-c是提取结果的边缘细节图；

图6为局部光学影像(鲁班水库)与水域提取细节示意图，其中6-a为水域提取结果图(图中白色框所示)，6-b是提取结果的高分辨率光学影像图，6-c是提取结果的边缘细节图；

图7为局部光学影像(桥梁及轮廓信息)与水域提取细节示意图，其中7-a为水域提取结果图(图中白色框所示)，7-b为所提取的桥梁的高分辨率光学影像图，7-c为所提取的桥梁的边缘细节图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合实施方式和附图，对本发明作进一步地详细描述。

本实施例中选取了两个典型区域，分别为中国西南丘陵合川地区三江交汇处(104°3′40″E——107°1′47″E，29°0′3″N——31°9′10″N)和西高原近横断山脉处(102°0′44″E——104°2′38″E，28°8′40″N——30°8′26″N)。其特征是前者地势平坦、水域遍布，河网纵横；后者平原与山脉地区接壤，平原与山地共处一景，地形起伏变化大。示例主要用前者进行流程展示，后者仅给出提取结果展示。本发明的处理流程如图1所示。

本实施例中采用的数据包括：哨兵1号(Sentinel-1)卫星的SAR数据(Sentinel-1ALevel-1级产品)和高精度地形网格数据(SRTM3Sec 90m)。

(1)数据预处理：

在遥感图像处理平台(ENVI)中，对Sentinel-1A level-1级数据进行预处理：

步骤一：对原始影像(即SAR原始数据)进行辐射校正，消除辐射偏差得到影像的地物雷达后向散射系数Sigma0，得到辐射校正后的SAR影像(即Sigma0波段图像)。

步骤二：对辐射校正后的SAR影像进行初步斑点噪声抑制，例如选用参数为Windowsize：7×7、Target window size：3×3、sigma：0.9的Lee滤波器进行斑点噪声滤波。

步骤三：对初步斑点噪声抑制处理后的SAR影像进行几何校正。

步骤四：将几何校正后的SAR影像从线性刻度转换到对数刻度，并保存；本实施例中，将转换到对数刻度的SAR影像以以ENVI标准格式导出，为下一步ENVI中的数据格式转换做准备。

步骤五：添加SRTM 3Sec数字高程模型波段。DEM Model：SRTM 3Sec，MapProjection：WGS84。对其进行几何校正与配准，并以ENVI标准格式导出，为在ENVI中的DEM拉伸做准备。

(2)ENVI中的数据准备：

步骤一：将Sigma0波段导出为GeoTIFF格式，此格式将数据层矩阵和地理信息模块两个部分整合为一个整体。

步骤二：将数字高程模型中的Elevation波段(高程波段)在遥感处理软件中进行拉伸处理，。为减少数据量，将拉伸后的Elevation波段导出为JPG格式，即仅有数据层不包含地理信息模块。

(3)水域提取与去噪

步骤一：利用MATLAB将GeoTIFF格式数据剥离为数据层矩阵和地理信息模块两个部分，仅使用数据层矩阵进行处理，最后将数据层矩阵和地理信息模块重新整合即可还原地理信息。

步骤二：计算Sigma0矩阵与Elevation矩阵直方图，读取Sigma0阈值参数ValueA与DEM阈值参数ValueB。

步骤三：将满足Sigma0矩阵小于ValueA且Elevation矩阵大于ValueB的坐标记录为前景，其余记录为背景完成初步提取。并且将矩阵转换为逻辑二值矩阵。

步骤四：计算四邻域连通区域像元面积，将小于30像元的对象剔除。

步骤五：使用圆盘状，半径为14的结构元对前景进行(形态学图像处理中的)膨胀操作。

步骤六：再次计算四邻域连通区域像元面积，将小于14000像元的对象进行剔除。

步骤七：使用圆盘状，半径为14的结构元对留下的兴趣区进行(形态学图像处理中的)腐蚀操作。

步骤八：将步骤七计算后的矩阵与步骤三计算所得矩阵进行运算恢复边缘细节以及桥梁信息。

其中，A_dilate_erode_(i,j)为步骤七矩阵。B_initial_seg_(i,j)为步骤三矩阵。D_(i,j)为细节恢复后的输出矩阵。

步骤九：将D_(i,j)矩阵与地理信息模块(即每个象元的坐标信息)重新整合写为GeoTIFF文件，完成水域提取，得到水域提取结果。

本发明得到的水域提取结果具有大范围、高精度、低噪声的特点。以下由图2至图7展示提取结果，其中图2为嘉陵江、培江、渠江交汇处水域提取结果展示。图3为川西高原近横断山脉处提取结果展示。图4为噪声抑制效果展示。图5为提取细节展示(三岔湖)。图6为提取细节展示(鲁班水库)。图7为提取细节展示(桥梁及轮廓信息)。

以上所述，仅为本发明的具体实施方式，本说明书中所公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换；所公开的所有特征、或所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以任何方式组合。

Claims

1.一种基于连通检测与噪声抑制的SAR影像高精度大范围水域提取方法，其特征在于，包括下列步骤：

步骤3：结合图像I_S与高程灰度图I_D进行基于阈值的水域提取处理，得到初始水域图I_{W_0}，其中水域提取处理为：对于同一地理位置的像元坐标(i,j)，若I_S(i,j)≤T_S且I_D(i,j)≥T_D，则将当前像元标记为表示水体的灰度值；否则标记为非水体的灰度值；I_S(i,j)、I_D(i,j)分别表示图像I_S、I_D中位于像元坐标(i,j)处的灰度值，T_S、T_D为预设的图像阈值；

步骤4：对初始水域图I_{W_0}进行去燥及边缘细节恢复处理：

404：采用圆盘状的结构元对步骤403处理后的图像的水体进行腐蚀操作，得到水域图像I_{W_1}，其中所述结构元的半径与步骤401中的结构元相同；

2.如权利要求1所述的方法，其特征在于，所述图像阈值T_S的优选设置方式为：对图像I_S的不同后向散射值进行像元数统计，记占比最多的后向散射值所对应的像元数为A，获取占比为0.03A的像元数所对应的后向散射值，取其中的最小后向散射值作为图像阈值T_S。

3.如权利要求1所述的方法，其特征在于，斑点噪声阈值的取值范围为：[T₁-δ,T₁+δ]，其中T₁表示SAR影像的总像元面积的两千万分之一，δ表示阈值偏差。

4.如权利要求1所述的方法，其特征在于，斑块噪声阈值的取值范围为：[T₂-δ,T₂+δ]，其中T₂表示SAR影像的总像元面积的四百万分之一，δ表示阈值偏差。

5.如权利要求1所述的方法，其特征在于，步骤401和404的结构元的半径的取值范围为：[T₃-δ,T₃+δ]，其中T₃表示SAR影像的长和宽的像元中的较大者的两千分之一，δ表示阈值偏差。