一种宽带合成孔径雷达子带干涉数据处理方法

摘要

一种宽带合成孔径雷达子带干涉数据处理方法，该方法有六大步骤：步骤一：主副图像的自适应精密配准及副图像的重采样；步骤二：子带频谱位置和带宽的优化；步骤三：主副图像的带通滤波；步骤四：子带干涉图和子带差分干涉图的生成；步骤五：参考面和地形高程(DEM)的相位补偿；步骤六：子带差分干涉图的堆栈处理。本发明在大多数情况下无需相位解缠，即可获取绝对雷达干涉相位，从而可以通过转换直接得到数字地形图或者地表形变图，消除了由于相位解缠错误引起的误差扩散；它特别适用于高山地区的数字地形测绘和大形变速率的滑坡及矿山塌陷监测。

说明书

技术领域

本发明涉及一种宽带合成孔径雷达子带干涉数据处理方法，属于合成孔径雷达干涉技术领域。它可以在大多数情况下无需相位解缠，即可获取绝对雷达干涉相位，从而可以通过转换直接得到数字地形图或者地表形变图。 

背景技术

合成孔径雷达干涉是一门根据高分辨率雷达图像的相位数据来提取地面目标三维空间信息的技术。无论是利用雷达干涉进行地形数字高程测量还是地表形变测量，都有一个相位解缠的过程。因为雷达干涉图的相位都被一个2π的角度所缠绕，所谓相位解缠，就是通过累积雷达干涉图像元间的相位梯度，从而获得像元绝对相位的过程。这一过程要能够得到正确结果，对于地形测量来说，一个必要条件是，在空间域上，相邻两像元间的相位差的绝对值不得大于π；对于利用堆栈差分干涉图的时序分析方法的形变测量，除了前面那个必要条件必须得到满足外，在时间域上，同一像元相邻两时间点上的相位差的绝对值亦不得大于π。在实际工程中，尤其是在高山地区的地形测量时，或者是对具有较大形变速率梯度的滑坡和矿山的形变场进行干涉测量时，上面提到的相位差的绝对值不得大于π的条件常常是不能得到满足的，因而造成测量结果出现很大误差，有时甚至会出现无解的情况。根据雷达干涉理论，要解决这一相位解缠的问题，就是要想办法提高干涉相位测量模糊度。对于地形测量而言，其模糊度定义为： 

$H_{2\pi }=-\frac{\mathrm{\lambda r}}{{2B}_{\perp }}\sin \left(\theta \right)---\left(1\right)$

上式中，λ是雷达载波波长，r是地面目标与雷达的距离，θ是雷达入射角，B_⊥为雷达干涉空间基线，H_2π即高程模糊度，它表示干涉相位差为2π时对应的地面高程。 

对于形变测量而言，其模糊度定义为： 

${\mathrm{\Delta r}}_{2\pi }=\frac{\lambda }{2}---\left(2\right)$

上式中，λ是雷达载波波长，△r_2π即形变模糊度，它表示差分干涉相位差为2π时对应的地表形变量。 

显然，从（1）和（2）式可以看出，增大雷达载波波长可以增大干涉测量的高程和形变模糊 度，但这会增大硬件制造难度和成本，且效果也并不显著，因为几倍或十来倍的增大还是不能满足工程实际需要。本专利申请利用宽带雷达信号的特性，在频域内，首先将雷达图像分裂为几个子带图像，然后将对应的雷达子带图像进行干涉运算，将高低两个子带干涉图再次进行干涉，这时的模拟波长被扩大了一个很大的倍数，这个倍数取决于原来的载波中心频率与子带中心频率的2倍的比值，从而可以将原来的厘米级的波长扩大到米级，雷达干涉技术中的难题—相位解缠的困难度将会得到缓解，甚至可以避免相位解缠这一过程。 

发明内容

1.目的：本发明的目的是提供一种宽带合成孔径雷达子带干涉数据处理方法，它克服了现有技术的不足，在大多数情况下无需相位解缠，即可获取绝对雷达干涉相位，从而可以通过转换直接得到数字地形图或者地表形变图，消除了由于相位解缠错误引起的误差扩散。 

2.技术方案：本发明是一种宽带合成孔径雷达子带干涉数据处理方法，该方法具体步骤如下： 

步骤一：主副图像的自适应精密配准及副图像的重采样 

由于子带干涉技术主要是应用在高山地区的高程测量或者是具有较大形变速率的滑坡和矿山塌陷测量，在主副图像进行配准时，不仅要考虑雷达干涉系统（成像几何，干涉几何，雷达运行轨道等等）引起的地面同一目标在主副图像上的位置的偏移，还要考虑地面目标因较大的高程或较大的形变位移引起的主副图像上的位置偏移。在本发明中，副图像利用Sinc函数进行重采样，Sinc函数的位移量由下式决定： 

Offset(i,j)=(Δi,Δj)=((Δi，Δj)_i,j)_System+((Δi,Δj)_i,j)_Object  (3) 

上式中，offset(i,j)表示了副图像中位于第i行第j列的像元相对于主图像上同一像元位置的偏移量，它是系统偏移和目标本身因高程或形变引起的偏移量的和。offset(i,j)由副图像s中位于(i,j)的像元的邻域与相应的主图像p中像元的邻域的协方差函数取最大值时的位移量决定，见下面的(4)和(5)式 

$c_{\mathrm{ps}}{(n,m)}_{i,j}=\underset{k=-N/2}{\overset{k=N/2}{\Sigma }}\underset{l=-M/2}{\overset{l=M/2}{\Sigma }}(p(k,l)-\tilde{p})·(s(k+n,l+m)-\tilde{s})---\left(4\right)$

$\left(\begin{array}{c}\mathrm{\Delta i}=n_{\max }\\ \mathrm{\Delta j}=m_{\max }\end{array}\right)---\left(5\right)$

上两式中，p和s分别表示主副图像，c表示在主副图像中以第i行第j列为中心的N行 M列的两个窗口的协方差函数，n和m分别表示行和列偏移量，k和l表示窗口中参与相关计算像元的行列位置，(Δi，Δj)表示当c取得最大值时的行列号，即（nmax，mmax）。 

步骤二：子带频谱位置和带宽的优化 

在理想的情况下，我们希望子带带宽尽可能的宽，这样能使雷达图像分辨率损失不要太大。另一方面，我们希望上下子带图像没有相同的频率成分，否则这些相同的频率成分会以噪声的形式出现在二次干涉图中。而根据信号处理理论，理想带通滤波器是不可实现的，也就是说，输入信号在带通滤波器的阻带内只能受到有限的衰减，因此，我们希望上下子带的位置相距越远越好。显然，两个理想状态是不可能同时得到满足的。一个折中的选择，是上下子带的带宽bw都等于雷达原始信号带宽B的三分之一，上下子带位于f_c±f₀的位置上，即 

$\left(\begin{array}{c}\mathrm{bw}=\frac{B}{3}\\ f_0=\frac{B}{3}\end{array}\right)---\left(6\right)$

在这样的优化选择下，模拟波长与载波波长的比值为： 

$\frac{{\lambda }^{\prime }}{\lambda }=\frac{{3f}_c}{2B}---\left(7\right)$

一般情况下： 

$\left(\begin{array}{c}{f}_0=\frac{B-\mathrm{bw}}{2}\\ {\lambda }^{\prime }=\frac{f_c}{B-\mathrm{bw}}·\lambda \end{array}\right)---\left(8\right)$

在上面三式中，B是雷达脉冲带宽，fc表示载波频率，λ表示载波波长，bw表示上下子带带宽，f₀表示上下子带中心频率对于载波频率的偏移量，λ′表示子带差分干涉图的模拟波长。 

步骤三：主副图像的带通滤波 

主副雷达SLC图像通过上下子带带通滤波器滤波，生成各自的上下子带复数图像。上子带带通滤波器的归一化中心频率为0.333，归一化带宽为0.333；下子带带通滤波器归一化中心频率为0.666，即-0.333，归一化带宽为0.333。为了减少频谱能量的泄漏，带通滤波器的频谱需加窗，窗口函数选用凯撒窗(Kaiser),其频谱表达式如下： 

$W(f;\beta ,F)=\frac{I_0\left(\beta \sqrt{1-{(2f/F)}^2}\right)}{I_0\left(\beta \right)}-\frac{F}{2}<f<\frac{F}{2}---\left(9\right)$

上式中，F为通带带宽，I₀是零阶贝塞尔函数，β是控制窗口通带宽和阻带衰减折中关系的一个常数，通常取3～6之间。 

图2,图3和图4分别显示了原始雷达SLC图像及通过带通滤波器滤波以后的上下子带图像的斜距(Slant Range)方向上的频谱。 

步骤四：子带干涉图和子带差分干涉图的生成 

主图像的上下子带图与对应的副图像的上下子带图进行干涉，产生上子带干涉图和下子带干涉图。将上下子带干涉图再次进行干涉，生成子带差分干涉图。记P和S分别为主副图像，PS为它们的子带差分干涉图，则 

$\mathrm{PS}=(P_{\mathrm{up}}·S_{\mathrm{up}}^{*})·{(P_{\mathrm{low}}·S_{\mathrm{low}}^{*})}^{*}---\left(10\right)$

上式中，P_up，P_low，S_up，和S_low分别表示由主副图像分裂出来的上下子带图，PS表示由这四副子带图像构成的子带差分干涉图。 

步骤五：参考面和地形高程(DEM)的相位补偿 

在本发明中，地球的WGS84椭球体模型作为干涉测量的参考面，如果要进行差分干涉测量，地球的WGS84椭球体表面还需叠加所测区域的DEM，参考面和地形高程(DEM)引起的模拟干涉相位要从子带差分干涉图中去除。为此，首先需计算主副图像中各像元与雷达的距离R，它满足下面的非线性方程组： 

$\left(\begin{array}{c}R=\sqrt{(\overrightarrow{R_s}-\overrightarrow{R_t})·(\overrightarrow{R_s}-\overrightarrow{R_t})}\\ \frac{2}{\mathrm{\lambda R}}(\overrightarrow{V_s}-\overrightarrow{V_t})·(\overrightarrow{R_s}-\overrightarrow{R_t})=0\\ \frac{x_t^2}{{(a+h)}^2}+\frac{y_t^2}{{(a+h)}^2}+\frac{z_t^2}{b^2}=1\end{array}\right)---\left(11\right)$

上式中，R_s为雷达位置向量，R_t为像元位置向量，即(x_t,y_t,z_t)，V_s和V_t分别为雷达和像元目标运动速度，h为像元高程，若仅补偿参考面相位，取h=0，a和b是地球椭球体模型的长短轴。计算出每一个像元分别在主副图像中与雷达的距离R以后，通过像元在主副中的距离的差ΔR，便可以计算出子带差分干涉图中的参考面和地形高程的模拟相位： 

$\phi =-\frac{4\pi }{{\lambda }^{\prime }}·\mathrm{\Delta R}---\left(12\right)$

上式中的波长应为子带干涉后的模拟波长。图5是从两张配准后的雷达图像开始到子带差分干涉图的形成的主要过程的示意图。利用(1)式和(2)式所表达 的高程和形变模糊度便可将子带差分干涉图转化为地表高程或地表形变场。 

步骤六：子带差分干涉图的堆栈处理 

a)数字高程图(DEM)的生成 

子带差分干涉图的高程模糊度虽然被提高了很多倍，但同时也将干涉相位的误差的方差亦扩大了几近相同的倍数，单张子带差分干涉图的相位会含有更大的噪声。因为干涉图的相位通过相位模糊度的因子正比于高程，故可以将多张子带差分干涉图组成一个堆栈，然后用最小二乘法估计每一个像元的高程： 

$\underset{h}{\mathrm{Min}}\Sigma {(h-\frac{\Delta {\phi }_i}{2\pi }{H}_{2\pi ,i})}^2---\left(13\right)$

上式中，Δφi表示像元i在堆栈中某一干涉图里的相位，H_2π,i表示该点的高程模 

糊度，h表示该点高程的估计值。 

b)地表形变场的测量 

由于子带差分干涉图的形变(位移)模糊度很大，故由差分相位转换为形变量时几乎不用进行相位解缠，但为了消除大气影响、高程误差及其它的各种噪声，还是有必要利用一组子带差分干涉图(堆栈)的时序分析，特别是短基线集法，来求得精确的地表形变场及其动态变化。子带差分干涉图的时序分析方法同传统的全带宽数据干涉时序分析法基本是相同的，唯一不同的是高相干点的选取方法。本发明中，利用了两个判据来确定候选的高相干点： 

$\left(\begin{array}{c}\frac{\sigma }{\mu }<c\\ \mathrm{SCR}>d\end{array}\right)---\left(14\right)$

上式中，μ和σ是像元在堆栈时间范围内强度值的均值和方差，c是对方差与均值的比值设定的一个阈值，它基本上能衡量该像元在堆栈时间范围内幅值的稳定与否，一般取值为0.2左右。SCR是像元的信杂比(Signal to Clutter Ratio)，在子带差分干涉图中，可以进行时序分析的像元的信杂比必须大于一个阈值d,一般取值25dB。选取了候选点后，利用一般的周期图法或者短基线集法都可以获得地表形变场的速率，线性与非线性形变量。如果借助人工反射器，大尺度（数十厘米甚至米级）的突变和非线性形变的测量在本发明中都不会受到限制。 

3.优点及功效：本发明提供了一种宽带合成孔径雷达子带干涉数据处理方法，其优点为： 

(1)本方法利用宽带雷达的子带信号干涉，成百倍地增大了模拟载波波长和干涉测量高程模糊度，在地形复杂的高山地区，不用极易引起错误的解缠过程即可获取雷达干涉绝对相位，从而得到地表的数字高程。 

(2)本方法利用宽带雷达的子带信号差分干涉，可以直接获取地表数十厘米甚至米级的形变，极大地扩大了传统差分干涉雷达的测量范围。 

(3)本方法利用宽带雷达的序列子带信号差分干涉，可以获取地表大尺度的非线性形变，这在山体滑坡形变监测中极为有效。 

(4)由于本方法构造了两个不同中心频率的子带信号，因而可以像双频GPS接收机一样，测量出电离层的电子密度变化，消除该变化引起的相位误差，提高干涉测量精度。 

附图说明

图1.宽带雷达SLC图像的频带及分裂子带示意图 

图1中，f_c为雷达载波频率，f₀为子带中心频率相对于雷达载波频率的偏移量，bw表示带宽。 

图2.宽带雷达SLC图像距离向频谱 

图3.上子带图像距离向频谱 

图4.下子带图像距离向频谱 

图5.子带干涉流程图 

图6.树坪滑坡地区的TerraSAR-X雷达图像 

图7.树坪滑坡体上的8号和9号反射器滑动位移量三种测量方法结果比较图 

图8.本发明流程框图 

具体实施方式

以长江三峡大坝地区树坪滑坡形变监测为例，说明本发明在实际工程应用中的具体操作步骤。见图8，本发明是一种宽带合成孔径雷达子带干涉数据处理方法，该方法具体步骤如下： 

步骤一：数据选取与主副图像的自适应精密配准 

1.子带干涉技术适用性考量 

根据工程目的和工作区的地形地貌，地面植被覆盖程度以及地表形变情况判断是否启用子带干涉技术。对于以生成DEM为目的，只适用于高山地区，对于形变监测，只适用于具有大形变速率的地区。无论何种目的，皆要求地面被监测(或测绘)目标能具有较高的雷达相干特性。对于用时序分析法的高相干点的形变监测，则要求地面被监测 (或测绘)目标能具有较高的雷达信杂比(SCR>25dB)，必要时选用和安装人工反射器。三峡树坪滑坡地形变化大，常年农作物和果树覆盖，滑坡非稳定区域滑动量很大，传统干涉雷达无法监测出该滑坡的形变信息，故宜采用安装人工反射器和子带差分干涉方法。在树坪滑坡体上安装了14个人工反射器，另有4个安装在滑坡体外，作为测量和验证的基准。 

2.数据选取 

选取多幅覆盖工作区的星载宽带合成孔径雷达SLC图像，带宽要求100MHz以上。为此目的，选取了德国TerraSAR-X星载X波段3米分辨率的2009年4月22号至2010年4月9号的29幅SLC数据。数据幅宽30公里，长60公里，带宽150MHz，数据覆盖时间跨距一年，卫星重复过境周期为11天。雷达主要数据参数和扫描日期见表1和表2所示： 

表1：选用TerraSAR-X雷达主要参数表 

载波频率 9.6500000+9Hz 载波波长 3.1cm 脉冲带宽 1.5000000e+8Hz 入射角 26.2542° 距离向像元大小 0.909403m 方位向像元大小 1.965121m 距离向采样率 1.6482919+8Hz 雷达扫描模式 条带 数据类型 单视复数

表2：雷达数据日期表 

20090422 20090503 20090514 20090525 20090605 20090616 20090627 20090708 20090719 20090730 20090810 20090821 20090901 20090912 20090923 20091015 20091026 20091117 20091128 20091209 20091220 20091231 20100111 20100122 20100202 20100213 20100307 20100409

3.自适应精密配准 

选取一幅图像为主图像，作为参考图像，原则上它应该位于时域(时间基线集)和空间域(空间基线集)的中心。由于本发明可以适用于大尺度形变测量，在树坪滑坡监测中又借助于人工反射器，故主图像的选取不受前面所叙的限制。为了时序分析方便，选取第一张图像，即2012年1月7号的数据为主图像。以(1)式为基础，利用(2)式和(3)式计算其余各副图像中所有像元相对于主图像的相对偏移向量，用移位的Sinc函数对各副图像进行自适应二维重采样，以达到精密配准的目的，精度为十分之一个像元。 

Offset(i,j)=(Δi，Δj)=((Δi，Δj)_i,j)_System+((Δi,Δj)_i,j)_Object  (1) 

$c_{\mathrm{ps}}{(n,m)}_{i,j}=\underset{k=-N/2}{\overset{k=N/2}{\Sigma }}\underset{l=-M/2}{\overset{l=M/2}{\Sigma }}(p(k,l)-\tilde{p})·(s(k+n,l+m)-\tilde{s})---\left(2\right)$

$\left(\begin{array}{c}\mathrm{\Delta i}=n_{\max }\\ \mathrm{\Delta j}=m_{\max }\end{array}\right)---\left(3\right)$

式（1）中的第一项表示主图像中位于第i行第j列的像元与副图像中的同名点的位置的偏差量，该偏差量是一个二维向量，即式(1)中第二项。偏差向量是由系统产生的分量和主图像中该像元对应的地面目标本身形变产生的分量构成，分别为式（1）中第三和第四项所示。 

在卫星轨道数和雷达成像时间等参数已知的情况下，可以先计算出系统产生的偏移分量。以该系统分量为基础，可以大致确定副图像中对应的像元位置。这样，对于主图像的中每一个像元（i,j），皆可以找到一个以该像元为中心的一个适当大小的窗口，比如512行和512列，同时亦可以确定副图像中对应的该窗口的位置。计算这两个窗口的协方差函数，如式（2）所示。M，N为窗口的行列数，p和s分别表示主副图像。使得协方差函数取得最大值时的m和n即为主图像中像元（i,j）在副图像中的偏移量，如式（3）所示。 

步骤二：带通滤波器的频域设计-子带频谱位置和带宽的优化 

在雷达载波频率f_c和雷达脉冲信号带宽B已知的情况下，确定带通滤波器的中心频率f₀，带宽bw以及滤波器的传递函数。首先，上下两个带通滤波器要相距尽量的远，即一个在雷达带宽频谱的最低端，另一个在最高端，这样可以最大限度地增大带通滤波器阻带衰减，从而减小子带差分干涉图中的噪声。中心频率的大小将决定子带差分干涉图的模拟波长，也将决定模拟模糊度的大小。一般取bw=B/3。由于子带带宽 缩小三倍，显然分辨率会增大三倍。bw决定以后，f₀=(B-bw)/2，子带差分干涉图的模拟波长即为雷达载波波长的k倍，k=f_c/(B-bw)。根据工作区的地形或形变信息最大值，检查子带干涉的高程模糊度或形变模糊度是否合适，否则，可以通过改变带通滤波器中心频率和带宽的大小来调整k的大小，从而确定合适的模糊度，以避免子带差分干涉图的相位解缠。在频域设计中，常将频谱对带宽B归一化，这样得到下面带通滤波器的优化设计值： 

$\left(\begin{array}{c}\mathrm{bw}=\frac{1}{3}\\ f_{\mathrm{up},0}=-f_{\mathrm{low},0}=\frac{1}{3}\\ k=\frac{{3f}_c}{2B}\end{array}\right)---\left(4\right)$

带通滤波器的窗函数，这里也即传递函数，选用(5)式所示的凯撒窗，将式（4）中的参数带入式（5），取β=6。 

$W(f;\beta ,F)=\frac{I_0\left(\beta \sqrt{1-{(2f/F)}^2}\right)}{I_0\left(\beta \right)}-\frac{F}{2}<f<\frac{F}{2}---\left(5\right)$

上式中，F为通带带宽，I₀是零阶贝塞尔函数，β是控制窗口通带宽和阻带衰减折中 

关系的一个常数，通常取3～6之间 

步骤三：主副图像的带通滤波与子带图像的生成 

子带图像的生成都是在频域里进行的。首先将所有29张配准后的图像（p，s₁,……,s₂₈）通过二维快速富里哀变换（2D-FFT）转换成频域图像（P,S₁,……,S₂₈）。这29张频域图像分别通过步骤四设计好的上下两个带通滤波器滤波，生成各自的频域里的上子带图像和下子带图像（P_up,P_low,S_{1_up},S_{1_low},……,S_{28_up},S_{28_low}）。在滤波前，子带图像要通过解调，也即频谱移位，将带通频谱移位到基带上，使得所有的上子带图像的模拟载波频率是f_c+f₀,所有的下子带图像的模拟载波频率是f_c-f₀，且所有子带图像的频谱都以其模拟载波频率为中心，以保证后续信号处理的频谱同一和结果的正确。然后再利用二维快速逆富里哀变换（2D-IFFT）将这58张频域图像转换成58张空间域图像（p_up,p_low,s_{1_up},s_{1_low},……,s_{28_up},s_{28_low}）。 

图1为宽带雷达SLC图像的频带及分裂子带示意图，图2，图3和图4分别显示了原始雷达SLC图像及通过带通滤波器滤波以后的上下子带图像斜距（Slant Range）方向上的频谱。 

步骤四：子带干涉图及子带差分干涉图的生成 

以2009年4月22号图像为主图像，其余皆为副图像的方式，生成28张上子带干涉图和28张下子带干涉图： 

$p_{\mathrm{up}}·s_{i\_\mathrm{up}}^{*},i=1,......,28---\left(6\right)$

$p_{\mathrm{low}}·s_{i\_\mathrm{low}}^{*}.i=1,......,28---\left(7\right)$

$(p_{\mathrm{up}}·s_{i\_\mathrm{up}}^{*})·{(p_{\mathrm{low}}·s_{i\_\mathrm{low}}^{*})}^{*},i=1,......,28---\left(8\right)$

如(8)式所示，每一干涉对产生一张子带差分干涉图，它是由该干涉对的两张下子带图像产生的干涉图像与两张上子带图像产生的干涉图像再次进行干涉而产生的。由于此时子带差分干涉图的带宽只有原来的三分之一(优化设计的情况下)，像元分辨率已经降低了三倍，故在做再次干涉时，应做至少3x3的多视处理，以降低斑点噪声。 

步骤五：参考面和地形高程的相位补偿 

参考面和地形高程的相位补偿过程与传统差分干涉中的相位补偿基本上是形同的。如果只做地形的数字高程(DEM)的生成，子带差分干涉图只需做参考面的相位补偿，如果是做地表形变监测，除了要做参考面的相位补偿外，还需做地形高程的相位补偿。为了与作为雷达载体的卫星的运行轨道数据相匹配，选取地球WGS84椭球体模型的表面作为参考面，椭球体表面和其上面的地形高程产生的模拟干涉相位根据卫星轨道数据，雷达成像相数据及(9)式和（10）式计算，与传统差分干涉不同，这时的模拟雷达载波波长是原始载波波长的k倍，k的值见式（11）。从子带差分干涉图中减去这个由参考面及地形高程产生的模拟干涉相位。如果是仅估计地形高程的话，令(9)式中的h为零。 

$\left(\begin{array}{c}R=\sqrt{(\overrightarrow{R_s}-\overrightarrow{R_t})·(\overrightarrow{R_s}-\overrightarrow{R_t})}\\ \frac{2}{\mathrm{\lambda R}}(\overrightarrow{V_s}-\overrightarrow{V_t})·(\overrightarrow{R_s}-\overrightarrow{R_t})=0\\ \frac{x_t^2}{{(a+h)}^2}+\frac{y_t^2}{{(a+h)}^2}+\frac{z_t^2}{b^2}=1\end{array}\right)---\left(9\right)$

上式中，R_s为雷达位置向量，R_t为像元位置向量，即(x_t,y_t,z_t)，V_s和V_t分别为雷达和像元目标运动速度，h为像元高程，若仅补偿参考面相位，取h=0，a和b是地球椭球体模型的长短轴。计算出每一个像元分别在主副图像中与雷达的距离R以后，通过像元在主副中的距离的差ΔR，便可以计算出子带差分干涉图中的参考面和地形高程的模拟相 位： 

$\phi =-\frac{4\pi }{{\lambda }^{\prime }}·\mathrm{\Delta R}---\left(10\right)$

上式中的波长应为子带干涉后的模拟波长: 

$k=\frac{{3f}_c}{2B}---\left(11\right)$

图5是从两张配准后的雷达图像开始到子带差分干涉图的形成的主要过程示意图。 

步骤六：子带差分干涉图的堆栈处理 

1.子带差分干涉图的滤波 

由于子带差分干涉图相位的高斯特性，宜采用二维高斯滤波器对子带差分干涉图的相位进行滤波，消除高斯噪声。 

2.子带差分干涉图的堆栈处理与地表形变场的测量 

由于子带差分干涉图的形变(位移)模糊度很大，故由差分相位转换为形变量时几乎不用进行相位解缠，但为了消除大气影响、高程误差及其它的各种噪声，还是有必要利用一组子带差分干涉图(堆栈)的时序分析，特别是短基线集法，来求得精确的地表形变场及其动态变化。子带差分干涉图的时序分析方法同传统的全带宽数据干涉时序分析法基本是相同的，唯一不同的是高相干点的选取方法。本发明中，利用了两个判据来确定候选的高相干点： 

$\left(\begin{array}{c}\frac{\sigma }{\mu }<c\\ \mathrm{SCR}>d\end{array}\right)---\left(12\right)$

上式中，μ和σ是像元在堆栈时间范围内强度值的均值和方差，c是对方差与均值的比值设定的一个阈值，它基本上能衡量该像元在堆栈时间范围内幅值的稳定与否，一般取值为0.2左右。SCR是像元的信杂比(Signal to Clutter Ratio)，在子带差分干涉图中，可以进行时序分析的像元的信杂比必须大于一个阈值d,一般取值25dB。选取了候选点后，利用一般的周期图法或者短基线集法都可以获得地表形变场的速率，线性与非线性形变量。如果借助人工反射器，大尺度（数十厘米甚至米级）的突变和非线性形变的测量在本发明中都不会受到限制。 

利用（12）式，可以确定所有18个人工反射器在每一张子带差分干涉图中的位置，读取这18个人工反射器在每张子带差分干涉图中的相位，并减去一个参考点的相位。在树坪滑坡中，6号人工反射器被选定为参考点。得到每一个人工反射器在每一张子带差分干涉 图中的相位后，按式（13），即可计算出每一个人工反射器在相应的干涉图时间段内的雷达射线方向上的位移量： 

$\rho =\frac{{\lambda }^{\prime }}{4\pi }·\mathrm{\Delta \phi }---\left(13\right)$

图6是树坪滑坡的TerraSAR-X雷达度图，在滑坡体上的14个人工反射器和4个滑坡体外的人工反射器清晰可见。为了验证本发明的测量准确性，这里选取了8号和9号反射器的子带差分干涉结果与地面水准和GPS的同步测量结果作比较。8号反射器位于滑坡体外稳定区域，9号反射器位于滑坡体较大形变区。图7是两个反射器的三种测量方法的结果的比较图。比较结果显示，三种测量结果十分吻合，这里要特别指出的是，5月份到7月份期间，9号反射器的位移达到近40厘米，这是用传统雷达干涉方法无法测量到的。 

3.子带差分干涉图的后续处理 

子带差分干涉图的解缠，相位-高程转换，相位-形变转换，子带差分干涉图堆栈时序分析，地理编码等处理过程完全与传统的全带宽雷达干涉数据处理相同，这几部分的数据处理不在本发明范围内。例外的仅是将(14)和(15)式模糊度表达式中的载波波长改成子带干涉模拟波长。对于地形测量而言，其模糊度定义为： 

$H_{2\pi }=-\frac{\mathrm{\lambda r}}{{2B}_{\perp }}\sin \left(\theta \right)---\left(14\right)$

上式中，λ是雷达载波波长，r是地面目标与雷达的距离，θ是雷达入射角，B_⊥为雷达干涉空间基线，H_2π即高程模糊度，它表示干涉相位差为2π时对应的地面高程。对于形变测量而言，其模糊度定义为： 

${\mathrm{\Delta r}}_{2\pi }=\frac{\lambda }{2}---\left(15\right)$

上式中，λ是雷达载波波长，Δr_2π即形变模糊度，它表示差分干涉相位差为2π时对应的地表形变量。