机场作为国家航天重要的基础设施，承载着飞机起降、维护等关键任务，在军事和民用领域具有重要的战略意义，对机场的自动检测也广泛应用于军事打击、灾害救援、飞机导航等实际问题。跑道作为识别机场的一个显著特征，可以通过检测跑道来实现对机场区域的检测。

目前机场跑道检测问题主要集中于光学遥感图像，已有的检测方法大致分为两步：第一步是跑道特征提取，常用的特征有跑道的灰度特性(艾淑芳 等，2017)，几何特征(Ebubekir和Hasan，2017)，直线段特征(Budak 等，2016；韩晓青和杜红悦，2016；艾淑芳 等，2017)，纹理特征(Aytekin 等，2013)，尺度不变特征变换SIFT(Scale-Invariant Feature Transform)(Budak 等，2016)等，其中直线段特征应用尤为广泛，但也会受到道路，河流，桥梁等同样具有直线段特征的地物的干扰，而纹理特征和SIFT特征虽然可能会得到更精确的检测结果，但同时也更加的耗时。第二步是选择合适的分类器如SVM(Budak 等，2016)，Adaboost(Aytekin 等，2013)进行最终的检测。基于光学图像的机场跑道检测主要问题是图像受到光照及天气的影响，不能保证全天时、全天候实时在线检测。

与光学遥感相比，极化合成孔径雷达PolSAR(Polarimetric Synthetic Aperture Radar)是一种先进的主动式微波成像系统，具有全天时、全天候、强透射性、远距离成像等特点。相对于传统单极化SAR，全极化SAR能够获取多种极化方式的数据，能够获取更多的地物极化信息，可以更好地表征地物细节和特征(Lee和Pottier，2009)。目前，全极化SAR图像越来越多地应用于地物目标检测、分类和识别等领域(Uhlmann和Kiranyaz，2014；张英海 等，2016；胡娇静 等，2018；车美琴 等，2016)。

在PolSAR的目标检测研究领域，目前大多集中在PolSAR舰船目标检测(Paes 等，2010；张嘉峰 等，2018；田巳睿 等，2007)，针对机场跑道检测的相关研究相对较少。已有研究主要针对单极化SAR功率图像，同样利用跑道的直线特征，通过边缘检测，霍夫变换，直线段检测器LSD(Line Segment Detector)等直线检测算法检测直线段，并结合先验知识(如机场跑道区域存在平行线对，且长宽必须符合机场跑道的建设标准等)筛选出符合条件的直线段，进而确定机场跑道区域所在的位置，再结合其他特征作进一步辨识(Liu 等，2018；乔寅骐 等，2016；贾承丽 等，2007；何勇 等，2004)；但是SAR功率图像中的边缘点容易被相干斑噪声淹没，导致边缘提取困难，直接影响检测结果，且图像场景中的道路、桥梁及河流容易影响跑道的检测结果。近几年提出的多极化SAR跑道检测算法中，往往采取如下流程：先利用目标极化特性并结合图像分类的方法得到感兴趣区域ROI(Region Of Interest)，然后再根据跑道的直线特征和其他特性对ROI进行筛选，最终得到真实的跑道区域(韩萍 等，2016；晋瑞锦 等，2014)。这类方法较好地解决了复杂大场景中河流、道路及桥梁对跑道检测的干扰，降低了虚警率，但是仍然存在以下问题：(1)对SAR图像分类时，通常需要设置分类的类别数，而图像中包含地物类别数目往往预先未知。已有的方法是根据经验人为设置(韩萍 等，2016)，或者根据某些特征硬性划分为固定的类别数(晋瑞锦 等, 2014；韩萍 等，2016)，不能自适应地确定。由于不同场景的SAR图像地物类别个数差别较大，导致算法推广性不强。(2)算法大多是基于单个像素进行处理，在不使用去噪算法对图像进行相干斑抑制的情况下，一些较强的噪声点依然会对分类性能造成干扰；并且，当SAR图像尺寸较大时，算法计算量大，处理效率低。

针对以上问题，本文提出一种基于自适应无监督分类的极化SAR图像机场跑道区域自动检测算法。通过超像素方法进行去噪并降低了原始图像维数，然后，采用无监督聚类方法自适应地估计类别。为进一步提高准确率，还利用了地物目标多数化特性作为分类特征。

VAT-DBE方法通过差异图RDI(Reordered Dissimilarity Image)显式反映数据间的差异度，其基本思想是通过构造数据间的差异矩阵，将差异矩阵沿主对角线排序并可视化，形成RDI，再利用投影信号的求导提取到潜在的聚类个数(Bezdek和Hathaway，2002)。具体步骤如下：

假设待测数据集为 ${{O}} = \left\{ {{o_1},{o_2},\cdots,{o_n},} \right\}$ ，n为数据个数。

(1)计算差异矩阵 ${{R}} = \left[ {{r_{pq}}} \right]$ ，式中， ${r_{pq}} = \left\| {{o_p} - {o_q}} \right\|$ ， $1 \leqslant p,q \leqslant n$ ， $\left\| \bullet \right\|$ 为向量的2范数。

(2)设 ${{K}} = \left\{ {1,2,\cdots,n} \right\}$ 为数据的索引集， ${{I}},{{J}}$ 为 ${{K}}$ 的子集，初始化 ${{I}} = {{J}} = \text{Ø} $ ；定义数组 ${{P}} = $ $ \left({0,0,\cdots,0} \right)$ ，用来存储排序后的数据索引。

(3)选择 $\left({i,j} \right) = \mathop {\arg \max }\limits_{p,q \in {{K}}} \left\{ {{r_{pq}}} \right\}$ ，表示所有 ${r_{pq}}$ 取得最大值的有序索引对 $\left({i,j} \right)$ 的集合；令 ${{P}}(1) = i$ ， ${{I}} = \left\{ i \right\}$ ， ${{J}} = {{K}} - \left\{ i \right\}$ 。

(4)对于任意 $t = 2,\cdots,n$ ，选择 $\left({i,j} \right) = \mathop {\arg \max }\limits_{p \in {{I}},q \in {{J}}} \left\{ {{r_{pq}}} \right\}$ ，令 ${{P}}(t) = j$ ， ${{I}} = {{I}} \cup \left\{ i \right\}$ ， ${{J}} = {{J}} - \left\{ j \right\}$ 。

(5)完成遍历后，令 ${\tilde{ R}}(i,j) = \left[ {{r_{P(i)P(j)}}} \right]$ ，即为重排序的差异矩阵；再对 ${\tilde{ R}}$ 进行归一化并转化为灰度图形式，得到RDI。

(6)用Otsu方法对RDI图分割形成二值图，形态学处理后通过二值图像距离变换将二值图转换为灰度图(距离度量为欧氏距离)。

(7)对灰度图所有像素值沿主对角线投影，即垂直主对角线的像素点分别相加，再作曲线拟合，得到差异度分布曲线。找到曲线极值点的个数即为数据集的估计类别数。

谱聚类算法是一种建立在谱图理论上的算法，与传统聚类方法相比具有能在任意形状的数据空间聚类并收敛于全局最优的特点(Ng 等，2001)。实际中常采用基于Nyström算法(Fowlkes 等，2004)的谱聚类以减小所花费的时间和空间代价。具体步骤如下：

假设待测数据集为 ${{O}} = \left\{ {{o_1},{o_2},\cdots,{o_n},} \right\}$ ，n为数据个数。

(1)计算相似矩阵 ${{{W}}_{n \times n}} = \left[ {{w_{ij}}} \right]$ ，式中， ${w_{ij}} = $ $ \exp \left({\dfrac{{ - {{\left( {d\left({{o_i} - {o_j}} \right)} \right)}^2}}}{{2{\sigma ^2}}}} \right)$ ，d为距离度量， $\sigma $ 为核函数尺度参数。

(2)计算拉普拉斯矩阵 ${{{L}}_{\rm{sym}}} = {{{D}}^{ - \frac{1}{2}}}({{D}} - {{W}}){{{D}}^{ - \frac{1}{2}}}$ ，其中 ${{D}}$ 为 ${{{W}}_{n \times n}}$ 的对角矩阵，其对角元素为 ${{D}}(i,j) = $ ${\sum\limits_{i = 1}^n w _{i,j}} $ 即 ${{{W}}_{n \times n}}$ 矩阵对应列的所有元素之和。

(3)计算拉普拉斯矩阵 ${{{L}}_{\rm{sym}}}$ 的特征值并取前k个最小特征值对应的特征向量，构造矩阵 ${{{V}}_{n \times k}} = $ $ \left\{ {{v_1},{v_2},\cdots,{v_k}} \right\}$ 。

(4)用k-means算法对 ${{{V}}_{n \times k}}$ 的行向量进行聚类，将n个k维的数据划分为k类，最终完成样本的分类。

本文算法主要包括预处理、分类、感兴趣区域筛选及后处理等部分组成。具体流程如图1所示。

图 1 算法流程图 Figure 1 Flow chart of the algorithm

预处理部分的主要目的是去噪及降维，为后续准确分类及快速处理奠定基础，一般情况下，PolSAR图像去噪主要是去除图像中的相干斑噪声，如典型的Lee滤波及近年来出现一些新方法(Mullissa 等，2017；行晓黎 等，2017；刘留 等，2017)。本文方法主要考虑图像中出现的虚影和伪影干扰去除问题，具体方法如下：

(1)特征值图计算。对原始PolSAR图像T矩阵进行特征分解，提取第3个特征值对应的图像，以去除图像中的虚影和伪影等(吴正升 等，2014；Wang 等，2010)便于后续分类使用。

如图2所示，图2(a)场景为美国Honolulu机场区域；图2(b)是其对应PauliRGB图，红框标出区域为受方位模糊干扰严重区域，可以看到在该区域由于方位模糊产生的竖直红线；在图2(c) ${\lambda _3}$ 特征值图；图2(d)为PauliRGB图和 ${\lambda _3}$ 特征值图受方位模糊干扰对应区域大图，可以看出前者有明显的红色竖直亮线，特征值图基本没有；因此，可以利用 ${\lambda _3}$ 特征值图生成超像素来减少干扰。

图 2 方位模糊在Pauli图和 ${\lambda _3}$ 特征值图中的表现 Figure 2 Representation of azimuth ambiguity in PauliRGB and eigenvalue image

(2)超像素生成。超像素指根据相邻像素间颜色、空间位置等特征将相似的像素合并形成的图像块。分割能够实现对图像相似像素点合并，降低图像冗余的同时保留图像有用细节，大大提高算法效率。本文对特征值图用SLIC算法(Achanta 等，2012)分割生成超像素图 ${I_0}$ ，根据边缘重合率通过实验设置超像素个数K=10000。

分类的目的是将图像中的不同地物有效的分开，为后续利用不同地物散射特性提取机场跑道区域打好基础，本部分包括特征提取、类别数估计及分类。

(1)特征提取及类别数估计。将超像素图 ${I_0}$ 中每个超像素包含的像素对应原始数据的平均极化信息作为该超像素的极化信息，超像素对应的极化相干矩阵特征 ${{\bar{ T}}_{si}}$ 表示为

式中，n为每个超像素中包含像素的个数， ${{{T}}_i}$ 为其中每个像素的极化相干矩阵。类似地，可以求得极化散射熵 ${\bar H_{si}}$ 、 ${\bar \alpha _{si}}$ 角以及功率 ${\bar P_{si}}$ 并将其作为特征。

将超像素 ${\bar P_{si}}$ 值作为类别数估计的输入，根据2.1节中VAT-DBE算法估计并提取出类别个数c。

(2)图像分类。利用上一步得到的特征分别构造 ${{{W}}_1}$ 和 ${{{W}}_2}$ 两个相似矩阵并构造 ${{W}}$ 作为最终谱聚类的相似矩阵，其中 ${{{W}}_1}$ 以 ${\bar T_{si}}$ 作为特征，距离度量采用Wishart距离； ${{{W}}_2}$ 以向量 ${{F}} = \left[ {{{\bar H}_{si}},{{\bar \alpha }_{si}},{{\bar P}_{si}}} \right]$ 作为特征，距离度量采用欧氏距离。

将最终谱聚类相似矩阵定义为

式中， $ \odot $ 表示矩阵对应元素相乘。得到相似矩阵 ${{W}}$ 后，结合类别数c及2.2节谱聚类算法步骤对图像进行分类，分类后的图像记为 ${I_1}$ 。

计算 ${I_1}$ 中各类别的平均功率，根据跑道的弱散射回波特性，提取平均功率最小的一类作为跑道疑似ROI区域，再利用形态学滤波去除孤立点，结合机场跑道尺寸特征、拓扑特征和平行线特征对疑似区域进一步筛选和辨识(韩萍 等，2014)，最终检测出真实跑道区域，具体步骤如图3所示。

图 3 机场跑道疑似ROI筛选过程 Figure 3 Suspected ROI screening process of airport runway

为验证本文所提算法的有效性，对多幅含有机场的PolSAR图像实测数据进行了实验，并与晋瑞锦等(2014)和Han等(2016)提出的方法(下文分别称为方法1和方法2进行了对比，鉴于篇幅有限在此给出其中两个代表性的实验结果进行分析。实验运行环境为：Dell2013 B6BVQZ1电脑，内存为4 G，CPU为Corei5-3470，3.2 GHz，实验所用软件为MATLAB 2016a。实验参数设置：尺寸特征阈值设置与韩萍等(2014)所提方法一致，超像素个数K=10000，形态学滤波阈值为M₀=459，欧拉数E₀=0，平行线对数阈值L₀=10。

实验一：数据采用美国UAVSAR系统在美国Honolulu地区上空采集的L波段全极化经过4视处理的数据。图像大小为600×900像素，图像距离分辨率为7.2 m，方向分辨率为4.9 m。该场景中除包括檀香山地区一个机场外，还包含有海洋、建筑物、草地、裸地等多种复杂地物。

图4(a)对应图2(c) ${\lambda _3}$ 特征值图生成的超像素图，超像素个数约为原始图像像素个数的百分之二，且形成超像素的大小均匀、地物类型边界保持良好，在有效降维的基础上，对后续处理不会造成影响。图4(b)表示自适应分类结果，其中类别数估计c=8，图中可以明显区分出浅绿色跑道区域，说明图像分类结果良好。图4(c)表示本文算法最终检测结果，对比图4(d)(方法1)，二者均能有效地检测出机场跑道，本文算法在细节上虽略有丢失但本文算法无需预先设置地物类别个数，减少了人为干预。图4(e)为方法2检测结果，从图中可以看出该方法未完整检出机场跑道且有虚警目标，是由于如图2(d)所示，极化SAR图像该区域受到方位模糊干扰，导致跑道断裂，在最终筛选提取时将跑道另一部分视为非跑道区域剔除。

图 4 Honolulu地区实验结果 Figure 4 Experiment results of Honolulu area

实验二：实验数据采用美国UAVSAR系统在美国Big Island地区上空采集的数据。图像大小为1151×1051像素，图像分辨率及实验参数设置同实验一，该场景中除包括科纳国际机场外，还包括建筑物、农田、森林、灌木、裸地等(图5)。

图 5 Big Island地区机场场景 Figure 5 The scenes of Big Island area airport

实验二分析过程同实验一。超像素个数大约为原始图像像素数的百分之一，c=6。最终检测结果如图6(d)所示。对比图6(e)和图6(f)，结果显示，方法1及方法2均能有效的检测出机场跑道，本文方法相对于方法2提取结果更加完整。但是，本文算法无需预先设置要分类的类别个数及迭代次数人为干预更少。

图 6 Big Island地区机场实验结果图 Figure 6 Results of Big Island area airport scenes

实验三：实验数据采用美国UAVSAR系统在美国San Andreas Fault地区上空，采集的数据大小为1051×1151像素，其他参数设置同实验一(图7)。

图 7 San Andreas Fault地区机场场景 Figure 7 The scenes of San Andreas Fault area airport

实验三的分析过程同实验一、二，对比图8(d)、图8(e)、图8(f)可知，本文算法在较少人为干预的情况下与两种对比方法都能够较完整地检测出机场跑道区域。

图 8 San Andreas Fault地区机场实验结果图 Figure 8 Results of Gulf San Andreas Fault area airport scenes

除了对给出的实验数据进行分析，在同一实验环境且未对算法进行优化的情况下，本算法耗时大约为方法1的1/10、方法2的1/5，如表1所示。综上可知，相对于方法1和方法2，本文所提方法能够自适应地确定地物类别个数，人为干预更少，算法运算速度更快、无虚警、漏警，所以实时性和适用性更好。

表 1(Table 1)

表 1 本章算法及对比算法时间统计 Table 1 Execution time of proposed method and comparative method /s 本章算法 方法1 方法2 数据一 178.455 2302.002 1020.125 数据二 380.410 4562.017 2599.766 数据三 399.606 7681.945 3316.572

表 1 本章算法及对比算法时间统计 Table 1 Execution time of proposed method and comparative method

在利用聚类算法提取疑似跑道区域时，提取结果往往受到初始类别数的影响，另外，基于像素点的检测易受到PolSAR图像中固有相干斑噪声的干扰，且当场景复杂时，算法计算量会非常巨大。本文提出的自适应无监督分类的机场跑道区域检测算法，具有以下特点和优点：

(1)利用SLIC算法生成超像素，一方面，生成的超像素会平滑强的噪声点，降低相干斑噪声的影响；另一方面，图像中超像素的个数将远小于原始像素点的个数，这样虽然牺牲了一定的分辨率，但是可以大幅降低算法的时间成本。

(2)利用VAT和DBE算法对待测图像中实际包含的地物类别数进行估计，提高了算法的自适应性。

(3)利用谱聚类的方法提取疑似跑道区域，与其他传统聚类算法相比，虽然谱聚类需具有更大的时间复杂度和空间复杂度，但是也同样具有更精确的分类性能，而且在将超像素作为处理单元的前提下，谱聚类稍高的计算成本对整体而言并没有太大的影响。

(4)本文算法充分利用了地物目标的极化特性作为分类特征，有助于提高分类的准确性。实验结果表明，本文算法虽然在丢失了一定分辨率的情况下，仍然得到了结构完整，轮廓清晰的检测结果，且大幅降低了算法耗时，提高了计算效率。

需要指出的是，虽然超像素可以有效降低时间成本，并对相干斑噪声有一定的抑制作用，但仍然会造成丢失细节的问题。在全部的检测结果中，小部分机场的细节略有损失，这可能设定与超像素个数有关，本文目前采取的做法是经验选取，后续或将结合图像本身信息研究超像素个数选取的相关方法。