一种SAR海岸图像中海岸线提取方法

摘要

一种SAR海岸图像中海岸线提取方法，属于图像处理领域；包括：获取SAR海岸图像；确定图像中初始海洋区域和初始陆地区域；计算初始海洋区域的几何中心；以初始海洋区域的几何中心为起点做射线；确定射线上的海岸边界点；将所有海岸边界点依次连接，得到海岸线；本发明利用G0分布对被测区域地表复杂程度敏感的特性，克服传统方法在进行海岸线提取时出现斑点噪声的问题；方法易于实现，运行速度快，且适用于大尺度图像。

说明书

技术领域

本发明属于图像处理领域，具体涉及一种SAR海岸图像中海岸线提取方法。

背景技术

海岸线是陆地和海洋的分界线，是海岸带的基本组成部分，同时也是划定海岸带范围的重要依据，由于人类活动和自然环境变化等因素的影响，海岸线不断变化，直接影响到海岸带的建设、管理和保护。因此，实现快速、准确、实时的海岸线监测对于临海国家的海岸带建设十分重要。合成孔径雷达(Synthetic Aperture Radar，SAR)系统不受气候条件及日照的影响，可以全天候全天时对海岸线情况进行监测，目前已经成为海岸线监测的有效工具。但由于SAR系统特殊的成像机制使得SAR海岸图像含有大量的斑点噪声，对实现SAR海岸图像中海岸线提取造成极大干扰。因此，针对SAR海岸图像进行海岸线提取是图像处理研究领域中的热点和难点。

目前，图像边缘检测方法被认为是SAR海岸图像中海岸线提取最为有效的方法，而以往的图像边缘检测方法大多是基于图像局部信息来进行提取的，如Canny算子、Sobel算子、Roberts算子、Prewitt算子和Laplacian算子等，虽然这些算子具有操作简单、运算速度快等优点，但是抗噪能力差，边缘定位不够准确，尤其是对含有大量斑点噪声的SAR图像，以上方法均不能取得较好的海岸线提取结果。经学者研究发现，G0分布是基于乘性噪声发展而来的统计分布模型，是目前SAR海岸图像数据建模的一个重要模型，利用G0分布对SAR海岸图像建模，可以有效的降低噪声对海岸线提取过程中的影响，并能够很好地描述SAR海岸图像中均匀、不均匀和极不均匀区域，对于海岸线提取有很大优势。

发明内容

针对上述现有技术存在的不足，本发明提供一种SAR海岸图像中海岸线提取方法。

本发明的技术方案：

一种SAR海岸图像中海岸线提取方法，包括以下步骤：

步骤1：获取SAR海岸图像；

步骤2：确定图像中初始海洋区域和初始陆地区域：

步骤2-1：将SAR海岸图像划分为大小相等的子块；

步骤2-2：根据G0分布求解每个子块的形状参数和尺度参数

$>\left(\begin{array}{c}\hat{\alpha }=-1+\frac{n\times {\hat{m}}_1}{n\times {\hat{m}}_1-(n+1/2){\hat{m}}_{1/2}^2}\\ \hat{\gamma }=(-\hat{\alpha }-1)\times {\hat{m}}_{1/2}\end{array}\right)$>

其中，n为等效视数，z_i为子块中的第i个像素点的灰度值，R为子块中像素点的个数；

步骤2-3：对设定阈值T，大于T的子块覆盖区域记为海洋，小于T的子块覆盖区域记为陆地，得到初始海洋区域和初始陆地区域；

步骤3：计算初始海洋区域的几何中心；

步骤4：以初始海洋区域的几何中心为起点，以θ为间隔角度向四周做射线，并将未经过初始陆地区域的射线删除；

步骤5：确定第p条射线上的海岸边界点：

步骤5-1：第p条射线所经过的像素点为其中j＝{1,2,…,m_p}，m_p为第p条射线所经过像素点个数，以及其周围的8个像素点作为子块，将该子块形状参数和尺度参数值作为像素点的形状参数和尺度参数

步骤5-2：求出第p条射线所经过的每个像素点的似然函数l(j)值：

$>l\left(j\right)=\underset{k=1}{\overset{j}{\Sigma }}ln\left(f_{G^0}\right(z_k^p;{\hat{\alpha }}_r^p,{\hat{\gamma }}_r^p\left)\right)+\underset{k=j+1}{\overset{m_p}{\Sigma }}ln\left(f_{G^0}\right(z_k^p;{\hat{\alpha }}_b^p,{\hat{\gamma }}_b^p\left)\right)$>

其中，表示第p条射线上所有像素点的值的集合，表示第p条射线上所有像素点的值的集合，为第p条射线上第k个像素点的灰度值，Γ()代表伽马函数；

步骤5-3：使|v_j1-v_j2|取最大值的第j个像素点为第p条射线上的海岸边界点，其中，

步骤6：重复执行步骤5，直到找出所有射线上的海岸边界点；

步骤7：将所有海岸边界点依次连接，得到海岸线。

有益效果：一种SAR海岸图像中海岸线提取方法与现有技术相比，具有如下优势：

(1)传统的边缘检测方法用于海岸线提取时难以克服SAR图像固有的斑点噪声，G0分布的参数对被测区域的地表复杂程度十分敏感，不仅适合对均匀区域、粗糙区域及极度粗糙区域进行精确建模，而且计算简单，是一种参数更少，更简单实用的统计分布模型；

(2)方法易于实现，运行速度快，且适用于大尺度图像。

附图说明

图1为本发明一种实施方式的SAR海岸图像中海岸线提取方法流程图；

图2为本发明一种实施方式的SAR海岸图像中海岸线提取方法的SAR海岸图像；

图3为本发明一种实施方式的SAR海岸图像中海岸线提取方法的海洋区域获取流程图；

图4为本发明一种实施方式的SAR海岸图像中海岸线提取方法子块划分示意图；

图5为本发明一种实施方式的SAR海岸图像中海岸线提取方法的初始海洋区域示意图；

图6为本发明一种实施方式的SAR海岸图像中海岸线提取方法的海岸边界点获取流程图；

图7为本发明一种实施方式的SAR海岸图像中海岸线提取结果图，其中(a)为SAR海岸图像中海岸线提取结果，(b)为海岸线提取结果和SAR海岸图像叠加图；

图8为本发明一种实施方式的基于Canny算子的SAR海岸图像中海岸线提取结果图，其中(a)为SAR海岸图像中海岸线提取结果，(b)为海岸线提取结果和SAR海岸线图像叠加图。

具体实施方式

下面结合附图对本发明的一种实施方式作详细说明。

如图1所示，本实施方式的一种SAR海岸图像中海岸线提取方法，包括以下步骤：

步骤1：获取SAR海岸图像，如图2所示；

步骤2：如图3所示，确定初始海洋区域和初始陆地区域：

步骤2-1：将待处理图像划分为大小相等的子块，如图4所示，将1024×1024的SAR海岸图像划分成大小为64×64的子块；

步骤2-2：根据G0分布求解每个子块的形状参数和尺度参数推导过程如下：

G0分布的概率密度函数为：

$>f_{G^0}(z,\alpha ,\gamma )=\frac{2n^n(n-\alpha )}{{\gamma }^{\alpha }\Gamma (-\alpha )\Gamma \left(n\right)}\frac{z^{2n-1}}{{(\gamma +{\mathrm{nz}}^2)}^{n-\alpha }}---\left(1\right)$>

其中，n为等效视数；α为形状参数；γ为尺度参数；Γ()代表伽马函数，z为像素点的灰度值；

式(1)的r阶矩表达式为：

$>m_{G^0}\left(Z^r\right)={\left(\frac{\gamma }{2n}\right)}^{r/2}\frac{\Gamma (-\alpha -r/2)}{\Gamma (-\alpha )}\frac{\Gamma (n+r/2)}{\Gamma \left(n\right)}---\left(2\right)$>

将r＝1和r＝1/2分别代入上式，得出1阶矩和1/2阶矩表达式为：

$>\left(\begin{array}{c}{m}_{1/2}={\left(\frac{\gamma }{n}\right)}^{1/4}\frac{\Gamma (-\alpha -1/4)\Gamma (n+1/4)}{\Gamma (-\alpha )\Gamma \left(n\right)}\\ m_1={\left(\frac{\gamma }{n}\right)}^{1/2}\frac{\Gamma (-\alpha -1/2)\Gamma (n+1/2)}{\Gamma (-\alpha )\Gamma \left(n\right)}\end{array}\right)---\left(3\right)$>

可推出：

$>\frac{m_{1/2}^2}{m_1}=\frac{{\Gamma }^2(-\alpha -1/4)}{\Gamma (-\alpha ){\Gamma }^2(-\alpha -1/2)}\frac{{\Gamma }^2(n+1/4)}{\Gamma \left(n\right){\Gamma }^2(n+1/2)}---\left(4\right)$>

为了求解α和γ，将1阶矩和1/2阶矩用1阶样本矩和1/2阶样本矩代替，即：

$>\frac{m_{1/2}^2}{m_1}=\frac{{\hat{m}}_{1/2}^2}{{\hat{m}}_1}---\left(5\right)$>

其中1/2阶样本矩和1阶样本矩分别为：

$>\left(\begin{array}{c}{\hat{m}}_{1/2}=\frac{1}{R}\underset{i=1}{\overset{R}{\Sigma }}{z}_i^{1/2}\\ {\hat{m}}_1=\frac{1}{R}\underset{i=1}{\overset{R}{\Sigma }}{z}_i\end{array}\right)---\left(6\right)$>

其中，z_i为子块中的第i个像素点的灰度值，即样本；R为子块中像素点的个数；

由上述公式得出子块的形状参数和尺度参数为：

$>\left(\begin{array}{c}\hat{\alpha }=-1+\frac{n\times {\hat{m}}_1}{n\times {\hat{m}}_1-(n+1/2){\hat{m}}_{1/2}^2}\\ \hat{\gamma }=(-\hat{\alpha }-1)\times {\hat{m}}_{1/2}\end{array}\right)---\left(7\right)$>

步骤2-3：对设定阈值T，大于T的子块覆盖区域记为海洋，小于T的子块覆盖区域记为陆地，得到初始海洋区域和初始陆地区域；

本实施方式中，阈值T为80；如图5所示，白色区域代表海洋，黑色区域代表陆地。

步骤3：计算初始海洋区域的几何中心；

本实施方式中，对初始海洋区域覆盖的所有像素点的横、纵坐标分别取平均值，并进行取整运算，得到初始海洋区域的几何中心C_h的坐标为(205，692)，如图5所示，白色海洋区域中的黑点即为海洋区域的几何中心。

步骤4：以初始海洋区域的几何中心为起点，以θ为间隔角度向四周做射线，第p条射线为Lp，其中p＝{1,2,…360/θ}，并将未经过初始陆地区域的射线删除；

本实施方式中，选取θ＝1°，第50-205条射线经过陆地区域。

步骤5：如图6所示，确定第p条射线上的海岸边界点：

步骤5-1：第p条射线所经过的像素点为其中j＝{1,2,…,m_p}，m_p为第p条射线所经过像素点个数，以及其周围的8个像素点作为子块，将该子块形状参数和尺度参数值作为像素点的形状参数和尺度参数

步骤5-2：求出第p条射线所经过的每个像素点的似然函数l(j)值：

$>l\left(j\right)=\underset{k=1}{\overset{j}{\Sigma }}ln\left(f_{G^0}\right(z_k^p;{\hat{\alpha }}_r^p,{\hat{\gamma }}_r^p\left)\right)+\underset{k=j+1}{\overset{m_p}{\Sigma }}ln\left(f_{G^0}\right(z_k^p;{\hat{\alpha }}_b^p,{\hat{\gamma }}_b^p\left)\right)$>

其中，表示第p条射线上所有像素点的值的集合，表示第p条射线上所有像素点的值的集合，为第p条射线上第k个像素点的灰度值，Γ()代表伽马函数；

步骤5-3：使|v_j1-v_j2|取最大值的第j个像素点即为第p条射线上的海岸边界点；似然函数l(j)单调递增，在海洋区域和陆地区域递增速度有所不同，在海洋区域增速慢，即相邻两点似然函数值差值较小，在陆地区域增速快，即相邻两点似然函数值差值较大，为第p条射线所经过第1个像素点的似然函数值l(1)到第j个像素点的似然函数值l(j)的平均增速；为第j个像素点的似然函数值l(j)到第m_p个像素点的似然函数值l(m_p)的平均增速。

步骤6：重复执行步骤5，直到找出所有射线上的海岸边界点；

步骤7：如图7(a)所示，将所有海岸边界点依次连接，得到海岸线；

海岸线和原SAR图像叠加结果如图7(b)所示。

利用现有Canny算子对SAR海岸图像提取海岸线，海岸线提取结果如图8(a)所示，海岸线和原SAR图像叠加结果如图8(b)所示。可以看出使用本发明方法得到的海岸线准确、连续，使用Canny算子得到的海岸线受陆地区域地物类型影响，有大量虚假海岸线，即噪声，且海岸线有不连续现象。