利用高分三号SAR影像进行双侧变化检测

蔡宣宣，张永红，崔斌，2，3

(1.中国测绘科学研究院，北京 100039；2.武汉大学 测绘学院，武汉430079；3.城市空间信息工程北京市重点实验室，北京 100089)

0 引言

利用遥感影像进行地表覆盖和变化检测已经越来越受到广泛的关注和研究，其应用方向也在不断增加，如民用领域的城市扩张研究[1-3]、农业及森林的调查监测[4]，军事领域中的打击效果评估、舰船位置变化监测等[5]。与传统光学遥感不同，SAR卫星作为主动式微波成像卫星，成像不受天气和光照条件限制，也正因这一特性，使利用SAR卫星进行地表覆盖变化检测有着更加广阔的利用前景和能力[6]。高分三号作为我国首颗C波段合成孔径雷达卫星，于2016年8月发射完成，是我国高分对地观测系统的重要组成部分，拥有聚束、条带成像、扫描成像等5大类12种成像模式，成像模式数量为世界上同类卫星之最，最高分辨率可达1 m[7-8]，其具体卫星参数如表1所示。在海洋、减灾、水利和测绘等领域都有着重要的研究意义和应用前景。利用高分三号进行SAR变化检测，达到了利用国产卫星进行全天时，全天候进行变化检测的目的，极大改变了我国民用星载SAR图像依靠进口的状态，为国内各行业用户提供高质量、高精度对地观测数据[9]。对国产卫星应用和SAR变化检测研究都有着极其重要的意义。

表1 高分三号属性信息

单极化SAR影像变化检测主要分为3个步骤，其中第二和第三步为当前的主要研究方向[10]。分别为：①图像预处理。包括影像的配准、定标、裁剪、滤波等，为下步生成差异图提供2个时相影像。②差异图生成。这步主要是通过对2个时相影像的差异运算，反应2时相影像间对应像素的差异度，主要有LR、领域比(neighborhood based ratio，NR)[11]、均值比(mean-ratio,MR)[12]、似然比(likelihood ratio，LLI)[10]等。③差异图分割。其目的是通过分割算法将差异图中的变化图斑与非变化图斑进行分离，被广泛使用的算法包括大津法(OTSU)[13]、最小误差法(KI)[14]等。

广义高斯模型及KI阈值分割(GKIT)阈值分割算法在分割LR差异影像上有着较为良好的应用效果，但由于GKIT算法假设差异图中只存在变化和非变化两类，就导致其只能发现单侧变化，极大削弱了其变化检测的能力[15]。针对这一缺陷，有学者[16]提出了利用DGKIT的方法，其主要思想是假设差异图中存在的后向散射增强，减弱和非变化类像素均服从广义高斯分布，在满足KI准则的条件下，确定出分割三类像素的阈值。但由于差异图的随机性，2种变化类像素之间不管是差异程度还是数量上都可能存在较大差距，导致在拟合和求取准则函数最小值的过程中，无法稳定获取正确的分割阈值，可能会导致在直方图中一侧确定出2个阈值而另一侧被忽略的情况，极大削弱了检测结果的可靠性。

本文正是针对进行双阈值变化检测这一目标，提出了在传统针对LR差异图进行单侧GKIT分割变化检测的基础上，改进为利用互为相反数的2幅LR差异影像，分别进行GKIT分割，得到双侧变化的初始分割阈值，这样可以稳定获取2个初始的双侧阈值，然后利用马尔科夫随机场(MRF)分割分别进行分割精化[17]，最后再通过差值阈值去除上步分割结果中伪变化区域，融合后得到最终变化检测结果。这样就避免了出现在直方图一侧的阈值被忽略的情况。文章首次以高分三号影像作为变化检测实验对象，对文中算法进行实验并验证了精度，最后简要分析了研究区发生变化的原因，证明了高分三号在变化检测尤其是舰船及近海区域进行变化检测的能力。

1 算法介绍

文中采用算法流程图如图1所示。

图1 高分三号数据进行变化检测算法流程图

1.1 图像预处理

进行变化检测之前，需要将2个时相原始单视复数据(single looking complex，SLC)影像经过配准、裁剪、强度提取、定标、多视处理等流程，之后为了降低SAR影像固有斑点噪声的影响，对影像进行滤波，用于生成差异图的强度影像[18]。

1.2 差异图生成

假设上步中用于生成差异图的2时相影像分别为I1和I2，由于SAR图像自身的噪声特性，应用LR不但能将乘性噪声转化为加性噪声，而且能较好地反映出2时相影像间的差异程度，所以本文中采用LR的方法。但与之前所不同，本文采用生成双向差异图的方法生成2幅对应像素互为相反数的差异影像Ds和Db：

(1)

和

(2)

假设Db中，图像中灰度值较大的像素代表后向散射可能发生增强的像素，且灰度值越大发生后向散射增强的程度就越大，而灰度较小的像素代表后向散射可能减弱的像素，且灰度值越小发生后向散射减弱的程度就越小；相对地，在Ds中，灰度较大的像素代表可能发生后向散射减弱的像素而灰度较小的像素代表可能后向散射增强的像素。

1.3 差异图分割

本文采用以2次单侧GKIT阈值作为初始分割阈值，再利用MRF分割进行精化的分割方法，并通过约束条件去除了对数比值较大但差值较小的伪变化区域，得到最终分割结果。

1)GKIT。假设在差异影像中，未发生变化的区域要远大于发生变化了的区域。GKIT的主要思想是假设广义高斯分布分布更适合描述LR差异图的概率分布，且差异影像中的变化和非变化像素均服从广义高斯分布。如公式(3)所示。

i∈(u，c)

(3)

式中：p(x|Ci)表示差异图像x中某一类Ci的概率密度函数；ai和bi是常数，可以通过与均值参数mi、方差参数σi和形状参数βi的关系求出，而均值、方差和形状参数都可以根据拟合的概率密度函数中求出相关参数[15]。

(4)

实际上由于假设中符合广义高斯分布的差异图直方图本身具有对称性，且尖峰拖尾形态中尖峰部分对应差异图中广大未变化区域的灰度值，而两侧拖尾形态对应的分别是差异图中后向散射系数增强与减弱区的灰度值范围。所以其在检测双侧变化(后向散射增强与减弱2种变化)时有很大缺陷。其拟合结果都将只是除拟合出最高尖峰外，准则函数较小的那一侧的阈值，而另一侧的变化阈值将被忽略。

所以具有近似对称性的广义高斯分布直方图在以下两种情况会有较好的效果，一是利用双阈值选取方法同时将差异图中的后向散射增强和减弱阈值确定，如DGKIT算法；二是差异影像中发生变化的变化类单一，均为单一后向散射增强或减弱变化。

本文在假设可拟合出阈值一侧的变化像素的均值大于非变化像素均值(可优先拟合出直方图尖峰右侧的阈值)的前提下，尝试通过对双向差异图分割来解决这一问题。对差异图Db应用GKIT得到阈值Td(>0)确定后向散射增强区域的初始分割结果；由于此时无法获得后向散射减弱区域的阈值，故只能大致获取初始阈值之后进行精化处理。本文就以对Ds应用GKIT获取的阈值Ts(<0)作为分割后向散射减弱区域的阈值，选取的阈值Td和Ts均为分割差异影像Db，最终得到分割后向散射增强区域的初始影像Gb和后向散射减弱区域的初始影像Gs。

2)MRF分割。由于GKIT是通过全局直方图选取阈值的分割结果，而全局阈值的缺点在于没有利用到局部领域信息，选取结果会受到局部噪声影响，且在无法拟合的尖峰一侧阈值选取结果不精确。故文章采用将上步GKIT分割结果作为初始分割，再利用MRF分割对2幅初始分割影像进行精化，一方面有效利用了领域信息，抑制了图像的噪声，减小孤立的变化像元与孤立非变化像元对最终分割结果的影响，另一方面也能精化无法拟合单侧阈值的一侧阈值选取不精确的问题。

MRF模型与Gibbs分布存在等效关系，其能较好描述图像空间的局部领域关系，使其在图像分割领域得到了广泛应用[17，19]。其目的是通过优化算法对图像进行递归求解系统能量函数最小值及涉及到的统计参数，获取对应的标号场，得到最终精化后的分割结果图Is和Ib。

1.4 分割后处理

文中实验数据采用的是未经拉伸的原始数据，灰度跨度较大，所以LR方法不可避免会遇到像素比值结果虚高的情况，比如研究区中出现的2组像素，第一组出现在陆地上且实际发生变化的2个对应像素灰度值为20和46.4，其对数比值为0.98，差值为26.4；另一组出现在海面上实际没有发生变化的2个对应像素分别为0.5和1.09，其对数比值同样为0.98，而差值只有0.59，这也从一个侧面反应了应用经典LR方法生成差异图的一个缺点[20]。为解决这一问题，在得到分割结果之后，将前后时相影像差值的绝对值影像S中阈值小于t，但分割结果是变化的像素分别从变化结果Fs和Fb上去除，只保留差值大于阈值t的部分，得到变化结果Fs和Fb，融合后得到最终的变化结果图F。

2 数据准备

2.1 前后时相影像数据

试验区原始数据基本信息如表2所示，本文中格式均为方位向在前，距离向在后，且未经过影像地理编码。研究区原始影像如图2所示。

表2 高分三号试验区数据基本信息

图2 宁波地区前后时相原始影像

2.2 实验区域概况

如图3所示，研究区为浙江宁波市区及周边海域。宁波位于浙江省东部，长江三角洲南翼，为浙江省经济中心，经济发展十分迅速。研究区的主要地物类型有城区、海域、滩涂，根据2016年数据，附近舟山港稳居全球第一大港口，海上船舶往来密集，可以用来测试高分三号的船舶检测能力。

图3 研究区范围

2.3 潜在可能发生真实地表变化

SAR影像变化可分为后向散射增强变化和减弱变化，如果根据地物后向散射强度将地物分为强、中、弱三类，则水面和裸地可认为是后向散射较弱和正常区域，船舶和人工建筑物则具有较强的后向散射强度。故反应在2时相SAR强度影像上有后向散射增强与减弱的变化之分。影像获取时间分别为2016-11-14和2017-03-10，期间由于植被变化和船舶移动等影响，引起影像上后向散射强度变化。

3 实验与分析

3.1 数据处理

1)预处理。预处理分为配准、裁剪、强度提取、多视处理和滤波等过程。多视处理的视数比为5∶4，裁剪后影像大小为4 071像素×3 753像素，分辨率为14 m左右。生成差异图时对影像进行了3×3增强LEE滤波处理，以减少斑点噪声对图像影响。滤波后结果如图4所示。

图4 滤波后前后时相强度影像

2)差异图生成。与传统的SAR变化检测中只生成一副差异图不同，本文不再以后时相除以前时相影像后取对数来确定差异影像，而是生成2幅互为相反数的差异影像Db和Ds，如图5所示。

图5 对数比差异图像

3)差异图分割。对经过中值滤波的差异图进行GKIT分割，假设Db中后向散射增强区域均值大于非变化区域均值，对Db进行分割得到Tb(=1.002 9)分割后向散射增强区域阈值，Db中大于Tb的像素是后向散射增强像素。对Ds进行分割得到Ts(=-1.097 9)分割后向散射减弱区域，差异图Db中小于Ts的像素被认为是后向散射减弱像素。初始分割结果如图6所示。

图6 GKIT初始分割结果

如图7所示，GKIT在直方图中只拟合出了尖峰右侧的阈值，而尖峰左侧无法拟合出阈值。故在选取初始阈值后采用MRF分割进行精化处理。

图7 差异图Ds各类型概率密度分布

经过MRF精化后的结果如图8所示。

图8 MRF精化分割结果

未经归一化原始影像生成LR差异图会带来比值较大而差值较小的伪变化区域，如图8(a)中大片白色海域变化区域，考虑通过设定差值差异阈值(本文t=5)，将伪变化区域去除，并将上步中后向散射增强变化结果和减弱变化结果进行融合，得到最终变化结果图。其中后向散射增强结果和减弱结果如图9所示，整体变化结果如图10。后向散射增强像素、非变化像素、减弱像素分别为255、128、0。

图9 伪变化区域去除后结果

图10 最终变化结果

3.2 典型变化与非变化图斑分析

图11是图10中区域1的2时相信息。变化主要原因是不同获取时间舰船位置及沿海区域发生变化。沿海地区发生了散射强度增强，主要原因可能是由于海产品养殖增加了反射地物和沿海含水量的变化；城镇地区发生此变化的原因主要是城市扩张使裸地变成城市建筑区。由于没有同时获取的光学影像，文中没有给出2时相对应的光学影像。图12是图10区域2的2时相SAR影像，由于地表植被及土壤含水量不同引起地表变化。图13所示为邻近时间光学影像，可以看出地表覆盖情况发生变化情况。

图11 区域1 2时相影像

图12 区域2前后时相影像

图13 区域2 2时相邻近时间光学影像

为计算方便，精度评价时，将后向散射增强和减弱像素均标记为变化像素(灰度值255)。其他情况下，仍然分别显示为后向散射增强和减弱区域。通过计算错检率(false alarms，FA)和漏检率(missed alarms，MA)，总错误率(overall error，OE)及Kappa系数4项作为评价指标进行精度评价。分别使用本文方法和只进行了去除伪变化区域而没有利用MRF分割精化的GKIT方法。由于缺少真实地表覆盖变化图，本文选取2个区域作为样本代替真实地表变化图。区域1，18 391个非变化像素，900个变化像素；区域2，2 284个非变化像素，912个变化像素。人工勾绘的变化图斑如图14所示。

图14 人工勾绘的变化图斑

为了验证本文方法的实用性，与DGKIT算法的检测结果进行了对比，DGKIT算法的检测结果如图15所示，其中2个分割阈值分别为-1.985 7和-1.093 7，可见DGKIT只在单侧确定了2个阈值而忽略了另一侧的变化情况，检测结果也反应出了这一结果，图中大部分没有用发生变化的区域都被错分为减弱区域。与真实地表变化相差较大，故没有对其进行精度评价，只是给出了其检测结果。

图15 DGKIT变化结果

区域1 2种方法分割结果和精度评价结果如表3和图16所示。可见经过MRF分割精化之后的结果，FA和MA分别减少了0.15%和0.16%，在本身局部已经具有较高精度的检测结果的基础上，考虑到邻域关系后，精度又有了部分提升，Kappa系数达到了0.889 1。而DGKIT算法的检测结果，如图16(c)中所示，与真实结果相差较大。而由本文算法检测结果图16(a)和图16(b)可见，利用适合的算法对高分三号影像进行处理，能在近海舰船检测(图16(a)中大部分黑白斑块)和沿海地区及海岸线检测方面产生较为良好的效果。

表3 区域1精度评价

图16 区域1变化检测结果

区域2分割结果和精度评价结果如表4和图17所示。在经过MRF精化分割结果后，尽管MA有小幅下降，但OE仍然有了1.15%的提升。最终的Kappa系数也达到了0.851 1。但DGKIT检测结果，检测结果全部为错误的黑色减弱像素，与真实地表变化相差较大。

表4 区域2精度评价

图17 区域2变化检测结果

3.3 变化结果分析

通过算法检测，提取了整个影像覆盖地区的变化情况，变化像素约占影像对应研究区大小的1.71%。从图中可以看出，变化主要分布在2个主要区域，一是集中在市内的鄞州区，二是北仑区和镇海区沿海。鄞州区由于区域内存在较多农业用地，相隔4个月后，土地植被种植情况及土壤含水量都发生了不同程度变化，但由城市建设引起的包含人工地物变化的情况则发生较少。相比之下，北仑、镇海区的沿海地区发生变化的主要原因是因为沿海滩涂，包括海带紫菜等近海养殖业及海面船舶位置变化引起的后向散射变化是造成这一区域变化的主要原因。

4 结束语

本文针对在SAR变化检测中应用DGKIT算法在增强和减弱像素变化存在较大差异的情况下可能会产生单侧出现2个阈值的情况，提出了利用双向差异图确定双侧初始阈值，并利用MRF分割精化，再去除结果中存在的伪变化区域，最终生成双侧变化检测结果的方法，确定了2时相影像中发生后向散射增强和减弱的变化区域。

通过对高分三号浙江宁波地区2016年11月14日和2017年3月10日2景的实验检测结果的分析，一方面证明了本文算法是一种有效的SAR变化检测方法，另一方面，首次利用高分三号SAR影像作为变化检测实验对象也显示出国产高分三号SAR影像在变化检测尤其是沿海滩涂及近海海产养殖、海岸线变化和舰船检测上有着较好的检测效果。