CN101881831A - 基于差分滤波的多波段InSAR相位解缠方法 - Google Patents

Abstract

一种基于差分滤波的多波段InSAR相位解缠方法，该方法由最长波长的干涉数据出发，在对长波长、低频率干涉图进行有效解缠的基础上，根据波长之间的比例关系，由解缠的低频率干涉图生成短波长、高频率干涉图的基准干涉图，再对短波长、高频率干涉图与基准干涉图进行差分处理得到差分干涉图，从而降低短波长、高频率干涉图的条纹频率；对差分干涉图进行滤波和相位解缠之后，再由基准干涉图和差分干涉图的解缠结果得到短波长、高频率干涉图的解缠结果；并依此类推，一直到最短波长的干涉图解缠完毕。本发明有效解决了较短波长情况下陡峭地区混叠干涉条纹的相位解缠难题，提高相位解缠的可解性和解缠精度。

Description

基于差分滤波的多波段InSAR相位解缠方法

技术领域

本发明涉及多波段InSAR(Multi-band Interferometric Synthetic ApertureRadar，多波段合成孔径雷达干涉测量)技术，特别涉及一种基于差分滤波的多波段InSAR相位解缠方法。

背景技术

利用多波段InSAR技术，可将不同波段(频率)的InSAR数据进行联合处理，提高对地形陡峭地区混叠干涉相位的解缠精度和可靠性，从而提高获取DEM(Digital Elevation Model，数字高程模型)的精度和可靠性。

对于多波段InSAR技术而言，相位解缠仍然是干涉处理的一个关键环节。传统的相位解缠方法大都基于如下假设：相邻像素的相位差不超过半个周期π。而实际上，对于地形起伏较大地区并采用较短波长、较长基线进行干涉成像的情况，很容易出现干涉相位的欠采样和干涉条纹的混叠。对于这种干涉条纹，很难甚至不能进行高效、高精度的相位解缠。

发明内容

本发明的目的是提供一种由最长波长开始的差分滤波多波段InSAR相位解缠方法，利用该方法可以由最长波长的干涉图解缠开始，把不同波段的干涉数据进行组合处理，逐步精化解缠结果，从而提高相位解缠的可靠性和精度。

为实现上述目的，本发明采用以下技术方案：

由最长波长的干涉数据出发，在对长波长、低频率干涉图进行有效解缠的基础上，根据波长之间的比例关系，由解缠的低频率干涉图生成短波长、高频率干涉图的基准干涉图，再对短波长、高频率干涉图与基准干涉图进行差分处理得到差分干涉图，从而降低短波长、高频率干涉图的条纹频率；对差分干涉图进行滤波和相位解缠之后，再由基准干涉图和差分干涉图的解缠结果得到短波长、高频率干涉图的解缠结果；并依此类推，一直到最短波长的干涉图解缠完毕。

本发明包括以下步骤：

①获取所有波长的干涉图Inter_i，其中，i＝1，2，…，N，N表示波段数，并对最长波长、最低频率干涉图Inter₁进行滤波处理；

②对滤波后的最长波长、最低频率干涉图Inter₁进行相位解缠，得到最长波长、最低频率干涉图Inter₁的解缠结果UnInter₁；

③根据较长波长、较低频率干涉图Inter_i的解缠结果UnInter_i构造较短波长、较高频率干涉图Inter_j的基准干涉图

其中，i＝1，2，…，N，N表示波段数；j＝2，…，N，且j＞i；当由最短波长、最低频率干涉图的解缠结果UnInter₁构造基准干涉图时，i＝1；

④对较短波长、较高频率干涉图Inter_j与基准干涉图

进行差分处理，得到差分干涉图DInter_j：

⑤对差分干涉图DInter_j进行滤波和相位解缠之后，再由基准干涉图和差分干涉图的解缠结果得到较短波长、较高频率干涉图的解缠结果；

⑥根据较短波长、较高频率干涉图的解缠结果，进而对更短波长、更高频率干涉图进行解缠，并依此类推，直到最短波长的干涉图解缠完毕。

所述的步骤①中，采用频率域低通滤波方法对最长波长、最低频率干涉图Inter₁进行滤波处理。

在所述的步骤①中，采用空间域矢量滤波方法对最长波长、最低频率干涉图Inter₁进行滤波处理。

在所述的步骤②中，采用质量图区域生长、直接相位积分和移动曲面相位拟合相结合的自适应相位解缠方法对最长波长、最低频率干涉图Inter₁进行相位解缠。

采用上述技术方案的本发明，能够较好将不同波段干涉相位进行组合处理，充分利用较长波长干涉图频率较低、相位解缠较容易的优势和较短波长干涉图频率较高、精度较高的优势，有效解决较短波长情况下陡峭地区混叠干涉条纹的相位解缠难题，提高相位解缠的可解性和解缠精度。

附图说明

图1为本发明的整体流程图；

图2为本发明中采用频域低通滤波方法对最长波长、最低频率干涉图Inter₁进行滤波处理的流程图；

图3为本发明中采用空间域矢量滤波方法对最长波长、最低频率干涉图Inter₁进行滤波处理的流程图；

图4为本发明中采用质量图区域生长、直接相位积分和移动曲面相位拟合相结合的自适应相位解缠方法对最长波长、最低频率干涉图Inter₁和差分干涉图DInter_j(j＝2，…，N)进行相位解缠的处理流程图；

图5为SRTM90米间距的DEM图；

图6为实验中波长＝0.18m的干涉图；

图7为实验中波长＝0.09m的干涉图；

图8为实验中波长＝0.06m的干涉图；

图9为波长＝0.18m干涉图的滤波结果；

图10为波长＝0.18m干涉图的解缠结果；

图11为根据波长＝0.18m干涉图解缠结果而构建的波长＝0.09m的基准干涉图；

图12为波长＝0.09m的差分干涉图；

图13为波长＝0.09m差分干涉图的滤波结果；

图14为波长＝0.09m干涉图的解缠结果；

图15为根据波长＝0.09m干涉图解缠结果构建的波长＝0.06m的基准干涉图；

图16为波长＝0.06m的差分干涉图；

图17为波长＝0.06m差分干涉图的滤波结果；

图18为波长＝0.06m干涉图的解缠结果。

具体实施方式

实施例1

对于不同波段的干涉数据而言，其干涉处理的难度和精度不同。当波长较长时，干涉图频率较低，干涉相位解缠容易，但获取的高程精度较低；反之，当波长较短时，理论上能够获取较高的高程精度，但此时容易出现干涉相位欠采样，其干涉质量较低，很难对干涉图进行有效滤波和相位解缠处理，无法进行有效地干涉测量和高程信息获取。

为了充分利用不同波段干涉数据的各自优势，本发明将不同波段的干涉数据进行组合处理，用较长波长干涉图的解缠结果指导较短波长干涉图的相位解缠，提高了干涉相位的可解性和解缠精度，提高获取DEM的能力和精度。

具体地说，本发明由最长波长的干涉数据出发，在最长波长、最低频率干涉图进行有效解缠的基础上，根据波长之间的比例关系，由最低频率干涉图的解缠结果构建较短波长、较高频率干涉图的基准干涉图，再对较短波长、较高频率干涉图与基准干涉图进行差分处理得到差分干涉图，从而降低较短波长、较高频率干涉图的条纹频率，并解决因地形起伏较大造成的相位欠采样问题；对差分干涉图进行滤波和相位解缠之后，再由基准干涉图和差分干涉图的解缠结果得到更短波长、更高频干涉图的解缠结果；并依此类推，一直到最短波长的干涉图解缠完毕，其具体的实施步骤如下，如图1所示：

①获取所有干涉图Inter_i，其中，i＝1，2，…，N，N表示波段数，并采用频率域低通滤波方法对最长波长、最低频率干涉图Inter₁进行滤波处理。其中，频率域低通滤波方法的流程图如图2所示，在滤波时首先对最长波长、最低频率干涉图Inter₁进行二维离散傅里叶变换，变换到频域后采用低通滤波器进行频域滤波，再进行二维离散傅里叶逆变换得到干涉图的滤波结果。

②对滤波后的最长波长、最低频率干涉图Inter₁，可采用质量图区域生长、直接相位积分和移动曲面相位拟合相结合的自适应相位解缠方法，先由质量图区域生长方法对高质量的干涉相位进行解缠，再采用移动曲面相位拟合方法完成整幅干涉图的相位解缠。该相位解缠流程如图4所示。

具体地说，为了对滤波后的最长波长、最低频率干涉图Inter₁进行高质量的相位解缠，需要先计算干涉质量图，干涉质量图的计算可采用伪相干值计算公式：

其中，ρ′(x，y)表示对应于干涉图像元(x，y)处的伪相干值；

表示相应处的干涉相位值；2k+1表示计算窗口的宽度。

对干涉质量图进行区域生长时，需首先给定干涉质量阈值，然后由人工给定的干涉质量较好的种子像元出发，根据质量图生长出干涉质量较好的区域，并对该生长区域采用直接相位积分方法进行相位解缠。

直接相位积分方法的步骤如下：

(i)计算干涉图相邻像元间的相位差分：

(ii)计算差分修正值：

(iii)直接相位积分，即：

在进行移动曲面相位拟合时，可采用二次曲面进行解缠相位的拟合：

移动曲面相位拟合过程如下：

(a)给定待拟合相位点P(x₀，y₀)；

(b)确定点P(x₀，y₀)的邻域范围。为了简化计算过程，可以点P(x₀，y₀)为中心，给定一个7×7的矩形窗口作为邻域范围；

(c)将选取的邻域范围内的点坐标(x_i，y_i)规划到以P(x₀，y₀)为原点的坐标系，即按下式计算：

$\left\{\begin{array}{c}{\tilde{x}}_i=x_i-x_0\\ {\tilde{y}}_i=y_i-y_0\end{array}\right.---\left(6\right)$

此时有：

(d)计算邻域内属于非屏蔽区的点数n，如果点数n＞6，列误差方程，求解移动曲面多项式的系数。移动曲面系数的计算公式为：

X＝(A^TPA)^-1A^TPL                                         (8)

其中：

$A=\left[\begin{array}{cccccc}{\tilde{x}}_1^2& {\tilde{x}}_1{\tilde{y}}_1& {\tilde{y}}_1^2& {\tilde{x}}_1& {\tilde{y}}_1& 1\\ {\tilde{x}}_2^2& {\tilde{x}}_2{\tilde{y}}_2& {\tilde{y}}_2^2& {\tilde{x}}_2& {\tilde{y}}_2& 1\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ {\tilde{x}}_n^2& {\tilde{x}}_n{\tilde{y}}_n& {\tilde{y}}_n^2& {\tilde{y}}_n& {\tilde{y}}_n& 1\end{array}\right];L=\left[\begin{array}{c}\mathrm{\φ}({\tilde{x}}_1,{\tilde{y}}_1)\\ \mathrm{\φ}({\tilde{x}}_2,{\tilde{y}}_2)\\ .\\ .\\ .\\ \mathrm{\φ}({\tilde{x}}_n,{\tilde{y}}_n)\end{array}\right];P=\left[\begin{array}{c}\frac{1}{\left|{\tilde{x}}_1\right|+\left|{\tilde{y}}_1\right|}\\ \frac{1}{\left|{\tilde{x}}_2\right|+\left|{\tilde{y}}_2\right|}\\ .\\ .\\ .\\ \frac{1}{\left|{\tilde{x}}_n\right|+\left|{\tilde{y}}_n\right|}\end{array}\right];X=\left[\begin{array}{c}A\\ B\\ .\\ .\\ .\\ F\end{array}\right];$

(e)根据移动曲面系数求解P(x₀，y₀)的解缠相位值：

(f)循环(a)-(e)步，直至屏蔽区内所有点都采用移动曲面拟合出相位解缠值。

需要说明的是，除此之外，也可采用枝切法、最小二乘法等其它方法进行相位解缠，这些均为本领域普通技术人员所熟知的技术。

③由于在基线长度、基线水平角均相同的条件下，较长波长、较低频率干涉图Inter_i，i＝1，2，…，N，和较短波长、较高频率干涉图Inter_j j＝2，…，N，且j＞i之间对应像素的解缠干涉相位微分存在如下关系：

$\frac{{\mathrm{d\φ}}_{\mathrm{ai}}}{{\mathrm{d\φ}}_{\mathrm{aj}}}=\frac{{\mathrm{\λ}}_j}{{\mathrm{\λ}}_i}---\left(9\right)$

即任意两幅干涉图上对应的解缠干涉相位微分值与波长成反比。由此，可利用较长波长、较低频率干涉图Inter_i的解缠结果UnInter_i按其相位微分关系构造较短波长、较高频率干涉图Inter_j的基准干涉图

其中j＝2，…，N；j＞i，逐像素的干涉相位值

为：

${\mathrm{UnInter}}_j^0=\frac{\mathrm{UnInt}{\mathrm{er}}_i}{{\mathrm{\λ}}_j}\·{\mathrm{\λ}}_i.---\left(10\right)$

为了显示基准干涉图的构建效果，可对其进行模2π处理。需要说明的是，在公式(9)和(10)中，i＝1，2，…，N，N为波段数，且当由最短波长、最低频率干涉图的解缠结果UnInter₁构造基准干涉图时，i＝1。

④在步骤③的基础上，对较短波长、较高频率干涉图Inter_j与基准干涉图

进行差分处理，得到差分干涉图DInter_j：

这样，通过差分处理，可将较短波长、较高频率干涉图Inter_j由较高频率转化为较低频率，从而降低干涉相位滤波的难度和干涉相位解缠的难度。

⑤对差分干涉图DInter_j进行滤波和相位解缠，得到解缠结果UnDInter_j，最后再由基准干涉图和差分干涉图的解缠结果UnDInter_j计算得到较短波长、较高频干涉图的解缠结果：

${\mathrm{UnInter}}_j={\mathrm{UnDInter}}_j+{\mathrm{UnInter}}_j^0---\left(11\right)$

⑥根据较短波长、较高频率干涉图的解缠结果，进而对更短波长、更高频率干涉图进行解缠，并依此类推，到最短波长的干涉图解缠完毕。

实施例2

本实施例与实施例1不同的是，在本实施例的步骤①中，采用空间域矢量滤波方法对最长波长、最低频率干涉图Inter₁进行滤波处理。

具体地说，如图3所示，它包括以下步骤：

(a)将干涉相位值映射为矢量空间中的单位矢量。设

表示干涉图中某像元处含噪声的干涉相位值，将其映射到矢量空间中的单位矢量后为：

(b)把矢量空间的干涉图分解为正弦分量和余弦分量

(c)分别对正弦分量sin(φ(x，y)+n(x，y))和余弦分量cos(φ(x，y)+n(x，y))进行均值滤波处理；

(d)由滤波后的正弦分量

和余弦分量

计算滤波后的干涉图，令

的计算公式为：

<math><mrow><mover><mi>&amp;phi;</mi><mo>^</mo></mover><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mfenced open='{' close=''><mtable><mtr><mtd><mfrac><mi>&amp;pi;</mi><mn>2</mn></mfrac></mtd><mtd><mi>a</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mn>0</mn><mo>,</mo><mi>b</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>></mo><mn>0</mn></mtd></mtr><mtr><mtd><mfrac><mrow><mn>3</mn><mi>&amp;pi;</mi></mrow><mn>2</mn></mfrac></mtd><mtd><mi>a</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mn>0</mn><mo>,</mo><mi>b</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>&lt;</mo><mn>0</mn></mtd></mtr><mtr><mtd><mi>a</mi><mi>tan</mi><mrow><mo>(</mo><mfrac><mi>b</mi><mi>a</mi></mfrac><mo>)</mo></mrow></mtd><mtd><mi>a</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>></mo><mn>0</mn><mo>,</mo><mi>b</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>&amp;GreaterEqual;</mo><mn>0</mn></mtd></mtr><mtr><mtd><mn>2</mn><mi>&amp;pi;</mi><mo>+</mo><mi>a</mi><mi>tan</mi><mrow><mo>(</mo><mfrac><mi>b</mi><mi>a</mi></mfrac><mo>)</mo></mrow></mtd><mtd><mi>a</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>></mo><mn>0</mn><mo>,</mo><mi>b</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>&lt;</mo><mn>0</mn></mtd></mtr><mtr><mtd><mi>&amp;pi;</mi><mo>+</mo><mi>a</mi><mi>tan</mi><mrow><mo>(</mo><mfrac><mi>b</mi><mi>a</mi></mfrac><mo>)</mo></mrow></mtd><mtd><mi>a</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mtext>&lt;0</mtext></mtd></mtr></mtable></mfenced><mo>,</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mrow></math>

整个矢量均值滤波可用上式表述，其中Angle()表示计算角度值：

其他技术特征与实施例1相同。

为验证本发明所达到的技术效果，特作以下实验：图5为SRTM 90米格网间距的DEM图，图6、图7、图8分别为由DEM仿真的不同波段的含噪声干涉图，表1为所采用DEM的相关参数及多波段InSAR系统参数。

                    表1

Claims (5)

1.一种基于差分滤波的多波段InSAR相位解缠方法，其特征在于：该方法由最长波长的干涉数据出发，在对长波长、低频率干涉图进行有效解缠的基础上，根据波长之间的比例关系，由解缠的低频率干涉图生成短波长、高频率干涉图的基准干涉图，再对短波长、高频率干涉图与基准干涉图进行差分处理得到差分干涉图，从而降低短波长、高频率干涉图的条纹频率；对差分干涉图进行滤波和相位解缠之后，再由基准干涉图和差分干涉图的解缠结果得到短波长、高频率干涉图的解缠结果；并依此类推，一直到最短波长的干涉图解缠完毕。

2.根据权利要求1所述的基于差分滤波的多波段InSAR相位解缠方法，其特征在于，它包括以下步骤：

①获取所有波长的干涉图Inter_i，其中，i＝1，2，…，N，N表示波段数，并对最长波长、最低频率干涉图Inter₁进行滤波处理；

②对滤波后的最长波长、最低频率干涉图Inter₁进行相位解缠，得到最长波长、最低频率干涉图Inter₁的解缠结果UnInter₁；

③根据较长波长、较低频率干涉图Inter_i的解缠结果UnInter_i构造较短波长、较高频率干涉图Inter_j的基准干涉图

其中，i＝1，2，…，N，N表示波段数；j＝2，…，N，且j＞i；当由最短波长、最低频率干涉图的解缠结果UnInter₁构造基准干涉图时，i＝1；

④对较短波长、较高频率干涉图Inter_j与基准干涉图

进行差分处理，得到差分干涉图DInter_j：

⑤对差分干涉图DInter_j进行滤波和相位解缠之后，再由基准干涉图和差分干涉图的解缠结果得到较短波长、较高频率干涉图的解缠结果；

⑥根据较短波长、较高频率干涉图的解缠结果，进而对更短波长、更高频率干涉图进行解缠，并依此类推，直到最短波长的干涉图解缠完毕。

3.根据权利要求2所述的基于差分滤波的多波段InSAR相位解缠方法，其特征在于：在所述的步骤①中，采用频率域低通滤波方法对最长波长、最低频率干涉图Inter₁进行滤波处理。

4.根据权利要求2所述的基于差分滤波的多波段InSAR相位解缠方法，其特征在于：在所述的步骤①中，采用空间域矢量滤波方法对最长波长、最低频率干涉图Inter₁进行滤波处理。

5.根据权利要求2所述的基于差分滤波的多波段InSAR相位解缠方法，其特征在于：在所述的步骤②中，采用质量图区域生长、直接相位积分和移动曲面相位拟合相结合的自适应相位解缠方法对滤波后的最长波长、最低频率干涉图Inter₁进行相位解缠。

Images