视频合成孔径雷达双域联合运动目标检测方法

丁金闪 仲超 温利武 徐众

王志豪, 李刚, 蒋骁. 基于光学和SAR遥感图像融合的洪灾区域检测方法[J]. 雷达学报, 2020, 9(3): 539–553. doi: 10.12000/JR19095
引用本文: 丁金闪, 仲超, 温利武, 等. 视频合成孔径雷达双域联合运动目标检测方法[J]. 雷达学报, 2022, 11(3): 313–323. doi: 10.12000/JR22036
WANG Zhihao, LI Gang, and JIANG Xiao. Flooded area detection method based on fusion of optical and SAR remote sensing images[J]. Journal of Radars, 2020, 9(3): 539–553. doi: 10.12000/JR19095
Citation: DING Jinshan, ZHONG Chao, WEN Liwu, et al. Joint detection of moving target in video synthetic aperture radar[J]. Journal of Radars, 2022, 11(3): 313–323. doi: 10.12000/JR22036

视频合成孔径雷达双域联合运动目标检测方法

DOI: 10.12000/JR22036
基金项目: 国家自然科学基金(62171358)
详细信息
    作者简介:

    丁金闪(1980-),男,江苏人,博士,教授,研究方向为视频雷达系统及信号处理技术、新体制雷达等

    通讯作者:

    丁金闪 ding@xidian.edu.cn

  • 责任主编:张润宁 Corresponding Editor: ZHANG Running
  • 中图分类号: TN957

Joint Detection of Moving Target in Video Synthetic Aperture Radar

Funds: The National Natural Science Foundation of China (62171358)
More Information
  • 摘要: 视频合成孔径雷达(SAR)具有高帧率成像能力,可作为地面运动目标探测的重要技术手段。经典SAR地面动目标显示(SAR-GMTI)依靠目标回波能量来实现动目标检测,同时动目标阴影亦可作为视频SAR动目标检测的重要途径。然而,由于动目标能量和阴影的畸变或涂抹,依靠单一方式难以实现稳健的动目标检测。该文基于目标能量和阴影的双域联合检测思想,分别通过快速区域卷积神经网络和航迹关联两种技术途径实现了视频SAR动目标联合检测,给出了机载实测数据处理结果,并进行了详细分析。该文方法充分利用目标阴影与能量的特征及空时信息,提升了机动目标检测的稳健性。

     

  • 遥感技术为地理勘探[1]、建筑测绘[2]、植被监测[3]、水域监测[4]、灾区探测[5]等提供了重要的技术支持。在洪灾区域检测中,光学遥感能够提供丰富的光谱信息和较高的图像分辨率,而雷达遥感则具备全天候的洪灾区域检测能力[6]。异质图像融合已经在复Contourlet域得到了实现[7],因此结合光学遥感数据和雷达遥感数据能实现洪灾区域的全天候、高时效检测[8]

    当前,学者们提出许多基于遥感图像的洪灾区域检测方法。其中,应用人工神经网络和自组织映射网络[9]方法处理洪灾前后的同质SAR图像时能大致区分洪灾区域,但是受到大量相干斑噪声的影响,检测结果仍然存在很高的虚警率。为了结合光学遥感图像和SAR图像的成像优势,一些学者提出了基于异质图像的洪灾区域检测方法,例如像素变换法[10]和特征空间变换法[11]。上述方法对洪灾前的SAR图像和洪灾后的光学图像有良好的检测效果,但是对洪灾前的光学图像和洪灾后的SAR图像的检测效果较差。

    针对SAR图像中洪灾区域的检测难点,本文提出了H-FCM算法。该算法的主要贡献如下:

    (1) 提出SAR图像的分级聚类模型;

    (2) 将洪灾前的近红外波段光学遥感图像中提取出的河流位置融合到洪灾后的SAR图像中,并作为洪灾区域检测的空间约束,进一步提升了检测结果的准确率。

    文章结构如下:第2节分析了近年来洪灾区域检测方法的发展和存在的问题。第3节介绍了分级C均值聚类算法Hierarchical Fuzzy C-Means (H-FCM)算法的原理、步骤以及参数设定对算法性能的影响。第4节应用包括H-FCM算法的4种算法对不同类型遥感数据进行洪灾区域检测,并分析了图像降噪处理对4种算法检测性能的影响。第5节系统总结了全文的主要工作和存在的不足。

    目前,学者们提出的基于雷达遥感图像的水体检测方法包括灰度阈值分割法、滤波法、机器学习法和结合辅助信息的提取方法等[12];基于光学遥感图像的提取水体方法包括基于像元分类的阈值法和基于目标分类的分类法等[12]。洪灾是由于大雨、暴雨或者持续降雨造成的,此时无法通过实时获取的光学遥感数据检测地面的洪灾区域。

    以2019年8月9日对利奇马台风实时追踪的GF4数据(见图1)为例,光学图像无法显示地表的洪灾区域。但是通过适当的卫星调度方案可以短时间内获取雷达遥感数据。以2019年8月8日对利奇马台风实时追踪的GF3数据(见图2)为例,SAR图像可以清晰地显示地面受灾情况。

    图  1  利奇马台风的GF4图像(资源卫星中心)
    Figure  1.  The GF4 image of Typhoon Lekima (Resource Satellite Centre)
    图  2  利奇马台风的GF3图像(资源卫星中心)
    Figure  2.  The GF3 image of Typhoon Lekima (Resource Satellite Centre)

    应用阈值分割法中的最大类间方差算法[13]对灾后的SAR图像进行洪灾区域检测时,由于受到相干斑噪声的影响,检测精度较低。为克服相干斑噪声的影响,学者们提出了滤波法[14]、边界追踪法[15]、Markov分割法[16]、基于邻域最小生成树的半监督分类法[17]、基于概率转移卷积神经网络的分类法[18]等。但是洪灾区域检测问题的核心在于保证高检测准确率的同时,尽可能缩短检测时间,而上述方法的计算复杂度普遍较高,并不能满足洪灾区域检测的高时效性需求。

    直接利用K-means算法[19]、C均值聚类方法(FCM)算法[20]等无监督聚类算法对洪灾后的SAR图像进行洪灾区域检测可以缩短检测时间,但是检测结果的分类准确率很低。具体原因如下:根据灰度阈值分割法的原理,在SAR图像中,水体的散射值较低,整体呈现为暗区。在没有发生暴雨或者持续降雨的情况下,可以通过分析SAR图像的整体像素值分布,应用FCM算法等聚类算法确定水体的分割阈值。根据分割阈值,将大于阈值的标记为非水体,小于阈值的标记为水体。但是洪水灾害发生时伴随着暴雨或者持续降雨,导致很多未受到洪灾影响的区域在SAR图像中也呈现为暗区。因此,直接通过阈值法获得的检测结果存在较高的虚警率。

    FCM算法的步骤总结为表1。以洪灾后的大小为n1×n2像素的SAR图像为例,将其所有像素点的像素值记作二维矩阵X。将X分成c类,则可以得到c个聚类中心。将每一个像素点作为单个样本xj,则xj对于第i类聚类中心的隶属度可记作uij, FCM算法的目标函数、约束条件的定义式分别为式(1)和式(2)

    表  1  FCM算法
    Table  1.  FCM algorithm
     (1) 给定聚类中心数目c,隶属度因子m,迭代次数t
     (2) 初始化和为1的隶属度矩阵U
     (3) 由式(4)确定c个聚类中心Ci
     (4) 计算目标函数J的值;
     (5) 将Ci代入式(3),更新隶属度矩阵U
     (6)重复步骤(3)—步骤(5),迭代次数为t
    下载: 导出CSV 
    | 显示表格
    J=ci=1nj=1umijxjCi2
    (1)
    ci=1uij=1
    (2)

    其中,Ci表示第i类聚类中心;m表示隶属度因子;n=n1×n2; i=1,2,···,c; j=1,2,···,n

    结合式(2)的约束条件,利用拉格朗日乘数法,将J分别对uijCi求导,并令其结果等于0,可得uijCi的迭代式分别为式(3)和式(4)

    uij=1ck=1(xjCixjCk)2m1
    (3)
    Ci=nj=1(xjumij)nj=1umij
    (4)

    本文提出的H-FCM算法致力于缩短SAR图像洪灾区域的检测时间、降低洪灾区域的错误检测概率以及抑制SAR图像中存在的大量随机分布的相干斑噪声对检测结果的影响。根据河流连续体(River Continuum Concept, RCC)[21]的概念,在已知的河流系统内,自上游到下游的物理变量的梯度变化是连续的[22]。因此发生洪水灾害之后,若地理环境没有发生剧烈变化,洪水的形态仍然符合RCC。并且,根据典型地物的光谱特征曲线图[23](见图3),水体对近红外波段的电磁波吸收能力比大部分地物强。

    图  3  典型地物的光谱特征曲线[23]
    Figure  3.  Spectral characteristic curves of typical features[23]

    H-FCM算法的应用前提为可获得洪灾前的光学图像以及洪灾后的SAR图像。该算法通过对遥感图像的像素强度划分进行分级聚类。以uint8格式的图像数据为例,应用FCM算法获得待测数据像素值的8个聚类中心,并且由低至高排列。根据图3, SAR图像中代表水体的像素点的像素值主要分布在0~64之间,即处于低像素值区间。但是由于发生洪灾时,SAR图像中存在大量随机分布的相干斑噪声以及低洼地区存在的积水,直接用阈值法进行洪灾区域检测的准确率较低。

    H-FCM算法采用反向聚类的方法,考虑到在SAR图像中,代表水体的像素点几乎不可能处于高像素值区间。H-FCM算法在聚类过程中将处于最高级的像素值区间的像素点记为1,其余像素点记为0,对所有零像素点进行区域聚类,获得初步的洪灾区域检测结果。又由RCC原理,洪灾发生后,主要受灾区域与洪灾前的河流相连通,利用灾害前的河流位置作为空间约束可以得到近似真实洪灾区域的检测结果。因此,应用H-FCM算法检测洪灾发生后的SAR图像的洪灾区域,可以降低相干斑噪声对检测结果的影响。

    H-FCM算法主要分为4步。首先,应用FCM算法获得洪灾后SAR图像的聚类中心,进而设置分类阈值,并确定分级聚类模型;其次,提取灾害前近红外波段光学遥感数据中的河流位置,并将其融合到洪灾后的SAR图像中;再次,以聚类模型为基础进行分级聚类,获得洪灾区域的初步聚类结果;最后,通过空间约束矩阵对初步聚类结果进行约束,进一步提升洪灾区域检测的准确率。

    利用FCM算法获得洪灾后SAR图像的八分类的聚类中心矩阵C和阈值矩阵T,并由此获得7个聚类模型

    C=[C1C2···C8]
    (5)

    其中,一维矩阵C表示聚类中心矩阵,包含8个按照升序排列的聚类中心

    T=[T1T2···T7]
    (6)
    Tq=(Cq+Cq+1)2
    (7)

    其中,T0=0,q=1,2,···,7

    以1999年英国格洛斯特洪灾后获得的SAR图像Ea[10]为例,其大小为n1×n2(见图4)。遍历Ea中的所有像素点,若其像素值大于Tq–1且小于Tq则在对应位置标记为1,否则标记为0,将所得到的矩阵记作聚类模型Gmq,如式(8)

    图  4  洪灾后的SAR图像Ea[10]
    Figure  4.  SAR image Ea after the flood[10]
    (Gmq)mn={1,Tq1<(Ea)mn<Tq0,otherwise
    (8)

    其中,“(Gmq)mnGmq”在本文中均代表(Gmq)mnGmq矩阵中第m行,第n列的元素,下文符号“∈”的含义与此处一致;(Ea)mnEa; q=1,2,···,7

    3.2.1   SAR图像洪灾区域检测的改进方法

    由于SAR图像的成像模式一般为宽幅扫描模式,同时受到SAR图像中相干斑噪声的影响,在没有灾害前河流信息的情况下直接采用无监督聚类的方法来检测Ea中的洪灾区域存在较高的虚警率,导致检测结果的准确率较低。H-FCM算法将洪灾前河流的空间信息融合到洪灾后的SAR图像中,通过局部空间约束,可以提高洪灾区域检测的准确率。但是在图4中,Ea的河流位置有一部分被背景地物或者噪声覆盖,另一部分和洪灾区域重合,导致大部分的河流位置已经无法识别。为了获得洪灾前的河流位置信息,H-FCM算法进一步检测洪灾前的光学图像。

    以1999年英国格洛斯特洪灾前的SPOT图像中提取的近红外波段图像Ei[10]为例(见图5)。由于不同河段的河流泥沙含量不同,造成不同河段对近红外波段电磁波的吸收能力不同[24],通过阈值法无法直接获得所有河流区域的空间位置。考虑到洪灾前的水体表面近似为平滑面,可以应用区域生长算法[25]提取洪灾前近红外波段遥感图像的河流位置。本文提出的区域生长算法为表2,根据图5中近红外波段的光学遥感图像可以利用该算法获得洪灾前的河流位置。

    表  2  区分河流和道路的区域生长算法
    Table  2.  Region growing algorithm to distinguish rivers and roads
     (1) 应用FCM算法获得图像的八分类聚类中心C
     (2) 获得阈值矩阵T,并获得小于T2的二值图Pt
     (3) 通过卷积消除Pt中的区块噪声,并确定初始生长点;
     (4) 设定初始生长方向矩阵Jdi,i=1,2,3,4,如图6图9
     (5) 每次迭代都以初始生长点为中心,取5×5区域中的最小值点
       作为下一个生长点,按同一方向迭代直至生长到全图边界;
     (6) 将每个方向的生长结果累加,得到初步生长结果;
     (7) 利用霍夫变换[26]提取初步生长结果中的道路;
     (8) 将初步生长结果去除道路部分,获得河流位置。
    下载: 导出CSV 
    | 显示表格
    图  5  洪灾前的近红外波段影像Ei[10]
    Figure  5.  Near-infrared image Ei before the flood[10]
    3.2.2   检测河流位置的区域生长算法规则
    图  6  方向矩阵Jd1(向下)
    Figure  6.  Direction matrix Jd1(down)
    图  7  方向矩阵Jd2(向上)
    Figure  7.  Direction matrix Jd2(up)
    图  8  方向矩阵Jd3(向左)
    Figure  8.  Direction matrix Jd3(left)
    图  9  方向矩阵Jd4(向右)
    Figure  9.  Direction matrix Jd4(right)

    由于步骤(2)获得的第2类阈值图Pt(见图10)中存在大量的区块噪声的干扰,无法从中直接获得灾害前的具体河流的空间位置。因此,需要筛除大面积区块性的非河流区域。像素点的面积定义为:将像素值为1的像素点的面积记作1,像素值为0的像素点面积记作0。连通区域的面积定义为:将连通区域内所有非零像素点的面积累加得到该区域的面积。4连通区域的定义为目标像素点(红色方块)的上、下、左、右4个像素点位置(见图11)。

    图  10  第2类阈值图Pt
    Figure  10.  Threshold graph Pt
    图  11  4连通区域的定义
    Figure  11.  Definition of four-connected regions

    为了充分消除区块噪声的影响,并保留河流区域的完整性,步骤(3)的卷积操作分为两步。首先,确定大小为α×α的全1卷积核Wα,将PtWα卷积得到初步消除区块噪声的结果P1,如式(9)和式(10);其次,将P1Wβ卷积得到基本消除区块噪声的结果P2,如式(11)和式(12)。为保证卷积之后的输出结果的尺寸与原图像大小保持一致,本文所有的卷积模式均采取matlab中conv2函数的‘same’模式进行卷积

    Pt1=PtWα
    (9)

    其中,“*”表示以‘same’模式进行卷积

    (P1)ij={1,(Pt1)ij<0.6α20,otherwise
    (10)

    其中,(Pt1)ijPt1; (P1)ijP1

    Pt2=P1Wβ
    (11)

    其中,β=α/2, Wβ为大小为β×β的全1卷积核

    (P2)ij={1,(Pt2)ij<0.6β20,otherwise
    (12)

    其中,(P2)ijP2, (Pt2)ijPt2

    3.2.3   区域生长算法的参数选择

    参数α一般为固定的大小,当选择的α较大或者较小时,在筛除区块噪声的同时,也筛除了河流。本文固定参数α=18的卷积核Wα初步筛除大面积聚点,获得P1;其次,固定参数β=9的卷积核Wβ进一步筛除小面积聚点,获得P2。利用式(13)至式(15)确定河流区域种子生长点的位置Po,在Po的每一个4连通区域中等间距选择10个像素点作为初始生长点(见图12)。

    图  12  初始生长点位置
    Figure  12.  Initial growth point location
    Po=naj=1(detjObjj)
    (13)

    其中,Objj表示P2中的第j个4连通区域;P2中4连通区域总数为na;每个连通区域的背景大小均为n1×n2; j=1,2,···,na

    detj={1,Areaj>5000,Areaj500
    (14)
    Areaj=n1m=1n2q=1(Objj)mq
    (15)

    其中,Areaj表示P2中第j个4连通区域Objj的面积。

    3.2.4   利用区域生长算法获得河流位置

    根据图12中获得初始生长点,按照步骤(5)进行4个方向的区域生长,最后将各个方向的生长结果累加得到初步生长结果(见图13)。但是由于本例中道路和河流交错,当道路和河流宽度相近且没有先验测绘信息时,直接通过机器识别提取并区分河流和道路的空间位置是非常困难的。所以在步骤(7)中利用霍夫变换[26]检测初步生长结果中属于公路的长距离直线段区域,进而识别道路所在区域(见图14)。最后,将初步生长结果去除道路区域后,得到河流位置Eb(见图15)。

    图  13  初步生长结果
    Figure  13.  Preliminary growth results
    图  14  道路的识别结果
    Figure  14.  Road recognition results
    图  15  洪灾前的河流位置Eb
    Figure  15.  River location Eb before the flood
    3.2.5   灾后SAR图像与灾前河流位置融合

    将洪灾前的河流位置Eb和洪灾后的SAR图像Ea融合。融合规则为将Ea中与Eb中灾前河流位置对应的所有像素点的像素值记作T1,其余像素点的像素值保持不变,如式(16),将融合后的灾后SAR图像记作Fu(见图16)

    图  16  融合后的灾后SAR图像Fu
    Figure  16.  Fusion post-disaster SAR image Fu
    (Fu)ij={T1,(Eb)ij=1(Ea)ij,otherwise
    (16)

    其中,(Fu)ijFu, (Eb)ijEb, (Ea)ijEa

    光学图像的水体和非水体的区分可以通过计算归一化水指数(Normalised Difference Water Index, NDWI)[27]检测,但是SAR图像中不含多通道图像信息,不能直接计算其NDWI来区分水体和非水体。通过观察融合后SAR图像的高像素值区间Gm7(见图17),可以发现真实洪灾区域的像素点在其中几乎全部表现为零像素值。同时,根据RCC原理,真实洪灾区域在整个图像空间呈现较高的区块性的稀疏度。即代表洪灾区域的非零像素值区域是一个连通区域,其余零像素值区域均为非洪灾区域,所以作为洪灾区域检测的聚类模型应当具备较高的稀疏度。Gm7是一个反向模型,洪灾区域基本呈现为零像素区域,非洪灾区域基本呈现为非零像素的散点分布。

    图  17  聚类模型Gm7
    Figure  17.  Clustering model Gm7

    在目标检测中,文献[28]定义图像稀疏度为非零像素点个数K1,考虑到洪灾区域检测的目标地物覆盖范围较广,本文进一步引入文献[29]中图像稀疏占比的定义。即模型的稀疏占比为待测图像中非零像素点个数和所有像素点个数的比例,稀疏占比K2的定义式为公式(17)。聚类模型的稀疏度分析,以Gm7聚类模型为例,检测结果如图18

    图  18  聚类模型的稀疏度分析
    Figure  18.  Analysis of sparsity of clustering model
    K1=n1i=1n2j=1(Gm7)ij,K2=K1n1×n2
    (17)

    通过分析聚类模型Gm7的稀疏度,得到其稀疏占比为0.10905,可以确认该模型本身具有较高的稀疏度。并且该模型的非零像素点代表最高级像素值的分布,实际洪灾区域在该模型中基本处于零像素值区域。对聚类模型Gm7进行卷积运算后,根据卷积结果Gtk7得到每一个像素点周围非零像素点的总数,并由此判断该像素点是否为洪灾区域。

    Gtk7=Gm7Wk
    (18)

    其中,Wk表示大小为k×k的全1卷积核。

    H-FCM算法对聚类系数φ和卷积系数k的选择较为敏感。当选择的φ过大时,聚类结果包含较多的非洪灾区域;当选择的φ较小时,聚类结果舍弃了较多真实洪灾区域。相应地,当选择的k过大或者过小时,聚类结果都将包含较多的非洪灾区域。为使这两类参数选择达到最优检测性能,应当在选择最佳聚类模型的同时,依据模型的稀疏度来确定相应的聚类系数和卷积系数。考虑到SAR图中大量随机分布的相干斑噪声,为减弱噪声对检测性能的影响,φ,k的选择要比实际的稀疏占比更小,φ,kK2的关系定义为式(19)。

    φ=K22,k=φn1×n2
    (19)
    (Gr7)ij={1,(Gtk7)ij<φ(k)20,(Gtk7)ijφ(k)2
    (20)

    其中,(Gr7)ijGr7; (Gtk7)ijGtk7

    根据卷积结果Gtk7,用阈值法进行二分类,得到初步聚类结果Gr7,即洪灾区域的初步检测结果(见图19)。观察图19的初步聚类结果可以发现,反向聚类的方法虽然减少了洪灾区域相干斑噪声的影响,但是由于洪灾发生时地表坑洼处存在积水,初始聚类的结果还存在大量区块的虚警。在3.4节中,根据融合后SAR图像中的灾害前河流位置信息,对洪灾区域进行空间约束,进一步提升检测性能。

    图  19  初步聚类结果Gr7
    Figure  19.  Preliminary clustering results Gr7

    H-FCM算法在对洪灾区域的空间约束上提出了一种参数自适应的空间约束方法,有效地筛除疑似洪灾区域,并提升检测性能。由于洪灾区域的水体满足RCC原理,定义洪灾区域主要范围为RtRt表示初步聚类结果Gr7的4连通区域面积最大的连通区域。

    由于真实洪灾区域在整个图像空间存在较高的稀疏度,即非洪灾区域处于零像素值区域,洪灾区域呈现为非零像素值的连通区域。此时,等间距增大空间约束范围会增加计算复杂度。引入约束矩阵CSβp作为初步聚类结果的空间约束,当约束系数βp不断减少时,整体的空间约束范围按照e指数级不断增大。由于真实洪灾区域在全局上具备较高的稀疏度,约束系数βp按照等间隔减小的过程,对应约束范围呈e指数级的扩大。因此,式(21)的搜索规则符合洪灾区域的区块稀疏分布特点,既能保证获得适应待检测图像的约束系数,又不增加额外的运算量。

    约束系数βp的自适应选择方法,主要是依据两个因素:

    (1) 根据RCC原理,洪灾发生时,洪灾区域在遥感图像上满足4连通区域的连通性;

    (2) 远离洪灾前河流位置的低洼地区在洪灾发生时与洪灾区域相邻接的概率较低。

    βp=βp1e1100e,β0=1
    (21)

    其中,e为自然常数,p=1,2,···,100

    (CSβp)ij={1,exp(dijdmax)>βp0,exp(dijdmax)βp
    (22)

    其中,(CSβp)ijCSβp; dij表示Ea中第i行,第j列像素点与灾害前河流位置(即标记为T1的像素点)的最小曼哈顿距离,dmax表示所有dij元素中的最大值。

    本文采用曼哈顿距离来表示像素点之间的空间距离,可以在保证检测性能的基础上,降低算法的运算量。洪灾区域的空间约束式exp(dij/dmax)的取值范围为(0.367, 1](见图20)。并且当dij越小时,exp(dij/dmax)的值越接近1,当dij越大时,exp(dij/dmax)的值越接近0.367,确保距离灾害前河流位置越近的像素点被检测为洪灾区域的概率越大。

    图  20  洪灾区域空间约束曲线
    Figure  20.  Spatial constraint curve of flooded area

    将洪灾区域主要范围Rt加上空间约束CSβp,可以进一步提升洪灾区域检测结果的准确率。由式(23)可以获得约束系数为βp时的洪灾区域检测结果Outβp

    Outβp=RtCSβp
    (23)

    其中,“A·B”表示矩阵A与矩阵B的点乘运算,即矩阵A和矩阵B对应点相乘。

    约束系数βp自适应选择的目的是判定疑似灾区S是否为真正灾区,以获得洪灾区域的最终检测结果。所以在βp的搜索过程中,若对应的疑似灾区的面积增长比率远小于聚类模型本身的稀疏度时,则判定疑似灾区为未受灾区域

    S=OutβpOutβp1
    (24)

    其中,“–”表示矩阵对应点像素值作差

    ratiop=n1m=1n2q=1(S)mqn1m=1n2q=1(Outβp1)mq
    (25)
    βp={βa,ratiopK2βb,ratiop<0.01×K2
    (26)

    其中,“≈”近似的误差为±0.001; p=1,2,···,100

    FF=lj=1[Sjdtj]+Outβa
    (27)

    其中,Sj表示S中的第j个4连通区域,其背景大小为n1×n2

    dtj={1,α1<Ajlj<α20,otherwise
    (28)

    其中,AjSj的面积,α1,α2为判定系数

    lj=n1m=1n2q=1(CHj)mq
    (29)

    其中,lj表示重合边界长度;(CHj)mqCHj; CHjSjOutβa的重合边界,其宽度为2个像素单元。

    将得到的约束系数βa, βb带入到式(24),获得待定的疑似灾区S。根据RCC原理,若疑似的4连通区域为真实洪灾区域,则其像素面积AjSjOutβa邻接的边长lj的比值应处于(α1,α2)区间内。当面积边长比过大时,说明疑似灾区与真实洪灾区域相邻接的边界较短,可判定为非灾区;当面积边长比过小时,说明疑似灾区面积较小,属于区块噪声的可能性比较大。最后,根据每一个疑似区域的判定结果dtj确定该区域是否为洪灾区域,并将判定结果与Outβa累加得到最终的洪灾区域检测结果FF

    (1) 正确检测概率。Righta表示所有像素点被正确检测为变化的概率与被正确检测为未变化的概率之和

    Righta=Nr_c+Nr_ucNt
    (30)

    其中,Nr_c表示正确检测为变化的像素点总数,Nr_uc表示正确检测为未变化的像素点总数,Nt表示像素点总数。

    (2) 遗漏检测概率。Missa表示所有实际发生变化的像素点没有被检测为发生变化的概率

    Missa=NmNc
    (31)

    其中,Nc表示实际发生变化的像素点总数;Nm表示实际变化的像素点未被检测为变化的像素点总数。

    (3) 错误检测概率。Wra表示所有像素点被错误检测的概率,即实际未发生变化的像素点被错误检测为发生变化的概率与实际发生变化的像素点被检测为未变化的概率之和

    Wra=Nf_uNun_c+Nf_cNc
    (32)

    其中,Nf_u表示实际发生变化的像素点被检测为未变化的总数;Nun_c表示实际未发生变化的像素点总数,Nf_c表示实际未发生变化像素点被检测为变化的总数。

    (4) 卡帕系数

    Kappa=pope1pe
    (33)

    其中,po表示总体分类精度。

    po=Nr_c+Nr_ucNt
    (34)
    pe=Nd_cNc+Nd_ucNr_ucNtNt
    (35)

    Nd_c表示检测为发生变化的像素点总数,Nd_uc表示检测为未变化的像素点总数。

    (5) 计算时间。以计算时间衡量算法的时间复杂度,计算时间定义为Time,其单位为“s”(秒)。

    (6) 迭代次数。以运算的迭代次数体现算法的运算复杂度。迭代次数定义为Iter,其单位为“次”,表示对应算法获得聚类中心的迭代次数。

    本文的实验数据是在确定目标灾害区域地理坐标的基础上,获得洪灾前遥感图像和洪灾后遥感图像的历史数据。实验数据的预处理过程主要分为4步。首先,利用ENVI5.3软件的Orthorectification模块对洪灾前后的遥感数据分别进行正射矫正;其次,利用ENVI5.3软件Registration模块中的Image to Image方式,以洪灾前的光学图像为基准,选择关键角点为控制点对同地区的其它遥感数据进行配准;再次,裁剪相同大小配准后遥感数据作为实验数据;最后,利用ENVI软件获得洪灾后光学图像的NDVI图像,通过人工标注的方法确定真实洪灾区域,作为各类算法检测结果性能指标评定的参照对象。

    H-FCM算法的参数选择,只需要确定判定系数α1, α2。卷积系数k、聚类系数φ可以通过分析聚类模型的稀疏度获得,约束系数βa, βb可以通过自适应选择获得。判定系数的范围是通过统计英国格洛斯特的历史洪灾数据的贝叶斯更新结果[30]获得的最优区间,在H-FCM算法中被设置为固定的区间。并且该区间定义在4.3节中实验1和实验2的结果中均获得了较好的检测性能。

    为了验证H-FCM算法在洪灾检测领域的普适性,本文通过收集案例中的河流历史洪灾信息,利用历史实测数据来确定不同河段发生洪灾的空间约束系数。基于该河流的历史灾害数据进行贝叶斯更新[30],可以发现某一河段附近地区受洪灾影响的范围与该河段的河流宽度、河床高度、河道转角系数密切相关。参考英国格洛斯特的历史灾害数据[31],可以得到洪灾区域的约束系数βa, βb与河流平均宽度lm的近似对应关系

    βa=exp(1002lmdmax)
    (36)
    βb=exp(2003lmdmax)
    (37)

    其中,lm通过计算图15的平均河宽得到。

    计算实际数据得到dmax=1338, lm=4.5,其单位均为像素单元。依据真实灾害数据统计的近似对应关系得到的βa=0.8997, βb=0.7719。对比自适应选择的约束系数,βa的选择结果与洪灾的统计历史数据较为吻合,但是βb的选择结果与洪灾的统计历史数据有一部分偏差。这是由于该流域的历史洪灾规模不同,本例所体现的洪灾区域要大于历史平均水平,所以本例的βb的自适应选择结果要小于根据历史统计得出的结果,对应实际更大的洪灾区域。

    为充分论证H-FCM算法对洪灾前光学图像和洪灾后SAR图像的变化检测性能,在对比实验中加入了图像分割算法—分水岭算法(Watershed Algorithm, WA)[32],边缘检测算法(Snake)[33], H-Kmeans算法等3种方法作为对比实验方法。此外,将1999年英国格洛斯特洪灾数据、2019年中国南昌附近洪灾的洪灾数据等两个场景作为实验数据。

    4.3.1   实验1

    以1999年英国格洛斯特洪灾前的近红外波段图像和洪灾后的SAR图像[10]为例,检验H-FCM算法的有效性。将融合后的SAR图像(见图16)记作Fu,其大小为2359×1318像素,采用的聚类模型为Gm7。设定判定系数α1=120, α2=160,计算获得其卷积系数k=85,聚类系数φ=0.054。通过自适应选择得到的约束系数βa=0.9199, βb=0.7199,最终的检测结果为图21

    图  21  最终检测结果FF
    Figure  21.  Final experimental result FF

    首先,将实验1中洪灾后的NDVI图[10] (见图22)进行人工标注获得洪灾区域的真值图(见图23)。WA算法的洪灾区域检测方法是根据整体像素值分布(见图24),将所有像素点进行二值化,然后求出零值与最近非零值的距离,画出分水岭,进而对图像进行二分类。Snake算法获得Ea的边缘轮廓(见图25),按照式(18)和式(20)对边缘轮廓检测结果进行卷积运算得到最终洪灾区域检测结果。H-Kmeans算法是在Kmeans聚类(八分类)的基础上,进行分级聚类,从而得到洪灾区域的检测结果。

    图  22  洪灾后的NDVI图像[10]
    Figure  22.  NDVI images after the flood[10]
    图  23  人工标注的真值图
    Figure  23.  Manually labeled truth map
    图  24  像素值分布三维图
    Figure  24.  3-D map of the distribution of pixel values
    图  25  Ea的边缘轮廓
    Figure  25.  Edge profile of Ea

    其次,将H-FCM算法、WA算法、Snake算法以及H-Kmeans算法获得的洪灾区域检测结果与真值图(图23)比较,得到4种算法的检测结果与真实灾区的比较结果(见图26)。在图26中,白色部分表示实际洪灾区域被检测为洪灾区域;绿色部分表示实际洪灾区域被检测为非洪灾区域;红色区域表示实际非洪灾区域被检测为洪灾区域;黑色部分表示实际非洪灾区域被检测为非洪灾区域。本文后续实验的检测比较结果与图26的表示方法一致。

    图  26  4种算法的洪灾区域检测结果
    Figure  26.  Detection results of flooded area based on four algorithms

    根据4.1节中对6项评价指标的规定,按照4种算法分别用matlab进行计算分析,得到每一种算法检测结果对应的性能指标(见表3)。观察表3的各项评价指标可以发现H-FCM算法有最高的正确检测概率、最低的错误检测概率以及最高的Kappa系数。H-Kmean算法虽然有着与H-FCM算法接近一致的检测性能,但是其计算成本较高。综合考虑洪灾区域检测的高准确检测率和高时效性的要求,H-FCM算法具有相对最优的检测性能。

    表  3  融合后SAR图像洪灾区域的检测结果
    Table  3.  Detection results of flooded area in fusioned SAR image
    方法H-FCMWASnakeH-Kmeans
    Righta0.89460.43460.67570.8837
    Missa0.15230.06800.16070.1519
    Wra0.25130.70010.50700.2634
    Kappa0.60920.18400.29640.5837
    Time(s)51.982.722.5373.27
    Iter(次)251125
    下载: 导出CSV 
    | 显示表格
    4.3.2   实验2

    实验2的数据来源于欧洲航天局的开源数据,分别为哨兵1号和哨兵2号卫星在2019年6月对中国南昌附近的观测数据(见图28)。按照H-FCM算法,首先,将图28(a)中提取的的洪灾前河流空间位置与图28(d)中洪灾后的SAR图像融合;其次,对融合图像进行分级聚类;最后,对聚类结果加以空间约束并获得检测结果。根据图28(c)的洪灾后的光学数据,人工标注真实洪灾区域范围(见图27)。比较4种检测方法对融合后SAR图像的洪灾区域检测结果(见图29)和性能指标(见表4),可以得出H-FCM算法具有相对最优的检测性能。

    图  27  人工标注的洪灾区域
    Figure  27.  Manually labeled flooded areas
    表  4  融合后SAR图像洪灾区域的检测结果
    Table  4.  Detection results of flooded area in the fusioned SAR image
    方法H-FCMWASnakeH-Kmeans
    Righta0.94140.62300.64260.9391
    Missa0.25200.05550.06080.2434
    Wra0.28860.46900.45180.2836
    Kappa0.69110.26290.27490.6848
    Time(s)30.592.823.1479.13
    Iter(次)251125
    下载: 导出CSV 
    | 显示表格
    图  28  洪灾前后哨兵1,哨兵2数据
    Figure  28.  Sentinel 1, 2 data before and after flood
    图  29  4种算法的洪灾区域检测结果
    Figure  29.  Detection results of flooded area based on four algorithms

    实验2的参数设定如下,将融合后的SAR图像记作NCa,其大小为1025×1025像素,采用的聚类模型为Gm7;设定判定系数α1=120, α2=160;卷积系数k,聚类系数φ和约束系数βa, βb均基于NCa由式(19)和式(26)获得。

    Refined Lee filter(RL)[34]在滤除相干斑噪声的同时能较好地保持SAR图像的边缘特征。为分析图像降噪处理对上述算法检测性能的影响,首先对实验1中洪灾后的SAR图像进行RL滤波降噪,再利用上述4种方法对降噪处理后的图片进行洪灾区域检测。定义新的检测方法分别为RL_H-FCM算法、RL_WA算法、RL_Snake算法、RL_H-Kmeans算法,其检测结果如图30所示,检测性能评价指标如表5所示。

    表  5  RL滤波后SAR图像的洪灾区域检测结果
    Table  5.  Detection results of flooded area in the SAR image after RL filtering
    方法RL_H-FCMRL_WARL_SnakeRL_H-Kmeans
    Righta0.89960.43910.76150.8899
    Missa0.17500.06600.25980.1351
    Wra0.26530.69330.49540.2418
    Kappa0.61470.18630.34600.6032
    Time(s)646.04593.72593.59745.62
    Iter(次)251125
    下载: 导出CSV 
    | 显示表格
    图  30  RL滤波后4种算法的洪灾区域检测结果
    Figure  30.  Detection results of flooded area based on four algorithms after RL filtering

    比较表3表5中上述4种算法对应的性能指标,可以发现RL_Snake算法的检测性能提升最为明显。这是由于Snake算法是先获得图像的边缘特征,再对边缘轮廓的内外进行区分。通过RL滤波可以减少相干斑噪声的影响,使得图像整体像素值的均值和方差更稳定,有利于提取图像的边缘特征。但是通过RL滤波对H-FCM算法、WA算法和H-Kmeans算法的性能提升有限,反而极大地增加了算法的计算复杂度。总体而言,利用H-FCM算法直接对洪灾前的光学图像和洪灾后的SAR图像进行洪灾区域检测时,其检测性能与降噪处理后的图像检测性能几乎一致,而直接检测的运算时间大幅度减少。所以应用H-FCM算法直接处理配准后的光学和SAR图像,既可以缩短洪灾区域的检测时间又可以保证相对较高的检测准确率。

    本文提出的H-FCM算法综合考虑洪灾区域检测的高准确率和高时效性的需求,既提升了SAR图像的洪灾区域检测精度又缩短了检测时间。H-FCM算法区别于传统的无监督聚类算法,没有直接对待处理的图像进行无监督聚类,而是提出了一种反向聚类的方法。首先,H-FCM算法将洪灾前光学图像的河流位置检测结果与洪灾后的SAR图像融合;其次,在融合图像的最高级像素值区间对零像素值点进行分级聚类;最后,H-FCM算法在对融合图像的洪灾区域进行分级聚类的基础上,将洪灾前的河流位置作为空间约束,进一步筛除疑似洪灾区域以提升检测性能。

    H-FCM算法将洪灾前的光学图像和洪灾后的SAR图像融合后进行洪灾区域检测,实现了全天侯、高准确率、高时效的洪灾区域检测要求。虽然H-FCM算法在一定程度上解决了SAR图像相干斑噪声给洪灾区域检测造成的问题,但是还有部分灾区的细节没有精确检测到,这有待在未来的工作中解决。

  • 图  1  美国Sandia实验室公布的视频SAR结果

    Figure  1.  SAR video released by Sandia laboratory

    图  2  SAR图像及距离多普勒谱示意图

    Figure  2.  Illustrations of SAR image and RD spectrum

    图  3  基于快速区域卷积神经网络的双域联合检测流程图

    Figure  3.  Flow chart of joint detection algorithm based on Dual Faster R-CNN

    图  4  跨域航迹联合的动目标检测算法流程图

    Figure  4.  Flow chart of joint detection algorithm based on JTA

    图  5  经典Faster R-CNN与Dual Faster R-CNN检测结果对比

    Figure  5.  Comparison results of classical Faster R-CNN and Dual Faster R-CNN

    图  6  图像域帧间关联与JTA算法检测结果对比

    Figure  6.  Comparison results of data association in image domain and JTA algorithm

    表  1  运动目标阴影及其能量所在帧

    Table  1.   Frame number of target shadow and energy

    动目标阴影所在帧能量所在帧
    T1[1, 14][1, 14]
    T2[1, 64][1, 64]
    T3[1, 24][1, 38]
    T4[9, 28][8, 42]
    T5[26, 64][20, 64]
    下载: 导出CSV

    表  2  动目标检测性能统计结果

    Table  2.   Statistical results of moving target detection in SAR imagery

    方法虚警漏警
    图像域帧间关联7563
    Faster R-CNN716
    Dual Faster R-CNN017
    联合域航迹关联(JTA)115
    下载: 导出CSV

    表  3  处理耗时对比结果

    Table  3.   Comparison of processing times

    方法单帧平均耗时(s)
    图像域帧间关联1.58
    Faster R-CNN1.55
    Dual Faster R-CNN2.89
    联合域航迹关联(JTA)2.21
    下载: 导出CSV
  • [1] DAMINI A, BALAJI B, PARRY C, et al. A videoSAR mode for the x-band wideband experimental airborne radar[C]. SPIE 7699, Algorithms for Synthetic Aperture Radar Imagery XVII, Orlando, USA, 2010: 76990E.
    [2] MILLER J, BISHOP E, and DOERRY A. An application of backprojection for video SAR image formation exploiting a subaperature circular shift register[C]. SPIE 8746, Algorithms for Synthetic Aperture Radar Imagery XX, Baltimore, USA, 2013: 874609.
    [3] WALLACE H B. Development of a video SAR for FMV through clouds[C]. SPIE 9479, Architecture/Open Business Model Net-Centric Systems and Defense Transformation 2015, Baltimore, USA, 2015: 94790L.
    [4] CHEN H C and MCGILLEM C D. Target motion compensation by spectrum shifting in synthetic aperture radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 1992, 28(3): 895–901. doi: 10.1109/7.256313
    [5] BARBAROSSA S and FARINA A. Detection and imaging of moving objects with synthetic aperture radar. Part 2: Joint time-frequency analysis by Wigner-Ville distribution[J]. Journals & Magazines, 1992, 139(1): 89–97. doi: 10.1049/ip-f-2.1992.0011
    [6] HUANG Penghui, LIAO Guisheng, YANG Zhiwei, et al. A fast SAR imaging method for ground moving target using a second-order WVD transform[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(4): 1940–1956. doi: 10.1109/TGRS.2015.2490582
    [7] LIGHTSTONE L, FAUBERT D, and REMPEL G. Multiple phase centre DPCA for airborne radar[C]. 1991 IEEE National Radar Conference, Los Angeles, USA, 1991: 36–40.
    [8] ZHANG Yuhong. Along track interferometry synthetic aperture radar (ATI-SAR) techniques for ground moving target detection[R]. AFRL-SN-RS-TR-2005-410, 2006.
    [9] 丁金闪. 视频SAR成像与动目标阴影检测技术[J]. 雷达学报, 2020, 9(2): 321–334. doi: 10.12000/JR20018

    DING Jinshan. Focusing algorithms and moving target detection based on video SAR[J]. Journal of Radars, 2020, 9(2): 321–334. doi: 10.12000/JR20018
    [10] JAHANGIR M. Moving target detection for synthetic aperture radar via shadow detection[C]. 2007 IET International Conference on Radar Systems, Edinburgh, UK, 2007: 1–5.
    [11] WANG Hui, CHEN Zhansheng, and ZHENG Shichao. Preliminary research of low-RCS moving target detection based on Ka-band video SAR[J]. IEEE Geoscience and Remote Sensing Letters, 2017, 14(6): 811–815. doi: 10.1109/LGRS.2017.2679755
    [12] ZHANG Ying, MAO Xinhua, YAN He, et al. A novel approach to moving targets shadow detection in VideoSAR imagery sequence[C]. 2017 IEEE International Geoscience and Remote Sensing Symposium, Fort Worth, USA, 2017: 606–609.
    [13] TIAN Xiaoqing, LIU Jing, MALLICK M, et al. Simultaneous detection and tracking of moving-target shadows in ViSAR imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59(2): 1182–1199. doi: 10.1109/TGRS.2020.2998782
    [14] ZHAO Baojun, HAN Yuqi, WANG Hongshuo, et al. Robust shadow tracking for video SAR[J]. IEEE Geoscience and Remote Sensing Letters, 2021, 18(5): 821–825. doi: 10.1109/LGRS.2020.2988165
    [15] DING Jinshan, WEN Liwu, ZHONG Chao, et al. Video SAR moving target indication using deep neural network[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58(10): 7194–7204. doi: 10.1109/TGRS.2020.2980419
    [16] ZHANG Yun, YANG Shiyu, LI Hongbo, et al. Shadow tracking of moving target based on CNN for video SAR system[C]. 2018 IEEE International Geoscience and Remote Sensing Symposium, Valencia, Spain, 2018: 4399–4402.
    [17] ZHANG Hao and LIU Zhe. Moving target shadow detection based on deep learning in video SAR[C]. 2021 IEEE International Geoscience and Remote Sensing Symposium, Brussels, Belgium, 2021: 4155–4158.
    [18] WEN Liwu, DING Jinshan, and LOFFELD O. Video SAR moving target detection using dual faster R-CNN[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 2984–2994. doi: 10.1109/JSTARS.2021.3062176
    [19] REN Shaoqing, HE Kaiming, GIRSHICK R, et al. Faster R-CNN: Towards real-time object detection with region proposal networks[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2017, 39(6): 1137–1149. doi: 10.1109/TPAMI.2016.2577031
    [20] ZHONG Chao, DING Jinshan, and ZHANG Yuhong. Joint tracking of moving target in single-channel video SAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5212718. doi: 10.1109/TGRS.2021.3115491
    [21] BAR-SHALOM Y, FORTMANN T E, and CABLE P G. Tracking and data association[J]. The Journal of the Acoustical Society of America, 1990, 87(2): 918–919. doi: 10.1121/1.398863
    [22] REID D. An algorithm for tracking multiple targets[J]. IEEE Transactions on Automatic Control, 1979, 24(6): 843–854. doi: 10.1109/TAC.1979.1102177
    [23] ROHLING H. Radar CFAR thresholding in clutter and multiple target situations[J]. IEEE Transactions on Aerospace and Electronic Systems, 1983, AES-19(4): 608–621. doi: 10.1109/TAES.1983.309350
    [24] SINGER R A and SEA R G. New results in optimizing surveillance system tracking and data correlation performance in dense multitarget environments[J]. IEEE Transactions on Automatic Control, 1973, 18(6): 571–582. doi: 10.1109/TAC.1973.1100421
    [25] RAMACHANDRA K V and DIVISION R C. Optimum steady state position, velocity, and acceleration estimation using noisy sampled position data[J]. IEEE Transactions on Aerospace and Electronic Systems, 1987, AES-23(5): 705–708. doi: 10.1109/TAES.1987.310865
  • 期刊类型引用(0)

    其他类型引用(2)

  • 加载中
图(6) / 表(3)
计量
  • 文章访问数: 1655
  • HTML全文浏览量: 730
  • PDF下载量: 258
  • 被引次数: 2
出版历程
  • 收稿日期:  2022-03-02
  • 修回日期:  2022-04-29
  • 网络出版日期:  2022-05-24
  • 刊出日期:  2022-06-28

目录

/

返回文章
返回