全极化探地雷达多极化数据融合分析研究

薛策文 冯 晅 李晓天 梁文婧 周皓秋 王 颖

①(吉林大学地球探测科学与技术学院 长春 130026)

②(地球信息探测仪器教育部重点实验室(吉林大学) 长春 130026)

③(近地面探测技术重点实验室 无锡 214035)

④(吉林大学国家发展与安全研究所 长春 130012)

1 引言

探地雷达(Ground Penetrating Radar,GPR)是一种无损的利用电磁波传播来探测地下的地球物理方法，其已经在工程[1]、考古[2]、水文[3]、地雷探测[4]、极地及冰层探测[5]、星球探测[6]等领域有了广泛的应用。大部分传统的商业探地雷达只能测量共极化数据(VV)，而全极化探地雷达(Full-Polarimetric GPR,FP-GPR)不但可以获得共极化数据VV和HH，还可以获得交叉极化数据VH。目标体通过极化散射机制不同可以分为单次散射目标体、二次散射目标体和体散射目标体。目前FP-GPR数据处理大体可以分为两个方面，一种是极化分解方法，它是通过遥感极化分析手段结合FP-GPR数据进行处理。例如，Pauli分解技术结合偏移处理[7]，H-α分解结合偏移处理[8]，Freeman分解结合3D偏移成像[9]。另一种就是数据融合方法，将VV,HH,VH 3种图像进行图像融合，使融合后图像具备VV,HH,VH的所有特点，再用传统探地雷达处理方法进行处理，得到分辨率更高更好的地下目标体图像。

图像融合是综合来自不同传感器的同一目标的图像信息，通过对多幅图像信息的提取和综合，从而获得对同一目标的准确全面可靠的图像描述[10]。Daily等人[11]首次将雷达图像和陆地卫星图像融合到一起并应用于地质探测上。最近20年，图像融合已经应用于智能机器人、医学、遥感等方面[12—14]。例如，Ding等人[15]提出了一种基于高斯过程模型的多传感器数据融合方法来进行复杂的表面测量；Wan等人[16]提出一种鲁棒的主成分分析方法用于解决多焦点图像的融合问题。但是，在FP-GPR数据融合领域，目前常用的加权平均融合方法，会掩盖全极化的优点且无法确定准确适合的权值。因此，需要一种或几种能够保留全极化信息的数据融合方法，且可以适应不同极化散射机制的目标体。为此，本文选择了几种已经应用于图像融合领域的方法(主成分分析、拉普拉斯金字塔、多尺度小波变换)，结合FP-GPR数据，进行对比，选择适合不同极化散射机制目标体的方法。

主成分分析(Principal Component Analysis,PCA)是一种基于统计特征的多维正交线性变换方法，通常用于信号的特征提取和数据降维[17]。PCA最初是Pearson[18]在1901年提出的，并已广泛应用于人脸识别，网络入侵检测和图像降噪。例如，Safayani等人[19]提出了扩展二维PCA来提高人脸识别的准确性和时间；Hadri等人[20]将PCA和模糊PCA结合，可以保留来自网络流量数据的最相关信息；Qiu等人[21]使用PCA来提取SAR图像目标分类中的特征。

拉普拉斯金字塔(Laplacian Pyramid,LP)是将原图像分解成不同空间分辨率的子图像，高分辨率子图像放到下面，低分辨率子图像放到上面，形成金字塔形[22]。Burt等人[23]提出的LP最早是用于图像的压缩处理，广泛应用于图像融合、图像检测等。例如，Mao等人[24]提出了一种多方向联合平均的LP算法，用来加强融合图像边缘信息；Liang等人[25]提出了由拉普拉斯锐化和LP组成的联合边缘检测器，加强了图像的边缘检测。

20世纪80年代发展起来的小波变换(Wavelet Transform,WT)技术，是一种图像的多分辨率多尺度分解[26]，可以将信号分解为多个高频信号和1个低频信号，如今已经广泛应用于图像检测、数字视频水印、语音增强等方面。例如，Kumar等人[27]用WT来分解高低频信号来进行图像边缘检测；Liu等人[28]将离散WT与奇异值分解结合方法进行了改进，在数字视频水印的不可感知性与鲁棒性方面有了显著提高;Tasmaz[29]利用双树复数WT对语音效果进行了增强。

本文在实验室中选择3个典型目标体获取其全极化探地雷达数据集，然后使用主成分分析、拉普拉斯金字塔、多尺度小波变换方法分别获得这3个典型目标体的融合图像，利用瞬时振幅为主、梯度为辅来比较加权平均融合和这3种融合的效果，并且找到适合未知散射机制目标体的方法。之后对冰裂缝进行实际数据测量，利用未知散射机制的融合方法与加权融合方法比较，最后得出结论。

2 全极化探地雷达数据融合理论

平面波传播时，极化是指电场强度方向与入射面的关系。电场方向位于入射面中，为平行极化(H)；电场方向垂直于入射面，为垂直极化(V)。传统商用探地雷达是一种单极化探地雷达，它可以传输垂直极化的波形并接收相同的波形(获取VV数据)。全极化GPR将在发射H极化和V极化波形之间交替进行，同时接收H极化和V极化信号，从而获得VV,HH,VH和HV数据。由于在GPR系统中发射天线和接收天线的作用可以互换，所以VH数据等于HV数据。天线配置如图1所示，其中T表示发射天线，R表示接收天线。

为了进行全极化探地雷达数据融合，我们可以使用VV,HH和VH探地雷达系统作为不同的传感器来探测同一地下目标，获得VV,HH,VH数据集。这3个数据集就是做数据融合的基础。加权平均融合公式为

其中，SVV,SHH和SVH是由FP-GPR获得的数据矩阵。ZVV,ZHH和ZVH是平衡SVV,SHH,SVH三者贡献的权重系数[11]。

2.1 主成分分析数据融合

主成分分析可以将n维数据转换为m维数据，其中m维数据称为主成分。主成分可以代表原始数据，因为它具有最大的方差，代表着原始数据的大多数特征。因此，我们可以使用寻找主成分的方法来去除3个极化数据的多余相同部分，并融合不同部分。

为了对FP-GPR数据S={SVV,SHH,SVH}使用主成分分析，首先需要将SVV,SHH和SVH3个矩阵转化为一个FP-GPR数据矩阵。全极化探地雷达PCA数据融合的步骤如下：

(1) 假设SVV,SHH和SVH矩阵维度都是M×N，而FP-GPR矩阵为X=[XVVXHHXVH]。其中，XVV,XHH,XVH分别是由SVV,SHH,SVH转化而得的，3个长度为M×N的一维向量，比如XVV是将SVV的N列数据按2～N的顺序依次放到第1列数据下，使其变成一维向量。

(2) 计算数据X的均值，其中是X的平均值

(7)C是一个长度为M×N的一维向量，再将其重新转化成一个维度为M×N的矩阵F。F就是数据融合后的FP-GPR矩阵。

2.2 拉普拉斯金字塔数据融合

图像的LP分解能够将原图像分解到各个不同的频带范围，这样数据融合就可以针对不同的分解层获得的不同频带，采用不同的融合规则，对VV,HH,VH图像分解获得的相同分解层进行融合，再对其进行图像重构，可以获得具备几个原图像的所有细节特点的融合图像。

图1 全极化探地雷达系统中的天线配置Fig.1 Antenna configurations used in a FP-GPR system

下面以VV数据为例进行LP分解。

同理，可以得到HH与VH的拉普拉斯金字塔。对同一层级的用相同融合规则进行数据融合，都是低频成分，而其余的LP成分是高频信息。高频部分代表细节信息，低频部分代表轮廓信息，所以对于融合法则的选取是对对应处的每个像素点，低频部分平均处理，保留3种极化方式的目标体反射情况，高频部分则是选取最大值，取3幅图像中最能代表该目标体的细节部分，增加了融合图像的细节反射。最后得到融合图像的拉普拉斯金字塔。

在数据融合之后，进行图像重构，得到融合后的图像。由LP重建融合后的图像，从LP顶层开始逐层由上到下进行逆推式(10)，就可以得到重建的高斯金字塔。由式(11)表示，最后，一直代入得到，就是所得LP数据融合图像。

2.3 多尺度小波变换数据融合

小波变换的多尺度特性，即在大尺度时，观察到图像的轮廓，在小尺度的空间里，则可以观察图像的细节[26]。根据mallat快速算法[26]，可以做二维WT，将VV,HH,VH进行图像分解，得到不同频带范围的图像，将对应频带的部分进行融合，再重构即可得到FP-GPR数据融合的结果。以VV为例，进行mallat小波变换方法分解，如式(12)所示

其中，j表示小波分解的层数，VVj表示小波分解j层的原始近似图像，j=0时为原始图像。h和g分别表示低通滤波器和高通滤波器，*表示共轭反转。V Vj+1为对VVj在水平方向降采样后做低通滤波，再在垂直方向降采样后做低通滤波，代表着图像的低频部分，也是下一级的原始图像；为对VVj数据先在水平方向降采样后做低通滤波，再在垂直方向降采样后做高通滤波，代表着图像在垂直方向的高频信息；为对VVj数据先在水平方向降采样后做高通滤波，再在垂直方向降采样后做低通滤波，代表着图像在水平方向的高频信息；为对VVj数据在水平方向降采样后做高通滤波，再在垂直方向降采样后做高通滤波，代表着图像对角线的高频信息。进行多层之后就是将图像的低频部分继续进行上述操作，最后分解J次之后，会得到关于VV的一个低频部分，和3J个高频部分。

之后与LP相似，对HH和VH进行相同层数的小波分解，取3组数据相同频带的分解部分进行融合，如进行融合，不同部分不同规则进行融合，对高频部分取最大值，低频取平均值，得到融合图像，最后进行小波逆变换进行重构，如式(13)所示，最后一直推到F0即为数据融合后的结果。

3 室内实验和讨论

为了检测3种FP-GPR数据融合方法的效果和对应3种极化散射机制地下目标体的成像效果，本文选取了3种典型的金属目标体，分别是对应单次散射的金属板，二次散射的金属二面角，体散射的金属多分支散射体，如图2所示。板的长宽分别是35 cm和20 cm。二面角是由两个摆成角度为90°的金属板构成的，在两个板的交际处具有最大的二次散射。多分支散射体的长度大约为40 cm。将这3种目标体分别埋进实验室的干沙中，深度分别为23 cm,32 cm,25 cm，并由实验室搭建的FP-GPR系统进行测量[30]，如图3所示。测线共有99个测点，测点间距为1 cm，样本点数为1024个，选择的频带范围为800 MHz～4.5 GHz。

本文对FP-GPR数据进行预处理，用傅里叶逆变换将频率域数据转变为时间域，都进行去直达波处理。处理之后3种目标体图像如图4—图6所示。对于单散射目标体，如图4所示，VV极化和HH极化可以清晰反映出目标体的位置，而VH极化虽然较杂乱，但仍可以看到目标体的位置；对于二次散射目标体，如图5所示，VV极化可以清晰看到目标体的整体反射剖面形状，而HH极化更突出的是二次散射的位置，即二面角的交界处，VH极化更多体现的是目标体反射剖面的左侧；对于体散射目标体，如图6所示，VV极化和HH极化的目标体反射剖面不清晰，而VH极化可以更清晰看到目标体的反射情况。这表明需要数据融合，将3种极化的不同特点进行融合，还需要找到更适合不同极化散射机制的方法。

图2 金属目标体Fig.2 Metallic targets

图3 FP-GPR系统Fig.3 FP-GPR system

图4 板的FP-GPR图像Fig.4 FP-GPR image of plate

图5 二面角的FP-GPR图像Fig.5 FP-GPR image of dihedral

图6 多分支散射体的FP-GPR图像Fig.6 FP-GPR image of multi-branch scatterer

最后，本文进行数据融合。将加权平均融合、PCA融合、LP融合和WT融合的适应性进行了比较，通过比较寻找更适应3种不同极化散射机制的方法。对于加权平均融合，为了方便对权值系数都选择了1做平均，而其余3种融合方法根据前面理论部分所述进行计算，LP和WT均选择分解4层。为了突出目标体位置，本文使用了偏移成像的方法来处理四种融合后图像，偏移公式[4,9]如式(14)所示

其中，vrms表示均方根速度，r和t表示测点(x,0)和图像点s(x,z)之间的距离和传播时间，θ表示传播方向与垂直方向夹角，SF表示融合后数据。本文使用瞬时振幅为主梯度为辅来判别这几种数据融合方法的效果，瞬时振幅由希尔伯特变换求得。瞬时振幅代表着整个图像中振幅的大小，可以突出异常体的位置。梯度公式由文献[10]中改变而得，如式(15)所示

图7 板偏移后融合图像Fig.7 Fusion images of plate after migration

图8 板偏移后融合瞬时振幅图Fig.8 Fusion instantaneous amplitude images of plate after migration

图9 板偏移后融合梯度图Fig.9 Fusion gradient images of plate after migration

其中，M与N表示行数与列数，SF为融合图像的数值。梯度公式表示图像值与附近其他值差异性的大小，可以突出目标体处的异常。结果如图7—图15所示。为了方便比较，色标以及梯度和瞬时振幅里的垂向轴是相同的。在肉眼无法辨别峰值大小的时候，做了3种目标体瞬时振幅和梯度的最大值表格，如表1、表2所示。

对于单散射目标体，如图7所示，4种方法都可以体现出地下目标体的位置。因为在单次散射目标体中，VV在FP-GPR中占有主体地位，所以另3种融合方法在图像上与加权融合方法无明显区别。它们都能看到目标体的反射剖面，以及极化信息是单次散射。为了更进一步比较，以瞬时振幅为主梯度为辅的方法，如图8、图9所示，我们看到PCA融合和WT融合瞬时振幅较大，LP融合次之。在梯度中WT融合比PCA和LP融合的值大很多，所以单次散射目标体选择WT融合方法。

对于二次散射目标体，如图10所示，加权平均融合只有较弱的反射剖面，其它3种有较强剖面反射，在二面角交际处有更突出能量，这个位置体现着二次散射的位置即二面角的极化信息。进一步比较，选择瞬时振幅和梯度，如图11、图12所示，瞬时振幅中，在二次散射位置即交际处，PCA融合有单一最突出的峰值，LP融合是两个突出峰值，而WT融合是多个峰值导致信息较杂乱，而加权平均融合在交际处和图像两侧有较弱的幅值；在梯度图像中，加权平均融合、PCA融合和LP融合的梯度峰值在二次散射的位置均较小，而WT融合的梯度峰值形成多个峰值无法突出目标体。以瞬时振幅为主的原则，二次散射目标体选取PCA融合。

图10 二面角偏移后融合图像Fig.10 Fusion images of dihedral after migration

图11 二面角偏移后融合瞬时振幅图Fig.11 Fusion instantaneous amplitude images of dihedral after migration

图12 二面角偏移后融合梯度图Fig.12 Fusion gradient images of dihedral after migration

对于体散射目标体，如图13所示，加权平均融合的目标体反射剖面不清晰，是因为融合中VV和HH掩盖了只有VH可以清晰得到的反射信号，而其余3种方法都有明显反射剖面。进一步比较，选择瞬时振幅和梯度，如图14、图15，在瞬时振幅中，加权融合有几个较小峰值，无法突出目标体信息，而PCA融合和LP融合有两个突出位置的峰值，如表1所示，与WT融合的峰值大小近似，WT融合却有多个峰值突出，都体现体散射目标体的位置；在梯度图像中，加权融合最小，无法反映目标体位置；PCA融合与LP融合较大，可以反映目标体的位置，如表2所示LP融合的梯度值稍高于PCA融合；而WT融合信息较为杂乱，无法反映出目标体的有效信息。综上，体散射目标体选取LP融合方法。

对于未知散射机制的目标体，本文通过瞬时振幅选取在不同散射机制中表现的效果均较好的数据融合方法。单次散射目标体中，PCA融合和WT融合在目标体位置具备最大的瞬时振幅，LP融合次之；二次散射目标体中，PCA,LP和WT融合均在二次散射位置具有最大瞬时振幅，三者相比PCA融合具有唯一突出峰值；体散射目标体中，PCA与LP融合瞬时振幅能量大小相似，WT融合略低。综上，面对未知目标体时我们优先选择PCA融合。

4 冰裂缝实验与讨论

在冰面上本文使用FP-GPR系统进行冰层裂缝探测如图16所示，得到冰层的FP-GPR数据，将其进行数据融合处理。测线垂直于冰裂缝摆放，位于测线中点处。测线共有101个测点，测点间距为1 cm，样本点数为1024个，选择的频带范围为1～4 GHz。进行预处理得到的FP-GPR图像如图17所示，由于天线馈电点过高，在天线内部传到冰表面有一定的距离，到8 ns处即为冰的表面。VV极化和HH极化在图像中间，即50 cm处均有反射双曲线，且由上至下能量逐渐衰减，这与冰裂缝存在位置吻合；且VH极化在50 cm处存在明显的相位断层，推断是由冰层裂缝引起的，本文需要数据融合的方法将3种不同极化方式进行融合，得到具备3种极化信息且分辨率更高的冰裂缝FP-GPR图像。

图13 多分支散射体偏移后融合图像Fig.13 Fusion images of multi-branch scatterer after migration

图14 多分支散射体偏移后融合瞬时振幅图Fig.14 Fusion instantaneous amplitude images of multi-branch scatterer after migration

图15 多分支散射体偏移后融合梯度图Fig.15 Fusion gradient images of multi-branch scatterer after migration

表1 瞬时振幅最大值Tab.1 Max of instantaneous amplitude

表2 梯度最大值Tab.2 Max of gradient

图16 冰裂缝实验Fig.16 Ice crack experiment

图17 冰裂缝FP-GPR数据图像Fig.17 FP-GPR image of ice crack

由于前文已经选择了PCA融合作为未知散射机制目标体的融合方法，所以本文对3种极化数据进行加权平均融合、PCA融合，并进行偏移成像处理。通过瞬时振幅和梯度进行对比，证明PCA融合的可靠性。为方便比较，图像色标均保持一致。结果如图18所示，PCA融合可以清晰得到冰裂缝的位置图像，在图像50 cm,8 ns处有较为明显的能量集中，数据融合突出了冰裂缝具备最强反射的位置，且随着深度的增加能量逐渐衰减，比加权平均融合有更好的图像清晰度。由瞬时振幅图像可见，如图19所示，PCA融合在冰裂缝位置具备明显的峰值，且明显高于加权平均融合。同时比较梯度图像，如图20所示，PCA融合在冰裂缝位置有更高的梯度值。由此得出，在实际数据及未知散射机制的目标体上，PCA融合可得到很好的融合图像，比原始3种极化图像具备更突出的反射特征，强于加权融合方法。

5 结论

本文将主成分分析、拉普拉斯金字塔和小波变换分别应用于全极化探地雷达数据融合中，均可以取得较好的融合效果，融合图像具备更清晰的几何图像特征。此外，引入了瞬时振幅和梯度对3种融合方法与已有的加权平均融合进行比较，得出了3种方法均强于加权平均融合，以及3种不同融合方法适应不同散射机制目标体的结论。对于已知散射机制目标体，WT融合更适合单散射目标体，PCA融合更适合二次散射目标体，LP融合更适合体散射目标体；对于未知散射机制目标体，选择了PCA融合。之后进行未知散射机制的冰裂缝实际数据采集，可以发现，PCA融合可以得到很好的实际数据融合图像，且强于加权平均融合方法。综上，融合方法集合了目标体的几何图像特征和不同的散射极化信息，可以得到更高分辨率的GPR图像。

图18 冰裂缝偏移后融合图像Fig.18 Fusion images of ice crack after migration

图19 冰裂缝偏移后融合瞬时振幅图Fig.19 Fusion instantaneous amplitude images of ice crack after migration

图20 冰裂缝偏移后融合梯度图Fig.20 Fusion gradient images of ice crack after migration