1 引 言
干涉合成孔径雷达(interferometric synthetic aperture radar,InSAR)能够提取大范围地区的数字高程模型(digital elevation model,DEM),是应用十分广泛的雷达对地测绘技术[1, 2]。常规InSAR系统受相位欠采样的影响,难以实现对陡峭地区的高程提取,为此国内外学者提出了多通道InSAR处理方法[3, 4, 5, 6],其中多频InSAR技术是近年来的研究热点之一。
在多频InSAR信号处理研究方面,大体有3类方法。第1类方法是对多频干涉相位进行融合处理,其中,文献[7]提出基于多频相干性分析和卡尔曼滤波的解缠相位融合和DEM生成方法;文献[8]提出的多频最小二乘相位解缠是对传统单频最小二乘解缠算法的推广; 而文献[9]则提出了利用相关信号子空间方法(coherent signal subspace method,CSSM)进行多频干涉相位融合的技术,这类方法能够降低相位解缠误差,但是难以解决干涉相位的频谱混叠问题。第2类方法是基于观测相位的算术性质进行解模糊处理,其中,文献[10]最早提出应用中国余数定理、投影方法和线性组合的方法来解算干涉相位的整数模糊周期的思想;文献[11]也设计了一种基于中国余数定理的多波段 InSAR 相位解缠方法,这类方法能够较好地解决干涉相位频谱混叠问题,但处理结果易受噪声影响。第3类是基于统计学的方法,文献[12]首先提出了一种基于最大似然准则的多频相位估计方法;文献[13]也提出了一种利用多通道干涉数据,不经过相位解缠直接进行地理编码和高程反演的最大似然估计方法;文献[14]则在其研究中将传统相位解缠的初步结果应用于多频干涉最大似然估计来重构DEM,这些方法因严格依赖于局部像素的多频观测相位而对噪声的适应性较差。文献[15—16]提出的基于马尔科夫随机场模型的多频干涉最大后验概率估计方法,虽然利用相位观测矩阵的结构特性提高了对噪声的适应性,但是算法运算量成倍增加。
为了提高多频InSAR系统在噪声背景下对陡峭地区的高程提取能力,本文提出了一种基于梯度重建的多频InSAR相位解缠方法,该方法利用目标像素的邻域统计特性,实现了对相位梯度的最大似然估计,并将梯度重建结果用于支持相位解缠。
2 相位梯度重建算法定义InSAR系统获取的二维相位观测矩阵位于大小为N×M的支撑域s上,即{i,j}∈s,且i=0,1,…,N-1,j=0,1,…,M-1。进一步地,定义s的子集s0为{i,j}∈s0,且i=0,1,…,N-2,j=0,1,…,M-2;并定义s的另外两个子集:s1为{i,j}∈s1,且i=0,1,…,N-2,j=0,1,…,M-1,以及s2为{i,j}∈s2,且i=0,1,…,N-1,j=0,1,…,M-2。
设由工作在l个不同载频的雷达传感器对同一目标场景获取SAR干涉图像对,并分别通过图像配准等处理,得到l个通道的InSAR干涉相位观测值φl[i,j],即
式中,W{·}表示相位缠绕算子;α是与InSAR系统有关的参数;h[i,j]表示目标像素对应的实际高程值;fl和ωl分别表示第l个通道的工作频率及该通道包含的相位噪声。在无其他先验知识的情况下,各通道纵向和横向相位梯度初始估计分别为 当且仅当相位梯度满足式(3)所示的无旋场性质时,才能获得无模糊的相位解缠结果[17] 根据式(1),在陡峭地区当相邻像素间的高程变化超过一定界限时,式(2)所示的相位梯度初始估计就无法正确反映像素间的绝对相位变化量,需要通过相位梯度重建估计出式(2)丢失的梯度模糊数,从而使干涉相位矩阵重新满足式(3)所示的无旋场特性。 2.1 多频InSAR相位梯度的最大似然估计原理在InSAR系统中,干涉信号一般服从联合圆对称高斯分布,干涉相位概率密度为[18]
该式反映了在一定相干系数γ及等效视数n条件下,实际观测相位φ与其理论值φ之差的概率分布情况,其中,Γ表示伽马函数,F表示高斯超几何函数,由于余弦函数的周期性,导致基于单通道干涉相位会得到存在模糊的高程估计结果。而干涉相位梯度为相邻像素间的相位差,不失一般性以纵向梯度为例,其概率密度函数可表示为[4] 式中,⊗代表卷积运算。需要注意的是,式(5)与式(4)的物理意义相类似,它反映了相位梯度观测值Δ1φ[i,j]与相位梯度理论值Δ1φ[i,j]之差的概率分布情况。对多频InSAR系统,假设各通道间载频满足f1>f2>…>fl,以最高载频的相位梯度为参考,其他通道相位梯度可根据式(6)进行换算,从而统一各通道间的相位梯度变化量
经换算后得到的相位梯度
,其概率密度函数转化为
式中,Δ1φref表示参考载频的相位梯度理论值。当通道间满足独立性条件时,多频相位梯度的联合概率密度函数即为
式中,
。分析式(8),在相干系数γ、等效视数n和梯度观测向量Δ1Φ[i,j]已知的情况下,多频相位梯度的联合概率密度只取决于Δ1φref[i,j]。而注意到Δ1φref[i,j]即为相位梯度重建中的待估计量,故对Δ1φref[i,j]作参数估计的似然函数可表示为
于是,可根据式(9)实现对参考载频相位梯度Δ1φref的最大似然估计
在信号处理中,
的具体估计方法如下:在梯度估计范围内,通过穷举法计算待估计量Δ1φref[i,j]对应的似然函数L(Δ1φref[i,j])取值,根据式(10),使似然函数取值最大的估计量即为相位梯度重建结果
。
分析式(4)~(10),参考载频的初始相位梯度概率密度函数以2π为周期,经换算后的其他载频相位梯度概率密度函数则扩展为以2ρkπ(k=1,2,…,l)为周期,而各通道联合概率密度函数的相位梯度模糊周期会进一步扩展,具体分析如下:对{ρk}k=1l中各载频比保留适当的小数位数q,并将各个ρk乘以10q得到整数集合{ρ′k}k=1l,ρ′k∈Z,则式(8)对相位梯度的无模糊估计区间将扩展为2mπ,m可由集合的最小公倍数(least common multiple,LCM)确定
2.2 基于邻点集的相位梯度重建方法由于干涉相位会受到噪声影响,根据式(10)进行相位梯度重建,需要的独立观测通道数目较多,而在系统体积、功耗以及天线设计等约束下,可增加的频点个数十分有限。对于宽带系统,可以通过频谱分割扩充频点数目[12],而对于非宽带系统,本文进一步研究了基于邻点集的最大似然估计方法,以解决在噪声背景下对有限观测频点系统的相位梯度重建。
分析2.1节中的最大似然估计方法,每个像素的梯度估计过程都独立于其邻域像素,故可通过对邻域信息的联合处理来提高估计的稳健性,定义像素(i,j)对应的邻点集Ωi,j为
式中,T1、T2为整数,决定了邻点集大小。以场景目标像素(i,j)为原点建立xOy局部平面坐标系,其邻点集Ωi,j对应的局部高程模型可由d阶二维曲面多项式来拟合,即 此时,局部曲面模型的多项式系数集合{βm,n}就决定了相位梯度大小,具体来说,纵向和横向梯度所依赖的系数子集分别为 故对各像素的相位梯度重建,可转化为对其邻点集局部曲面多项式系数的最大似然估计,仍以纵向梯度估计为例,可表示为 在信号处理中,通过穷举法寻找系数组合的最优估计
,将导致运算效率严重下降,本文3.1节将针对此问题给出一种快速估计方法。
由于局部曲面坐标系原点定义在待估计像素所处的位置,故结果向量
中的一次项系数β1,0就是待估计像素的纵向相位梯度重建结果
。同理,可求得横向相位梯度重建结果
。关于参考载频的相位梯度模糊数可由式(16)求得
对于多频InSAR系统提供的多通道独立观测量,采用本文提出的基于邻点集的相位梯度重建算法实现对参考载频的梯度模糊数估计,并将其用于修正式(2)求得的初始相位梯度,得到新的相位梯度先验矩阵如下
由于相位解缠算法大多是基于相位梯度信息实现的,故可将式(17)的结果直接用于支持相位解缠处理,从而可解决陡峭地区的干涉相位欠采样问题,其基本流程如图 1。
|
| 图 1 基于梯度重建的相位解缠方法基本流程图 Fig. 1 Algorithm flow chart of phase unwrapping based on multi-frequency phase gradient reconstruction |
步骤1:在状态空间内随机产生一组初始解Β1,并令最优解
=Β1,计算当前的似然函数值L(
),设初温T=T0(事先均匀抽样一组状态,并将各状态似然函数值的方差作为初温)。
步骤2:在当前解的邻域内,按照均匀分布的概率密度随机采样,产生一组新的候选解Β′1,计算其对应的似然函数值L(B′1),从而得到似然函数增量ΔL=L(B′1)-L(
)。
步骤3:根据Metropolis准则[20]确定状态转移概率p,即
接下来产生[0, 1]之间的随机数ε,并作如下判定,若p>ε,则令最优解
=Β′1;否则最优解
保持不变。
步骤4:在当前温度T,重复一定次数的扰动和接受过程,即反复执行步骤2~步骤3。
步骤5:采用对数退火降温方式逐渐降低温度T,即T(i)=T0/lg(1+i),i表示循环次数。
步骤6:重复执行步骤2~步骤5,直至算法搜索到的最优值L(
)连续若干步保持不变,即可输出最优系数组合
,从而得到纵向相位梯度的最优解。
与文献[13, 21]中提出的利用多频InSAR数据直接进行高程估计的传统方法相比,本方法对各通道间绝对相位误差不敏感,因而可以在无控制点情况下完成相位梯度重建。相比之下,多频InSAR直接高程估计虽无需相位解缠,但必须通过定标手段来消除通道间的绝对相位误差,而且往往需要外部输入的粗精度DEM[13],才能取得比较理想的高程反演结果。
为分析该方法对相对相位误差的敏感性,下面以双频InSAR系统为基础研究其参数选择问题。相对相位误差不仅会导致梯度估计错误,还可能使估计区间内出现两个或多个似然函数值十分接近的估计量。对不同的载频比及梯度估计区间,通过加入不同程度的噪声,分别统计系统能够容忍的最大噪声标准差以及出现模糊估计的概率,所得结果如图 2。
|
| 图 2 不同载频比情况下双频InSAR系统的解模糊性能和噪声适应性分析 Fig. 2 Ambiguity and noise compatibility analysis for dual-frequency InSAR under different frequency ratio |
可见,双频InSAR系统在某些载频比条件下,出现模糊估计的概率会达到或接近100%,此时虽然系统对噪声的适应性也较强,但是在载频比设计上要尽量避开这些点,而应该优先选择出现模糊估计概率较低,且对噪声适应性较强的点,如图 2(a)中载频比为1.67,以及图 2(b)中载频比为1.75的点,对于具有更多频点的InSAR系统,在载频选择上也应使两两频点之间符合上述原则。此外,梯度估计区间大小也会影响估计性能,估计区间范围设定越大,梯度解模糊能力越强,但是对噪声的适应能力也会降低。考虑到实际InSAR系统较高的信号采样率,随地形陡峭程度不同,干涉相位梯度的模糊数一般在±2~±4之间,因而可以通过限定梯度解区间达到抑制噪声的目的,而且相位解缠算法对梯度重建环节产生的个别错误估计也有一定的适应能力,从而可以进一步提高该方法的稳健性。
4 试验验证及分析 4.1 仿真数据首先,利用仿真DEM数据来验证本文提出的多频InSAR相位梯度重建算法的有效性。设仿真高程模型为陡峭多山地形,且InSAR系统采用编队飞行方式,由主卫星发射雷达信号,主、从卫星同时接收形成干涉图,该系统分别工作在C波段(载频5.39 GHz)以及X波段(载频9.65 GHz),其他主要系统参数如表 1。
| 轨道高度/km | 下视角/(°) | 基线长度/m | 干涉图分辨率/m |
| 500 | 36 | 400 | 20×20 |
在上述参数下生成的SAR图像平均信噪比约8dB,图 3分别给出了C/X两个波段的干涉相位,其中方框标注为相位欠采样区,此时由单频干涉条纹将无法得到正确的相位解缠结果。
|
| 图 3 仿真陡峭地区的不同波段干涉相位图 Fig. 3 Interferometric phase of different wave band for simulated rugged terrain |
以X波段载频为参考,分别采用逐点估计法和基于邻点集的估计方法,对该双频干涉相位进行联合处理。以横向(距离向)相位梯度为例,图 4(a)显示了原始相位梯度的欠采样现象,图 4(b)为对这组数据的逐点最大似然估计结果,可见相位梯度基本得到正确重建,但是受相位噪声影响而产生了诸多病态估计点,而采用基于邻点集的梯度估计方法所得结果如图 4(c)所示,其中病态估计量大大减少,表明该方法增强了对相位噪声的适应性。此时,将重建的相位梯度用于支持相位解缠,得到如图 4(d)所示的正确相位解缠结果。
|
| 图 4 参考载频原始相位梯度及其重建结果 Fig. 4 Original range phase gradient and corresponding reconstruction result for reference wave band |
接下来,利用实际DEM数据进一步分析该算法的高程提取精度。所用地形数据来源于美国科罗拉多地区Isolation Peak附近山脉,为了便于与高程提取结果作比较,图 5(a)给出了将试验地区参考高程值重采样到SAR成像斜平面后的结果。设InSAR系统采用与之前相同的工作模式和工作频率,其他参数如表 2。图 5(b)和图 5(c)分别为这组参数下的低频和高频干涉条纹图,而图 5(d)和图 5(e)则分别给出了参考载频(X波段)的横向(距离向)和纵向(方位向)相位梯度的欠采样部分。考虑到InSAR系统中各种非理想因素产生的相位噪声,干涉相干系数分布情况如图 7(e),其中相干系数主要分布在0.4~0.9之间,在少数发生叠掩的局部地区,相干系数仅为0.1~0.25左右。可见所选试验地区存在大量分散及连续的地形突变区域,而这部分区域所对应的相干系数总体偏低,随之产生的较为严重的相位噪声,进一步加剧了干涉条纹的模糊程度,从而增加了对该地区的高程提取难度。
|
| 图 5 Isolation Peak地区的干涉量及相位欠采样区 Fig. 5 Interferograms and under sampling profile for mountainous area around Isolation Peak |
| 轨道高度/km | 下视角/(°) | 基线长度/m | 干涉图分辨率/m |
| 600 | 30 | 360 | 3.75×7.5 |
|
| 图 7 Isolation Peak地区数据仿真结果 Fig. 7 Simulation result for DEM of Isolation Peak |
此时,仅根据参考载频的干涉条纹将无法完成相位解缠,而采用文献[7]中提出的多频相位解缠算法进行处理,得到如图 7(a)所示的结果,可见受相位欠采样影响,解缠结果中存在相位不连续现象,并引起轻微的高度模型失真。接下来,利用高、低频干涉相位采样值,通过逐点最大似然估计实现对参考载频的相位梯度重建,再经后续处理会得到如图 7(b)所示的高程提取结果,此时由相位欠采样所引起的高度模型失真基本被消除,但是由于梯度重建结果中存在较多病态解,因而会影响高程提取精度,并导致图 7(b)中出现相位解缠无法处理的孤岛。最后,通过对局部高程模型的二维曲面建模,采用基于邻点集的估计方法实现多频InSAR相位梯度重建,再经后续处理所得结果如图 7(c),可见该方法有效抑制了孤岛现象,并获得了较为理想的高程提取结果。
对文献[7]方法与本文方法的相位解缠结果,分别评估其绝对相位误差,统计结果如图 6。由于本文方法通过梯度重建有助于解决相位欠采样问题,使相位解缠过程中的误差积累得到有效控制,因此大部分像素的绝对相位误差都比较小,相位解缠精度得到明显改善。
|
| 图 6 相位解缠的绝对相位误差比较 Fig. 6 Comparison over absolute phase error for phase unwrapping |
多频InSAR对地观测作为对传统InSAR技术的拓展,能够提高陡峭地区的高程提取能力,针对多频InSAR联合处理易受相位噪声影响的问题,本文提出了一种基于梯度重建的多频InSAR相位解缠方法。该方法通过基于邻点集的局部曲面建模来增强相位梯度估计的稳健性,并用梯度重建结果支持相位解缠,从而可解决陡峭地区的相位欠采样问题。文中进一步从频点选择、估计区间限定等方面研究了改善噪声适应性的途径。试验表明,该方法对噪声具有较好的稳健性,能够有效抑制局部病态解的出现,适于实现对陡峭地区的高程提取。
| [1] | LI Deren, YANG Jie. Principle and Applications of Extracting DEM from SAR[J]. Journal of Geodesy and Geodynamics, 2002, 22(2): 1-6. (李德仁, 杨杰. 从卫星雷达提取地面高程信息的原理与应用[J]. 大地测量与地球动力学, 2002, 22(2): 1-6.) |
| [2] | LU Lijun, ZHANG Jixian, WANG Teng. A DEM Mapping Method Assisted by External DEM with High Resolution InSAR Data in Complex Terrain Area[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(4): 459-463. (卢丽君, 张继贤, 王腾. 一种基于高分辨率雷达影像以及外部DEM辅助的复杂地形制图方法[J]. 测绘学报, 2011, 40(4): 459-463.) |
| [3] | JIN Guowang. Research on Key Technologies of InSAR Topographic Surveying and Mapping[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(5): 668. (靳国旺. InSAR地形测绘若干问题研究[J]. 测绘学报, 2011, 40(5): 668.) |
| [4] | FORNARO G, PAUCIULLO A, SANSOSTI E. Phase Difference-based Multichannel Phase Unwrapping[J]. IEEE Transactions on Image Processing, 2005, 14(7): 960-972. |
| [5] | FERRALUOLO G, MEGLIO F, PASCAZIO V, et al. DEM Reconstruction Accuracy in Multichannel SAR Interferometry[J]. IEEE Transactions on Geoscience Remote Sensing, 2009, 47(1): 191-201. |
| [6] | ZHANG Hongmin, JIN Guowang, XU Qing, et al. Application of Chinese Remainder Theorem in Phase Unwrapping for Dual-baseline InSAR[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 770-777. (张红敏, 靳国旺, 徐青,等. 中国余数定理在双基线InSAR相位解缠中的应用[J]. 测绘学报, 2011, 40(6): 770-777.) |
| [7] | LANARI R, FORNARO G. Generation of Digital Elevation Models by Using SIR-C/X-SAR Multifrequency Two-pass Interferometry: the Etna Case Study[J]. IEEE Transactions on Geoscience and Remote Sensing, 1996, 34(5): 1097-1114. |
| [8] | VINOGRADOV M V, ELIZAVETIN I V. Phase Unwrapping Method for the Multifrequency and Multibaseline Interferometry[C] //IEEE International Geoscience and Remote Sensing Symposium. Seattle:[s.n.], 1998: 1103-1105. |
| [9] | ZHANG Sen, TANG Jinsong, YANG Hailiang, et al. Interferometric Phase Fusion and DEM Reconstruction Method with CSSM for Multi-frequency InSAR Systems[J]. Geomatics and Information Science of Wuhan University, 2010, 35(3):328-333. (张森, 唐劲松, 杨海亮, 等. 利用CSSM进行多频InSAR相位融合及DEM重构[J].武汉大学学报:信息科学版, 2010, 35(3): 328-333.) |
| [10] | XU W, CHANG E C. Phase Unwrapping of SAR Interferogram with Multi-frequency or Multi-baseline[C] //IEEE International Geoscience and Remote Sensing Symposium. Pasadena:[s.n.], 1994: 730-732. |
| [11] | JIN Guowang, ZHANG Hongmin, XU Qing, et al. Phase Unwrapping Algorithm with CRT for Multi-band InSAR[J]. Journal of XiDian University, 2011, 38(6): 110-117. (靳国旺, 张红敏, 徐青, 等. 多波段InSAR的 CRT相位解缠方法[J]. 西安电子科技大学学报, 2011, 38(6): 110-117.) |
| [12] | PASCAZIO V, SCHIRINZI G. Estimation of Terrain Elevation by Multifrequency Interferometric Wide Band SAR Data[J]. IEEE Signal Processing Letters, 2001, 8(1): 7-9. |
| [13] | EINEDER M, ADAM N. A Maximum-Likelihood Estimator to Simultaneously Unwrap, Geocode, and Fuse SAR Interferograms from Different Viewing Geometries into One Digital Elevation Model[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(1): 24- 36. |
| [14] | WANG Wei. Height Reconstruction Using Multifrequency InSAR Data[D]. Harbin :Harbin Institute of Technology, 2006. (王伟. 多频率InSAR高程信息重建算法研究[D]. 哈尔滨:哈尔滨工业大学, 2006.) |
| [15] | FERRAIUOLO G, PASCAZIO V, SCHIRINZI G. Maximum a Posteriori Estimation of Height Profiles in InSAR Imaging[J]. IEEE Geoscience and Remote Sensing Letters, 2004, 1(2): 66-70. |
| [16] | BIOUCAS-DIAS J M, VALADAO G. Phase Unwrapping via Graph Cuts[J]. IEEE Transactions on Image Processing, 2007, 16(3): 698-709. |
| [17] | GHIGLIA D C, PRITT M D. Two-dimensional Phase Unwrapping: Theory, Algorithms, and Software[M]. New York: John Wiley & Sons, 1998. |
| [18] | JONG-SEN L, HOPPEL K W, MANGO S A, et al. Intensity and Phase Statistics of Multilook Polarimetric and Interferometric SAR Imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 1994, 32 (5): 1017-1028. |
| [19] | KANG Lishan, XIE Yun, YOU Shiyong, et al. Non-Numerical Parallel Algorithms (Vol.1): Simulated Annealing Algorithm[M]. Beijing: Science Press, 1994.(康立山, 谢云, 尤矢勇, 等. 非数值并行算法(第一册): 模拟退火算法[M]. 北京:科学出版社,1994.) |
| [20] | LI S Z. Markov Random Field Modeling in Image Analysis[M]. New York: Springer-Verlag, 2001. |
| [21] | PASCAZIO V, SCHIRINZI G. Multifrequency InSAR Height Reconstruction through Maximum Likelihood Estimation of Local Planes Parameters[J]. IEEE Transactions on Image Processing, 2002, 11(12): 1478-1489. |


