核估计的最佳阈值SAR图像海岸线提取

(辽宁工程技术大学测绘与地理科学学院, 辽宁阜新 123000)

0 引言

海岸线的实时监测对国民经济发展来说是非常重要的，而当今遥感技术是完成此项任务的最佳手段[1]。在众多传感器中，合成孔径雷达(Synthetic Aperture Radar, SAR)相比于光学传感器可以在不受时间和天气的限制下连续工作，因此其更加适合作为海岸线监测的原始数据。但SAR图像固有的斑点噪声和海面的不均匀散射特性极大地增加了海岸线提取的难度，因此准确、迅速、自动的SAR图像海岸线提取仍然是一个具有挑战性的研究领域。

到目前为止，SAR图像海岸线提取方法主要分为3大类：基于边缘、基于区域和基于阈值的方法[2]。基于边缘的方法着重研究具有灰度跳变的边缘，故其具有高效率和高定位精度的优点[3-5]。但它们最致命的缺陷是它们检测出的海岸线是不连续的。更重要的是，有许多位于陆地区域的虚假海岸线也被检测了出来，这无疑会对后续的真实海岸线提取和判读造成不小的困难。

为了克服这些缺陷，基于区域的方法就应运而生了。它们通常以图像全局或某些同质区域为研究对象，故它们具有抗噪性好以及检测边缘连续的优点。其中它们又可再细分为：基于统计模型的方法、水平集法和其他法[6]。其中，其他法具有耗时长、难以适用于工程实践的缺点[7]，此外，其还需要预处理操作(如滤波)来进行辅助[8-9]，但这就使得边缘的细节信息不可避免地遭到了一定程度的破坏，进而使提取结果的精度下降。而由于基于统计模型的方法中参数的数值求解较为繁琐且耗时较长，故需要进一步研究与优化才能被广泛应用[10]。而赵泉华等[11]所用G0分布的参数估计虽然简单，但他们仅仅根据参数的变化而确定水陆边界点使得误差很大，提取的精度较低。水平集法的参数变化对SAR图像中较强斑点噪声的敏感性高，因此仍有待改进[12-14]。

由于图像中只有海陆两类地物，因此将图像分类是海岸线提取的最直接手段，而阈值法则是将图像分为两类的最快方法，且其可以不需要使用滤波等降低精度的预处理，故可以保留精细的海岸线细节。Al Fugura等[15]将Lee滤波后的图像的所有像素的灰度均值作为阈值，但这种取值原则并不具有普适性。而对于著名二值化算法OTSU法，很多学者都应用它进行了海陆分割[16-18]。但其并非是针对海陆分割而设计的，故其获得的阈值并非最佳阈值。为了寻找最佳阈值，陈祥等[19]在OTSU分割的基础上获得了最终阈值，但因OTSU获得的阈值本身距最佳阈值的差距就较大，且其只考虑了海水的分布而没有考虑陆地，故分割效果仍然不理想。为了兼顾海水与陆地的分布特性，Liu和Jezek[20]假设整幅图像服从双峰高斯分布，并认为双峰中间的“谷”所对应的灰度值就是最佳阈值。但在很多情况下，由于陆地上复杂的地物分布与不同的散射机制，使得陆地区域是不服从高斯分布的，进而使得到的阈值也是错误的。为了顾及不同的散射机制与统计分布，Nunziata等[21-22]组成的研究小组根据极化SAR数据的物理模型确定最佳阈值，对阈值分割法作出了较大贡献[23-24]。但该方法仍然受不同统计分布的限制，无法获得在任意情况下的最佳阈值。

为了摆脱分布的限制，本文借助非参数核估计为基本工具[25]，在分析了海水和陆地的分布规律后利用跳变检测的基本思想估计出了海水和陆地分布的临界点[26]，即最佳阈值。可为阈值法提取海岸线提供一种新思路。

1 差值核估计

对于n个数据点组成的点集{(xi,yi),i=1, 2,…,n}来说，每一个xi和yi都具有一一对应的关系，且有不确定因素所造成的微小偏差。那么可以把它们的关系描述为如下的模型：

yi=f(xi)+εi,i=1,2,…,n

(1)

式中，f表示xi和yi对应的关系，i为数据点的索引，ε表示不确定因素对第i个点的值所造成的偏差。假设某一点的y值相对于其他点发生了明显的跳变，那么可以利用差值核估计来进行跳变检测，分别估计出跳变点的所在位置与相应的跳变幅度。

1.1 核函数

(2)

式中，k(>0)用以调整xi+k所获得的权重幅度，这里取k=0.618。采用式(2)可以将点i的左右两个邻域内的数据点的y值进行加权，使越接近xi的xi+k获得的核函数值(权重)越大，反之越小。随着xi逐渐接近跳变点的所在位置，会使得两个邻域内数据点的加权均值之差的绝对值增大的速度呈现越来越快的非线性增长。最终可使估计出的跳变幅度在靠近真实跳变位置时更大，远离时更小，从而估计出的跳变位置更加接近真实跳变位置。图1为邻域内不同数据点的所在位置获得的权重示意图，其中黑点代表xi，红点代表xi+k，而绿线代表权重大小，越长的权重越大，反之越小。

图1 邻域内不同数据点的权重

1.2 跳变检测

以点i为中心点，其差值核估计为

(3)

(4)

M(xi)=|M1(xi)-M2(xi)|

(5)

式(3)和式(4)分别表示点i左邻域内和右邻域内数据点y值的加权平均值；式(5)表示两个邻域内加权均值之差的绝对值。

假设数据点集中只有一处跳变，那么相对于给定尺寸的邻域，对远离跳变位置的数据点，其左右邻域内的数据点值一般都是平稳变化的，因此其M值较小且在误差范围内恒定不变；对靠近跳变位置的数据点，存在某一邻域内包含跳变前与跳变后的数据点，使得M值变大；而对于跳变点，会使得两个邻域内分别全部包含的是跳变前的数据点和跳变后的数据点，这时M取得最大值。由此可知最大M值所对的位置(x值)即为跳变位置i*，而跳变位置对应的最大M值即为跳变幅度s：

(6)

s=M(xi*)

(7)

以100个数据点为例，它们的横坐标x分别是它们的序数。而x与纵坐标y的关系满足：当x≤50时，y=x；而当x>50时，y=150。因此该例中的第50个点就是跳变点，而跳变的幅度为150-50=100，如图2(a)所示。对每点进行差值核估计后的结果如图2(b)所示，这里带宽取h=10，可以看出在跳变点附近时M值开始增大，而跳变点处M值则达到了最大(接近100)，因此较好地估计出了跳变位置与跳变幅度。

(a) 原数据点集

(b) 经跳变检测后的数据点集图2 差值核估计的跳变检测

2 海岸线提取

2.1 最佳阈值分割

在高分辨率雷达成像后的SAR强度图像中，海面杂波的分布呈现出一种有不同程度“拖尾”特性的近似高斯分布的现象[19]。由于强散射与斑点噪声，陆地上灰度值的统计特性则是均值和方差都较大。这就使得陆地中某些相对较暗的像素的灰度值与海面中某些相对较亮的像素的灰度值很接近或相同，表现在分布直方图上则为：隶属于陆地分布的线段与隶属于海面分布的线段相交在了陆地与海面中互相接近的像素灰度值附近。当陆地整体相对较亮时，该相交处通常位于近似高斯分布中的“尾”处；而当陆地整体相对较暗时，较暗像素的数目更多，因此通常在“尾”之前就相交了。以图3(a)中256×256像素的SAR图像为例，该现象可见于图3(b)中的直方图，且红色箭头所指处即为相交处。

(a) SAR图像

(b) SAR图像的直方图图3 SAR图像与其直方图

而直方图的分布特性为：绝大多数位于海面的像素的灰度值都小于该相交处所对应的灰度值(横坐标)，而绝大多数位于陆地的像素的灰度值都大于该相交处所对应的灰度值。因此可以认为该相交处是海面和陆地分布的临界点，那么该临界点所对应的灰度值就可以作为海陆分割的阈值。为了自适应地确定该阈值，需要分析该临界点周围的几何特性。通过观察直方图可知，在该临界点左侧一定范围的邻域内由线段组成的“坡”通常较陡，而右侧的“坡”通常较缓。即该点是“坡”从陡到缓的转折点，而纵观整个直方图，以一定范围的邻域来衡量每个点的几何特性时，没有其他点有该临界点的几何特性。所以可认为该直方图中只有一处坡度从大到小的跳变点。为了衡量坡度，可以以直方图中的某点为中心点，计算某侧邻域内所有点的纵坐标与中心点纵坐标之差的加权平均值。若该加权平均值的绝对值较大，则某侧邻域内坡度较陡，反之坡度较缓。设{(xi,Yi),i=1, 2, …, 256}为组成直方图的点集，x和Y分别为横坐标与纵坐标，那么利用差值核估计来检测从大到小的坡度跳变点：

(8)

(9)

M(xi)=M2(xi)-M1(xi)

(10)

其中，式(8)表示直方图中点i的纵坐标与点i左邻域内所有点的纵坐标之差的加权平均值；式(9)表示点i右邻域内所有点的纵坐标与点i的纵坐标之差的加权平均值；式(10)表示前两个加权平均值之差。

那么相对于给定尺寸的邻域，当远离坡度从大到小的跳变点以及远离近似高斯分布的峰时，其M1与M2相差不大，因此其M较小且波动不大；当靠近近似高斯分布的峰时，存在某一邻域内包含该峰两侧的“坡”的点，即某一邻域内会同时包含纵坐标比中心点大的与纵坐标比中心点小的点，使得M减小；而对于近似高斯分布的峰，会使得两个邻域内分别包含的是峰左侧的“上坡”的点与峰右侧的“下坡”的点，即右邻域内所有点的平均纵坐标都小于中心点很多，而中心点的纵坐标都远大于左邻域内所有的点平均纵坐标，故这时M取得最小值；当靠近坡度从大到小的跳变点时，会使得左、右邻域内分别包含的是“陡坡”的点与“缓坡”的点，其右邻域内绝大多数点的平均纵坐标通常会与中心点相差不大，而中心点的纵坐标都远小于左邻域内点的平均纵坐标，这时M取得最大值。由此可知最大M所对应的灰度值即为坡度从大到小的跳变位置i*，也即最佳阈值T：

(11)

以图3(b)中的直方图为例，对其进行从大到小的坡度跳变点检测后的结果如图4(a)所示。这里带宽取h=27，且由于SAR图像没有绝对精确的最佳阈值，位于坡度跳变点附近的位置都可以认为是理论上的最佳阈值，因此可以通过改变带宽来调整最佳阈值。从图4(a)可知，最大M值所对应的横坐标(估计的阈值)与理论上的最佳阈值还是很接近的。用估计的阈值对SAR图像进行海陆分割的结果如图4(b)所示，可以看出陆地和海面都存在少量未被正确分割的像素(即杂散点)，这些未被正确分割的像素就对应于前文所述的陆地中某些相对较暗的像素与海面中某些相对较亮的像素，这个现象同样证明了估计的阈值与理论阈值相契合。

(a) 坡度跳变检测后的M值

(b) 经阈值分割后的结果图4 最佳阈值估计与海陆分割结果

2.2 分割结果后处理

为了去除分割后产生的杂散点，需要借助形态学中寻找最大连通域的方法。首先，将连通标准设为4邻域连通，然后寻找到值为1的最大连通域后将其保留，那么值为1的非最大连通域则被自动去除。对应到二值图像中：大片的白色陆地区域被保留了，而海面中白色的杂散点则被去除了。接下来为了去除陆地区域中的黑色杂散点，需要将图像取反，然后再次保留最大连通域，然后再次将图像取反以还原视觉效果。如果图像中存在多片陆地区域，那么则需标定所有的连通域并设置阈值，大于一定阈值的连通域才能被保留，这样则可以保留多片陆地区域且去除杂散点。

至此，虽然已经获得了理想的分割结果，但理想的海岸线目前还不能在该结果中直接提取。这主要是因为经初始阈值分割后海陆交界处有许多的“毛刺”以及不规则形状出现。为此本文设计了一种用于消除“毛刺”及不规则形状的滤波手段。设H为去除杂散点后的二值分割结果，则消除过程可表示为

(12)

Hf={Hf(x,y):(x,y)∈H}

(13)

式中，N(xi,yi)表示以H或Hf中值为1的像素(xi,yi)为中心的3×3邻域，s表示N(xi,yi)中值为1的像素数目，N0(xi,yi)表示N(xi,yi)中所有的像素的值都为0。将式(12)的操作遍历完H后得到式(13)中的Hf，再将式(12)遍历Hf后得到新的Hf，以此重复迭代，直到Hf不再发生变化为止。

当N(xi,yi)经过海陆交界的规则形状处时，会使其中超过半数的像素都是值为1的像素，故此时需保留；反之，当N(xi,yi)经过海陆交界的不规则形状处时则会将少于半数的值为1的像素的值都置为0，起到填充孔洞从而将其改变为规则形状的作用。

图5为消除不规则形状示意图，图5(a)为图4(b)去除杂散点后的二值分割结果H；图5(b)为图5(a)中蓝色方框所在区域的局部放大，且其中的红框表示N(xi,yi)正位于不规则形状处，而绿框表示N(xi,yi)正位于规则形状处；图5(c)为图5(b)去除了不规则形状后的结果，且为图5(d)中橙色方框的局部放大；图5(d)为图5(a)去除了不规则形状后的结果。至此已经得到了最终的分割结果，然后再利用边缘检测算子就可以得到最终的海岸线。

(a) 去除杂散点后的结果 (b) 去除不规则形状

(c) 去除不规则形状后的结果 (d) 最终分割结果图5 去除不规则形状

3 实验结果与讨论

为了证明提出方法的有效性，将其分别应用到了由RADARSAT-2和Sentinel-1A卫星获得的4幅不同特点的大尺度SAR强度图像中。第1个实验图像尺寸为2 066×2 066像素，其特点是陆地区域有较为均匀的纹理结构；第2个图像的尺寸为1 872×2 419像素，其特点为有许多山位于陆地上以致陆地区域忽明忽暗；第3个图像的尺寸为1 617×2 005像素，其陆地区域中有很多建筑物与山；第4个图像的尺寸为2 331×2 224像素，其陆地区域主要由山与城区所组成。图6为该4幅实验图像，且在4个实验中，带宽h皆为17。

(a) 实验图像1 (b) 实验图像2

(c) 实验图像3 (d) 实验图像4图6 4幅大尺度实验SAR图像

图7(a)～(e)分别为本文方法应用到第1个至第4个图像后的直方图、最佳阈值估计结果、初始阈值分割结果、最终二值分割结果与最终提取的海岸线。为了目视效果更好，将提取的海岸线(红线)进行了加粗并叠加到了原图上。

为测试提出方法的准确性，采用Modava和Akbarizadeh[17]提出的一种基于邻域像素的评价标准。首先将手画的海岸线作为用于评价的正确海岸线，然后以正确海岸线为中心，9个像素为半径的范围作为评价区，计算提出方法所提取的海岸线像素分别落入不同半径评价区的累加百分比，精度评价结果如表1所示。

表1 提取的4幅SAR图像的海岸线像素分别落入不同半径评价区的累加百分比 %

(a) 直方图

(b) 最佳阈值估计

(c) 初始分割结果

(d) 最终分割结果

(e) 海岸线图7 4个实验的相关结果

由表1可知，本文方法对4幅大尺度SAR图像的提取精度均达到了90%以上，可以满足工程上的应用需求。

为了证明提出算法相对其他阈值分割法的有效性，将OTSU算法、陈祥等[19]提出的方法、Al Fugura等[15]提出的方法与Liu和Jezek[20]提出的方法作为对比方法来分割上述4幅大尺度SAR图像。其中依次用上述几个对比方法进行初始分割的结果与保留最大连通域后的最终分割结果分别如图8～图11所示。

从对比方法的结果可以看出，OTSU算法因为并非是针对海陆分割而设计的，因此陆地上的大量暗区也被误分为海域，导致后面的保留最大连通域的步骤中将大量的陆地区域都去除了，根本无法直接用来实现海陆分割；而陈祥等提出的方法由于是借助OTSU来确定最佳阈值的，故也无法较好地实现海陆分割；Al Fugura等提出的方法由于事先使用了Lee滤波操作来去噪，所以使得利用图像的灰度均值作阈值的分割结果好于前两种对比方法，但这无疑会降低原始图像的精度与最终提取的精度，且仍然无法避免陆地上暗区的错分现象；Liu和Jezek使用双峰高斯分布无法较好地拟合SAR强度图像中的分布，所以初始的与最终的分割结果也很不理想。因此直接利用这些阈值分割法无法获得理想的海岸线。

相比之下，本文提出的方法无须事先对图像进行滤波等预处理，也并不依赖于一种或几种特定的统计分布。而只是根据SAR强度图像的直方图中海域与陆地分布交界处有转折的特点，寻找到该转折处所对应的灰度值，因此更接近理论上的最佳阈值。最佳阈值分割结果中有少量的杂散点，结合形态学处理与本文设计的去除不规则形状方法就可以获得连续、光滑的理想海岸线。

(a) 初始分割

(a) 初始分割

(b) 最终分割图9 陈祥等提出的方法的初始分割结果与最终分割结果

(a) 初始分割

(b) 最终分割图10 Al Fugura等提出的方法的初始分割结果与最终分割结果

(a) 初始分割

(b) 最终分割图11 Liu和Jezek提出的方法的初始分割结果与最终分割结果

4 结束语

针对现有的基于阈值的SAR图像海岸线提取方法具有的难以确定最佳阈值的问题，本文首先介绍了基于差值核估计的跳变检测方法，然后分析了包含有海域和陆地的SAR强度图像的直方图分布特点。结合这种一维直方图独特的分布特点，再借助于跳变检测方法从而估计出了SAR强度图像海陆分割最佳阈值，进而利用该阈值完成了海陆初始分割并利用设计的后处理方法实现了较为理想的最终分割，从而得到了海岸线。最后以RADARSAT-2和Sentinel-1A卫星获得的SAR强度图像作为实验数据并进行了相关实验，实验结果表明，本文提出的方法可以快速有效地提取大尺度复杂SAR强度图像中的海岸线，且对实验结果进行的精度评价也显示提出方法的提取精度基本都达到了90%以上，从而验证了该算法在工程上的可行性与有效性。