朱振宇,周乃恩,贺少帅
(1.航天彩虹无人机股份有限公司,北京 100074;2.青岛海洋科学与技术国家实验室发展中心,山东 青岛 266235)
合成孔径雷达(synthetic aperture radar,SAR)是一种全天时、全天候、广覆盖、高重访、中高分辨率的主动成像传感器[1]。随着观测角、频率、极化、空间和时间分辨率的不断提升,卫星SAR在满足陆海、军民需求上将会发挥越来越多的应用价值[2]。海上固定设施的金属、混凝土成分和动态船只的金属成分的后向散射系数远大于海洋的后向散射系数,因此在海洋SAR幅度图上海洋像素的灰度值一般很小,而海上固定设施和动态船只的灰度值一般很大,因此卫星SAR具有海洋目标识别、监测的工程化应用潜力。应用价值主要体现在两个方面:一是动态船只识别可以用于海洋伏季休渔期非法渔船捕捞监督、非法移民监控、海洋救援等民用领域;二是海上固定目标和动态船只的识别可以用于获取海上军事设施的位置,敌方舰船的位置、类型、航向等情报信息,夺取海上战场主动权[3]。
二十世纪八九十年代,加拿大、美国、挪威科研人员使用SAR卫星(如Seasat-1、ERS-1)成功提取了海上船只目标[4],而后国内外科研人员陆续开展了海洋SAR卫星数据的船只识别算法研究。如基于K分布的恒虚警概率(constant false alarm rate,CFAR)船只识别算法因简单易行、适应性强而被广泛应用[5],现已被用于加拿大海洋监测工作站[6]。随着SAR成像及其数据处理技术的发展,SAR影像海上固定目标和动态船只目标的识别研究得到了科研人员的关注[7-9]。针对不同应用场景和SAR数据特点,科研人员对CFAR算法进行了改进,其中双参数CFAR算法在良好海况下能够自动适应海杂波背景噪声的变化且算法简单、实用效果好[10-12],成为SAR目标检测中最为常用的算法。
当遇到台风、洋流、潮汐等较为复杂的海况时,SAR影像将出现局部海杂波,影响海上目标的检测、识别和提取的精度。另外,以往的研究大都是使用单期SAR单极化数据的幅度信息或者多极化数据的极化和幅度信息,通过开发、改进算法将船只与海杂波区分开,但未考虑海上固定设施的影响,尤其海上风力发电场的大小与船只接近的问题。本研究利用多期SAR数据的干涉相干性和幅度信息,使用双参数CFAR算法和众数滤波算法[13],分别开展海上固定设施和动态船只识别、提取的工程化应用研究。
选取江苏省南通市如东县近海区(位于中国黄海)作为研究区,地理坐标为121°04′ E~121°37′ E,32°26 N~32°59′ N,如图1中绿色矩形框所示。研究区内建有风力发电场、码头等固定设施和动态船只。在研究区内,划定了两个重点分析区,“重点分析区1”是近海风力发电场,可以看到成排成列分布的风力发电设施;“重点分析区2”是近海码头,码头附近可以看到船只,如图1红色矩形框所示。
实验辅助数据是日本基金会最新发布的15″格网(450 m)全球海陆高程数据(the GEBCO_2020 grid),作为数字高程模型(digital elevation model,DEM)数据,用于5期Sentinel-1A数据之间的精配准和地理编码。
总体技术方法流程如图2所示。使用5期Sentinel-1A SAR数据之间的干涉相干性和双参数CFAR算法与众数滤波算法(四邻域和八邻域窗口)提取近海海上固定设施。将海上固定设施、海杂波强反射区进行掩膜,使用双参数CFAR算法提取动态船只,并与未考虑掩膜信息的双参数CFAR算法、迭代加权多元变化检测算法提取的动态船只结果进行对比。
1)SAR影像配准和相干计算。同一地区多期重访SAR影像之间一般存在10~102量级像素的位置偏移。SAR干涉处理的必要条件是主辅SAR影像的精配准优于0.25个像素。一般可以使用卫星轨道参数估算主辅SAR影像之间的初始偏移量进行粗配准。在此基础上,本研究使用主辅影像相关性测度的相干系数法选取主辅影像的同名像点[14],使用最小二乘法分别在距离向和方位向拟合同名像点的偏移量,求得双线性多项式系数,将辅影像上各像素分别重采样到主影像坐标空间。
假设M和S分别是精配准的主、辅影像,选定大小为k×l的滑动窗口,计算滑动窗口内像素相干系数的平均值作为窗口中心像素的相干系数值[15]。
2)双参数CFAR算法。双参数CFAR算法应用在卫星SAR海上目标识别方面是海杂波(如海洋内波、洋流、海浪等引起的强反射)和海上目标分别服从不同的幅度统计分布概率密度函数。设Pfa为虚警概率,表示海上目标检测的错误概率,那么检测阈值T的表示如式(1)所示。
(1)
式中:μB和σB分别为局部滑动窗口幅度值的均值和标准差。在Pfa下,当幅度值I大于阈值T时,则该像素判定为海上目标;否则,判定为海面。
双参数CFAR算法适用条件有两个:一是海上目标的SAR像素具有较大和相近的幅度值,具有连续的空间分布;二是海杂波的SAR像素在幅度值和空间分布上是随机的、无规律的。
3)Frost滤波和众数滤波算法。Frost滤波算法是目前最为通用的基于局部统计特性的自适应滤波算法。Frost滤波算法假设斑点噪声是乘性噪声,并假设SAR影像是平稳过程。Frost滤波的冲击响应为双边指数函数[16]。
在海杂波较为极端复杂时,双参数CFAR算法提取的海上目标虚警率较高。因此,本实验对双参数CFAR算法获取的二值图进行众数滤波处理。在GIS应用中,众数滤波是指某一像素及其邻域像素的所有像素值组成了一个离散数值集合{x1,x2,…,xn},那么,出现频率最高的数值xi为该像素的像元值。本实验离散数值集合为{0,1},Sentinel-1A多视后空间分辨率为15 m,海上船只大小一般在2个像素及以上,因此滑动窗口选定四邻域和八邻域。
4)迭代加权多元检测(iteratively reweighted MAD,IR-MAD)算法。为了改进多元变化检测(multivariate alteration detection,MAD)算法[17]的去噪效果和实现变化阈值的自动获取,Nielsen[18]提出了IR-MAD算法。IR-MAD算法通过迭代加权计算自动获取变化阈值,是一种非监督变化检测方法。
SAR像对相干性受多种因素影响,如时空基线、传感器工作频率、天气和视场环境等。时空基线越小,传感器工作波长越大,则相干性越高。云雾雨天气以及视场内植被、河湖等地物会降低相干性,而稳定的建筑物则会提高相干性。因此,如果数据量丰富,则可通过控制时空基线进行SAR像对组合。考虑成像时间应避开海洋伏季休渔期和结冰期,本研究选出2020年9月27日—2020年11月14日的5期Sentinel-1A影像,数据量适中。考虑相干性受成像环境(大气、海况等)影响,因此不再控制时空基线优选干涉像对,而是对这5期SAR影像进行两两组合,共形成10个干涉像对。
分别对5期覆盖研究区的SAR影像进行感兴趣区(region of interest,ROI)裁剪、影像粗配准和精配准、干涉组合和相干计算、双参数CFAR算法处理、众数滤波、多视处理、海杂波强反射边界提取等数据处理。
首先提取近海海上固定设施。裁剪后的5期VV、VH双极化SAR影像两两组合,共可以形成10个VV和10个VH极化干涉对。对20个干涉对组合进行精配准和相干计算。使用双参数CFAR算法提取具有强干涉相干性的固定设施,得到二值化相干图。使用四邻域和八邻域众数滤波器交互式迭代处理二值化相干图,直到固定设施像素的变化个数小于固定设施像素总个数的1%时停止,则得到固定设施提取结果。图3和图4第一列为VV极
化干涉对提取的固定设施,第二列为VH极化干涉对提取的固定设施,第三列为双极化干涉对提取的固定设施,即对VV和VH极化干涉对提取的固定设施取并集;白色像素代表海上固定设施,黑色像素代表海域。对10个双极化干涉对提取的固定设施矢量边界取并集,得到“总干涉对”提取的固定设施结果,如图5所示。
之后提取动态船只目标。对裁剪后5期SAR SLC影像进行方位向1×距离向4的多视处理,得到空间分辨率为方位向13.95 m×距离向14.72 m的VV、VH双极化SAR幅度图。对5期幅度图进行Frost滤波处理,滤波窗口大小为7×7,得到降噪后幅度图。使用GEBCO_2020数据将SAR坐标系的降噪后幅度图地理编码到WGS1984 UTM坐标系。
IR-MAD算法提取动态船只。由于IR-MAD算法提取的结果是两期影像的变化目标,不考虑噪声情况下,则是两期的动态船只。因此在多时相SAR组合设计上,首先,使用第1期影像分别与第2~5期影像进行变化检测,分别得到“1期+2期”“1期+3期”“1期+4期”和“1期+5期”动态船只;然后,使用第2期影像与第3期影像进行变化检测,得到“2期+3期”动态船只;接着,使用OTSU算法和众数滤波对变化检测结果进行自动阈值分割和滤波处理,提取矢量边界;最后,用“1期+2期”矢量边界对“2期+3期”矢量边界做擦除(erase),得到第1期动态船只,分别用“1期+2期”“1期+3期”“1期+4期”和“1期+5期”矢量边界对第1期动态船只矢量边界做擦除,得到第2期、第3期、第4期、第5期动态船只。IR-MAD算法最大循环次数设置为30,停止迭代的最小变化率阈值设置为0.001。分别对每期VV、VH幅度图进行以上处理,再对提取的结果取并集作为每期动态船只,如图6所示。
双参数CFAR算法提取动态船只目标。设置船只大小为2~134个像素,直接使用双参数CFAR算法和众数滤波从5期UTM坐标系VV、VH幅度图提取船只,对每期VV、VH幅度图提取的船只目标取并集,结果如图7第二列所示,而第一列表示VV极化幅度图。将固定设施和人工提取的海杂波强反射区进行掩膜,设置船只大小为2~134个像素,使用双参数CFAR算法和众数滤波从5期UTM坐标系VV、VH幅度图提取船只,对每期VV、VH幅度图提取的船只取并集,结果如图7第三列所示。
图3和图4横向对比,所有干涉对VV、VH提取的固定设施可分为三类:第一类是VV提取,而VH未提取;第二类是VH提取,而VV未提取;第三类是VV、VH同时提取,该类数量远大于第一类和第二类数量之和。实验结果表明:海上固定设施的结构和所处环境会共同决定其对不同极化方式雷达波的后向散射特性;同极化(VV)与交叉极化(VH)SAR数据提取的海上固定设施可互为验证和补充,因此在海上固定设施检测、识别方面极化信息具有重要价值。图3和图4纵向对比,不同时相组合的干涉对提取的固定设施在数量和分布上存在差别,这是因为不同时相组合的影像对在雷达成像时受到海况、天气和时空基线的影响不同。
如图8所示,在“重点分析区1”内,海上风力发电场为固定设施,在时间序列上能保持较好稳定性,因此在任意两期影像之间具有较高相干性。将SAR干涉相干性提取的94个近海风力发电场与2020年天地图影像进行比对,风力发电场被全部提取(如图8红色轮廓所示),错提个数为9个(如图8蓝色轮廓所示)。因此,SAR干涉相干性提取近海风力发电场的正确率达到90.4%,虚警率为9.6%,漏检率为0。
由图7第一列可知,5期SAR影像VV幅度图都存在不同程度的海杂波,VH幅度图亦然,因篇幅原因,未在文中展示。由图7第二列可知,直接使用双参数CFAR算法和众数滤波从VV和VH幅度图提取了动态船只并取并集,得到的双极化动态船只结果存在较大虚警率,主要原因有两个:海上固定设施和动态船只都表现出了较强的后向散射;复杂海况引起的海杂波。由于动态船只在时间序列上位置不定,所以其任意两期SAR的相干理论值为零,因此可通过使用SAR幅度和相干性提取近海动态船只。将海杂波强反射区和SAR干涉相干提取的固定设施掩膜,使用双参数CFAR算法和众数滤波从VV和VH幅度图提取了动态船只并取并集,如图7第三列所示。SAR幅度图提取的动态船只基本消除了海上固定设施和海杂波强反射的干扰。
图6是直接对多时相SAR双极化幅度图使用IR-MAD算法变化检测算法提取的船只目标。对比图6和图7可知,IR-MAD算法提取的船只目标存在较大虚警率,这主要是因为受到海杂波强反射的影响。
在“重点分析区2”内,为了有效分析输入掩膜的双参数CFAR算法和IR-MAD算法提取的动态船只目标的精度,本文使用船舶自动识别系统(automatic identification system,AIS)数据进行对比、统计。通常船只会安装AIS系统,通过该系统船舶能够在公用无线信道上向附近船舶和岸上主管机关持续发送其身份、位置、航向、航速等信息。本文所用AIS数据为船队在线HiFleet提供,包括MMSI(maritime mobile service identify)数字码、接收时间、船只航向和方位、经纬度坐标、船只类型等信息。2020年11月14日24小时的AIS数据与“重点分析区2及周边区域”叠加展示如图9所示。SAR干涉相干提取的近海海上码头设施和IR-MAD算法、输入掩膜的双参数CFAR算法提取的码头附近船只如图10所示,其中红色轮廓为SAR干涉相干提取的海上码头设施,白色轮廓为船只目标。
分别提取2020年9月27日、10月9日、10月21日、11月2日和11月14日的9点55分时刻的AIS数据与本文算法所提取的5期船只目标比对分析,如图10所示,图中白色不规则多边形为本文算法提取的船只目标,黄色点为AIS数据显示船只。通过统计分析可知,输入掩膜的双参数CFAR算法和IR-MAD算法都存在一定的漏检现象,输入掩膜的双参数CFAR算法平均漏检率为3.34%,IR-MAD算法平均漏检率为7.65%,IR-MAD算法的虚警率为16.39%,远高于输入掩膜的双参数CFAR算法2.29%的虚警率。由此可知,输入掩膜的双参数CFAR算法提取的船只目标结果优于IR-MAD算法变化检测算法提取结果,详见表1。
表1 输入掩膜的CFAR算法、IR-MAD算法提取船只统计分析
针对近海海上固定设施和动态船只较难通过SAR幅度信息有效区分的问题,尤其是海上风力发电场的大小与船只接近,本研究以江苏省南通市如东县近海区(位于中国黄海)作为实验区。首先,利用5期SAR数据之间的干涉相干性,采用双参数CFAR算法和众数滤波算法,有效提取了海上固定设施,如风力发电场和码头等。之后,人工提取VV、VH极化幅度图出现的海杂波强反射区,将海上固定设施和海杂波强反射区掩膜,使用双参数CFAR算法和众数滤波器从VV和VH极化强度图提取了动态船只并取并集,得到了5期SAR数据海上船只,并与未考虑使用掩膜信息的双参数CFAR算法和众数滤波器直接提取的海上船只进行了对比。最后,分析和讨论了本研究实验结果,在“重点分析区1”,SAR干涉相干提取的近海风力发电场的正确率为90.4%,虚警率为9.6%,漏检率为0;在“重点分析区2”,给出了码头附近5期AIS数据,通过与输入掩膜的CFAR算法、IR-MAD变化检测算法的动态船只提取结果对比可知,输入掩膜的CFAR算法提取的船只平均漏检率为3.34%,平均虚警率为2.29%,IR-MAD算法变化检测算法平均漏检率为7.65%,平均虚警率为16.39%,输入掩膜的CFAR算法提取的船只目标结果优于IR-MAD算法变化检测算法。
本实验所用总体技术方法适用于海域场景的海上固定设施和动态船只的识别,同时为陆域场景的固定设施(如建筑物)和动态目标(如车辆、飞机等)的识别研究提供参考。后续研究重点,计划开展SAR极化分解模型的海杂波强反射区提取研究,从而实现近海海上固定目标(如风力发电场、码头等)和动态船只识别和提取的自动化、工程化应用。
致谢:感谢欧洲航天局提供的Sentinel-1A数据,感谢日本基金会提供的全球海陆高程数据,感谢“天地图·江苏”提供的2020年光学影像。