1 引言

复杂多变的海冰是极区最重要的自然现象之一,海冰及冰山在外部驱动力的作用下会发生漂移,直接影响极区船只的安全航行。北极海冰范围持续萎缩,大大提升了北极航道开通的可能性^[1]。北极由于其丰富的石油、天然气资源,备受各国瞩目^[2],在北极开展海洋运输以及近海活动已是当今的趋势^[3]。然而,我国在极区的研究资料有限,自2014年“雪龙”号的遇难事件以来,我国也在加强对极区海洋环境的监测^[4]。因此,提高对极区海冰运动的分析,不仅有重要的战略意义,也对我国乃至北半球的气候分析大有裨益^[5]。

极区由于其恶劣的环境条件,极大地限制了实测资料的获取,而早期的浮标观测和极地科考站都有其局限性,仅能观测小面积海冰漂移。然而,为了能够有效地观测海冰运动,模拟和预测海冰漂移,必须在一定的时空尺度上获取海冰的运动特征。随着卫星遥感技术的发展,这一问题得到有效地解决,目前最常用的卫星传感器包括微波辐射计(AVHRR)^[6]、散射计(NSCAT)^[7]、微波成像仪(SSM/I)^[7-8]、合成孔径雷达(SAR)^[9-10]等。其中,微波辐射计和散射计由于其分辨率的限制,常用于大范围的海冰漂移监测,无法得到高分辨率、重点区域的海冰漂移产品。合成孔径雷达具有全天时、全天候、多视角及对地物有一定的穿透能力等优点,近年来越来越多的国内外学者将其应用于海冰漂移检测。本文旨在选用SAR图像从空间上检测海冰的漂移,并对其进行量化。

本文提出一种基于三角剖分的特征跟踪与模式匹配结合算法,并将其有效地应用于高分辨率SAR图像海冰漂移检测。利用特征跟踪产生粗定位,通过建立三角网络,不仅为模式匹配提供搜索区域,而且可以自适应地确定模板窗口,较好地提高模式匹配算法的计算效率。该算法不仅能提高模式匹配的检测速度,而且不易受尺度变换、光照和角度的影响,其检测结果对图像噪声也有较高的稳健性,能进一步提高检测精度和均匀度,弥补了传统的特征跟踪算法在局部区域特征点缺失的不足。

2 材料与方法

2.1 卫星数据及研究区域介绍

Sentinel-1卫星是欧洲航天局2014年4月发射的对地观测卫星,其分辨率最高达5 m,幅宽达到400 km,支持双极化(HH+HV、VV+VH,其中H表示水平,V表示垂直),与其他卫星相比,Sentinel-1卫星数据不仅具有超高的辐射分辨率(1 dB/3σ,3σ表示辐射值的误差),有效地提升了雷达图像参数反演的精度,而且具有优良的覆盖和重访性能,能够在一天内覆盖整个欧洲和加拿大区域及周边海域^[11],这就为海冰漂移的实时检测提供了有力支持。

选用的Sentinel-1数据为超宽幅(EW)模式的地距影像(GRD),其为中等分辨率(40 m×40 m),该模式主要用于海上、冰川、极地等需要大范围覆盖和短重访周期的区域。相对于单视复数影像(SLC),GRD 数据有效地消除了热噪声并在一定程度上提高了图像质量^[12]。选择的图像对位于波弗特海域,覆盖区域如图1所示,其中,红色区域表示2017年10月10日的数据,绿色区域为2017年10月12日数据,虚线为本文截取的感兴趣区域,卫星数据如表1所示。

图 1. 地理信息图

Fig. 1. Geographic information map

下载图片 查看所有图片

2.2 算法介绍

2.2.1 Sentinel-1数据预处理

对于每个像素p,Sentinel-1数据Level-1产品的校准查找表(LUT)会提供一个数字量化值 DpN和一个归一化系数A_p,根据该归一化系数A_p,本文对图像进行了辐射校准,其定量关系为平方级^[13],校准公式为

vvalue(p)=|DNp|2Ap2。(1)

针对Sentinel-1图像存在的相干斑噪声,进行滤波处理。采用5×5的自适应滤波进行平滑处理,既在一定程度上抑制了斑点噪声,又保留了图像的细节信息,有效地提高了特征点检测的精度。

2.2.2 算法流程

基于特征跟踪和模式匹配相结合的算法,流程如图2所示。

图 2. 流程图

Fig. 2. Flow chart

下载图片 查看所有图片

1) 选择特征跟踪算法,对时间跨度为t的两幅图像进行初步检测。以尺度不变特征变换匹配(SIFT)算法为例,使用高斯差分算子检测极值点,其计算公式为

D(x,y,σ)=[G(x,y,kσ)-G(x,y,σ)]·I(x,y)=L(x,y,kσ)-L(x,y,σ),(2)

式中:k为相邻两层尺度相差的比例因子;变尺度高斯函数G(x,y,σ)= ×exp - ,其中m、n为高斯模板的维度;L(x,y,σ)为图像的高斯尺度空间,其中(x,y)为图像的像素位置,σ为尺度空间因子,其值越小表明图像的平滑度相对较小; I(x,y)为原图像。将每一个像素都与相邻点在图像域和尺度域进行比较,这样在二维空间和尺度空间都能检测极值点,且该方法对不同角度的卫星图像都有较高的适应性。同时,算法利用特征点邻域像素的梯度方向分布特征确定方向参数,使算子具备方向不变性。梯度的模值和方向的计算公式为

z(x,y)=[L(x+1,y)-L(x-1,y)]2+[L(x,y+1)-L(x,y-1)]2,(3)θ(x,y)=arctan{[L(x,y+1)-L(x,y-1)]/[L(x+1,y)-L(x-1,y)]}。(4)

然而,主方向的确定过于依赖局部区域像素的梯度方向,使得找到的主方向不准确,而特征向量提取以及匹配都严重依赖于主方向,即使微小的偏差角度也会造成特征匹配的较大误差,从而造成不可避免的误匹配。为了提高后续的检测效率及精度,本文对特征跟踪结果进行误匹配剔除。Thomas等^[14]的离群点剔除方法为数据正则化提供了一种有效的手段;同时还可利用一个中值滤波器剔除特征跟踪的异常矢量。

2) 针对剔除后的离散特征点,采用模式匹配算法,得到更均匀、规则的漂移矢量集。具体步骤为:

(1) 基于特征点的位置构造三角网络。为了得到最优三角网络,选择Bowyer-Watson算法进行二维delaunay三角剖分。由于delaunay三角剖分具有空圆特性和最大化最小角特性,实验得到的三角网络不仅具有唯一性,而且每一个网格都是趋于规则的,即相对于其他三角网络构造法,在delaunay三角网络中对每一个三角形的最小角进行升序排列得到的数值最大,这为模式匹配提供了更稳定的区域选择。图3所示为算法示意图,IMG₁和IMG₂分别为参考图像和待匹配图像,参考图像中的点A₁、B₁、C₁分别对应待匹配图像中的A₂、B₂、C₂,W为搜索窗口,点P为插入的预测点,通过在搜索窗口里进行模式匹配确定匹配点P'。

图 3. 算法示意图

Fig. 3. Diagram of algorithm

下载图片 查看所有图片

(2) 利用三角网络对特征点进行划分。为了有效逼近准确位置且减少不必要的运算,计算每一个三角形的面积,选择合适的阈值滤除高密集度的区域。

本文选用最大互相关(MCC)算法,确定最佳匹配区域。MCC算法的主要思想是匹配相关系数的计算,可表示为

ρ(x,y)=cov[f(x,y),g(x,y)]Dx,yD,(5)

式中:cov为两者的协方差;f(x,y)、g(x,y)分别为参考图像和待匹配图像的计算模板;D_x,y为f(x,y)的方差;D为g(x,y)的方差。由于所选择的卫星图像精度相对较高,因此对于样本窗口的选择需要遵循两方面的原则:样本窗口不能过大,以保证模板内海冰运动均一;样本窗口不能过小,以保证模板内有足够多的特征,并提高区分度。

海冰漂移的速度v可以表示为

v=dx2+dy2×40/(3600×t),(6)

式中:d_x和d_y分别为海冰在水平和垂直方向的偏移分量(以像素为单位),分辨率为40 m,d_x=x₂-x₁,d_y=y₂-y₁;(x₁,y₁)和(x₂,y₂)分别为时刻t₁和t₂的目标在两景影像中的位置;t为监测海冰漂移的时间间隔t=t₂-t₁。

3 实验与结果讨论

3.1 参数调整

特征跟踪算法剔除后有142个特征向量和249个剖分三角形,但部分三角区域较小,进行模式匹配有两点不足:1)匹配模板相对较小,限制了MCC算法的匹配精度;2)计算得到的漂移矢量与初始的漂移矢量临近,检测结果冗余。为了提高计算效率,选择面积大于300 pixel的区域进行模式匹配,同时为了保证算法的准确性与精度,还设置阈值为0.95的相似度,滤除匹配系数相对较低的矢量,在129组运算中,准确匹配的有126组,准确率高达98%。

3.2 实验

对HH极化和HV极化的Sentinel-1图像分别进行海冰漂移检测,如图4所示,其中HH极化的特征跟踪结果为164组,经过误匹配剔除后为142组,最终的漂移矢量为268个,而HV极化下的特征结果为234组,经过误匹配剔除后为215组,最终的漂移矢量为367个,约为HH极化的1.37倍。其结果与Komarov等^[15]对于Radarsat-2 ScanSAR的结论相吻合,相比于HH极化数据,HV极化数据在海冰漂移检测上具有更加可靠的检测效率,本文算法对两种极化方式的检测结果有明显的提高;相比于Muckenhuber等^[16]的结论(HV极化数据提供的海冰漂移矢量大约是HH极化数据的4倍),本文算法对HH极化数据的检测结果提升更明显,可见该算法适用性不受极化方式限制。

从主观角度评价,无论是HH极化还是HV极化下,特征跟踪的结果存在明显的信息缺失,而本文算法很好地弥补了这些缺失,有效地提高了海冰漂移检测的空间覆盖率。

本文的另一个目的是提高计算效率,特征跟踪的用时约为5.04 s,模式匹配的用时约为9.61 s,共用时14.65 s,而单独使用模式匹配算法,选取10个分布均匀的海冰样本,用时为125.34 s,算法时间缩短了约88%,为海冰漂移检测的实时性研究提供了有力的保障。

图 4. 海冰漂移矢量图

Fig. 4. Sea-ice drift vectors

下载图片 查看所有图片

3.3 结果讨论

为了验证本文测出海冰漂移速度的准确性,根据参考图像和待匹配图像的时间信息,选择相同时间对应的浮标数据,使用国际北极浮标计划(IABP)的浮标数据集进行标定和验证。通过浮标的经纬信息,确定其速度,与本文的海冰漂移结果相比较,并计算误差。由于浮标数量有限,本文针对感兴趣区域寻找最近的浮标(表2)。

将提取的海冰漂移速度用散点图的形式表示。图5为SIFT算法与本文算法的海冰漂移速度散点图,可以看出本文算法除了有少数小于0.080 m·s^-1,大部分都在0.080~0.085 m·s^-1之间,与浮标速度相一致。本文又比较了两种算法的均方根误差(RMSE),其中SIFT算法的RMSE为0.0036,本文算法的RMSE为0.0031。结果表明,本文算法具有较好的稳健性和较高的检测精度。

图 5. 海冰漂移速度散点图。(a) SIFT;(b) SIFT+MCC

Fig. 5. Scatter diagrams of sea-ice drift speed. (a) SIFT; (b) SIFT+MCC

下载图片 查看所有图片

从客观角度评价,通过盒状图对本文结果进行分析,如图6所示。本文算法的离群值明显降低,基本控制在40 pixel内,其中HH极化的降低效果相对更好,箱体明显减小,漂移矢量分布更加均匀,而下4分位数和平均值的略微提高并非与实验目的相违背,本文的参数设置也在刻意避免矢量的相互重合,从而提高了海冰漂移结果的可观性。

通过Voronoi图比较特征点的空间密度,如图7所示。颜色映射带由蓝到红分别表示由小到大的区域块,图7(a)和 (c)分别显示了SIFT算法在HH、HV极化下的特征跟踪结果的空间分布,而图7(b)、7(d)分别显示了本文算法的空间分布。可以明显看出对于两种极化方式,本文算法的检测结果具有更好的空间分布,仅在边缘区域存在较大缺失,其原因为目标区域海冰在单位时间跨度内,已经漂离待检测图像的区域。对于大部分区域,本文算法提供了更高的空间覆盖率,即更完善的漂移数据集。

图 6. 矢量距离盒状图。(a) HH极化;(b) HV极化

Fig. 6. Box charts of vector distance. (a) HH polarization; (b) HV polarization

下载图片 查看所有图片

图 7. 空间密度图。(a)(b) HH极化;(c)(d) HV极化

Fig. 7. Spatial density diagrams. (a)(b) HH polarization; (c)(d) HV polarization

下载图片 查看所有图片

4 结论

基于三角网络,提出一种特征跟踪和模式匹配相结合的新算法。在双极化Sentinel-1图像的海冰漂移检测应用中,本文算法表现出较高的空间覆盖率和较好的均匀性,相比模式匹配算法,算法时间降低了约88%,有效提高了检测效率。本文算法的优越性也在一定程度上依赖于预处理,其中辐射校准、热噪声处理以及利用自适应滤波器降低相干斑噪声对图像的干扰,在一定程度上提高了特征识别的效率。为了提高海冰漂移结果的准确率,没有直接在特征跟踪结果上进行模式匹配,而是采用有效的剔除算法,对离群点进行剔除,降低了误匹配率。

通过建立delaunay三角网络将两种算法进行结合,有以下优势:

1) 相较于简单的插值,小面积网格的滤除,使模式匹配结果更加均匀;

2) delaunay三角网络具有唯一性和规则性,为模式匹配提供了稳定可靠的搜索空间;

3) 利用网格区域大小自适应地确定模式匹配的样本窗口与搜索范围,提高了算法的稳健性。

本文算法对于HH极化和HV极化都有良好的检测效果,其检测结果也与浮标数据保持较高一致性,相对于传统特征跟踪方法,其误差更小,对噪声稳健性更高,可以有效地应用于海冰漂移检测的定量分析。