CN105182337A - 一种基于曲面后向投影算法的形变反演方法 - Google Patents

Abstract

本发明公开了一种基于曲面后向投影算法的形变反演方法，包括步骤一、将成像网格设置为地平面，进行BP成像得到形变前的主辅图像，并进行干涉处理，得到粗DEM；步骤二、将粗DEM即H₀(N,M)作为新的成像网格所在平面进行BP成像得到形变前的主辅图像，并进行干涉处理，得到新的DEM，进行多次迭代，得到精确的目标的DEM；步骤三、将迭代得到的目标的精确DEM即H(N,M)作为成像网格所在平面进行BP成像得到形变前后的主图，差分处理得到目标的形变信息△H(N,M)；本发明利用曲面BP算法进行成像，得到目标较精确的DEM，在此基础上进行二轨法的差分处理，可得到目标的形变信息，处理更为方便。

Description

一种基于曲面后向投影算法的形变反演方法

技术领域

本发明涉及雷达技术领域，具体地说，是指一种地基差分干涉合成孔径雷达(Ground-BasedDifferentialInterferometricSyntheticApertureRadar，简称GB-DInSAR)基于曲面后向投影算法(BackProjection，BP)的形变反演方法。

背景技术

GB-DInSAR是一种常用的形变监测方法，具有区域性、全天候、全天时、定点、连续、高精度监测的优点，具有很好的灵活性和可操作性，其非接触的测量方式可以在安全距离内获取被监测危险区域的形变数据，同时采集所得信息为区域性大面积的形变信息比单点的形变信息更有助于灾害的理解和预测。地基SAR的这些特点使其在形变监测领域有广阔的应用前景，获取的形变信息对研究监测目标形变机理，预测和防止灾害的发生有积极的意义。

BP成像算法是一种应用模式广泛的SAR成像算法，同时在进行干涉处理时可具备高保相精度和自配准的性能。

对于一般的形变监测目标，其较精确的DEM(数字高程模型)获取并不容易，因此在采用BP成像进行干涉测量时，一般将成像网格置于水平面，导致反演目标DEM存在移位，进而得到的形变位置等形变信息也存在畸变。

发明内容

本发明的目的是为了解决上述问题，在曲面BP算法成像的基础上能够得到较为准确的形变量等形变信息，本发明首先采用常规BP算法即成像网格设置于水平面，成像后进行干涉处理得到目标DEM，基于此DEM进行曲面BP并进一步得到新的DEM，进行迭代最后得到精度较高的目标DEM，再在此基础上进行差分处理，得到目标形变，实现一种基于曲面BP的形变反演方法。

一种基于曲面后向投影算法的形变反演方法，包括以下几个步骤：

步骤一、将成像网格设置为地平面，进行BP成像得到形变前的主辅图像，并进行干涉处理，得到粗DEM；

步骤二、将粗DEM作为新的成像网格所在平面进行BP成像得到形变前的主辅图像，并进行干涉处理，得到新的DEM，进行多次迭代，得到比较精确的目标的DEM；

步骤三、将迭代得到的目标的精确DEM作为成像网格所在平面进行BP成像得到形变前后的主图，差分处理得到目标的形变信息。

本发明的优点在于：

本发明利用曲面BP算法进行成像，得到目标较精确的DEM，在此基础上进行二轨法的差分处理，可得到目标的形变信息，该方法使得GB-DInSAR对目标进行形变监测时不必依靠目标的外部DEM，处理更为方便。

附图说明

图1是本发明的方法流程图；

图2是本发明仿真中使用的目标地形图；

图3是本发明仿真中地基SAR系统的工作示意图；

图4是本发明仿真中目标形变图；

图5是本发明仿真中得到的主辅天线回波距离向压缩后的结果图；

图6是本发明的采用常规BP算法时的BP成像网格；

图7是本发明的常规BP算法得到的主辅天线成像结果；

图8是本发明的对主辅图像进行干涉处理后得到的目标粗DEM；

图9是本发明的第一次曲面BP成像时得到的主辅图像；

图10是本发明的第一次曲面BP成像后干涉处理得到的目标DEM；

图11是本发明的第十次曲面BP成像后干涉处理得到的目标DEM；

图12是本发明的主天线形变前后曲面BP处理得到的成像结果；

图13是本发明的仿真得到的形变反演结果；

图14是本发明的理想形变与反演形变对比图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种基于曲面后向投影算法的形变反演方法，流程如图1所示，主要包括以下几个步骤：

步骤一、将成像网格设置为地平面，进行BP成像得到形变前的主辅图像，并进行干涉处理，得到粗DEM。

具体为：

第一步，BP成像，输入主辅天线回波距离压缩后的结果，输出主辅SAR图像。

(1)网格划分

在已有的主天线回波距离压缩后的结果S_r1(n,m)，以及辅天线回波距离压缩后的结果S_r2(n,m)的情况下，对于已知的方位向分辨率ρ_a和距离向分辨率ρ_r，估计目标在方位向和地距向的范围，其中，n为回波方位向采样点数，m为距离向采样点数，设置位于地平面方位向长度L和地距向长度W的BP网格，使网格覆盖目标范围，且网格高度均设置为零，网格点数为设置为N×M，其中，N为方位向网格点数，M为地距向网格点数。

(2)计算网格点斜距，进行相位补偿；

对于网格上的每个网格点A_ij，i＝1，2，…N，j＝1，2，…M，计算在每个方位时间t_p到雷达主天线的斜距R_p1及对应时延τ_p1，p＝0，1，…n，并根据方位时间t_p和时延计算其在主天线回波距离压缩后的结果S_r1(n,m)中的对应位置a_p1，b_p1，其中，a_p1为A_ij在t_p时刻对应S_r1(n,m)中的方位向点数，b_p1为其对应的距离向点数,则在t_p时刻网格点A_ij对应主天线的成像结果S_ap1[i,j]为S_r1[a_p1,b_p1]与补偿相位exp{-j2πf_cτ_p1}的乘积，同理，计算在每个方位时间t_p到雷达辅天线的斜距R_p2及对应的时延τ_p2，并根据方位时间t_p和时延计算其在辅天线回波距离压缩后的结果S_r2(n,m)中的对应位置a_p2，b_p2，其中，a_p2为A_ij在t_p时刻对应S_r2(n,m)中的方位向点数，b_p2为其对应的距离向点数，则在t_p时刻网格点A_ij对应辅天线的成像结果S_ap2[i,j]为S_r2[a_p2,b_p2]与补偿相位exp{-j2πf_cτ_p2}的乘积，具体表达式为：

$S_{apq}\[i,j\]=S_{rq}\[a_{pq},b_{pq}\]*\exp \{-j2{\mathrm{\πf}}_c{\mathrm{\τ}}_{pq}\}---\left(1\right)$

其中，f_c为雷达工作中心频率，q＝1，2，q＝1表示主天线，q＝2表示辅天线。

(3)相干累加

将上述步骤(2)中网格点A_ij所有方位时刻t_p时，相位补偿之后的主天线成像结果S_ap1[i,j]和辅天线成像结果S_ap2[i,j]分别进行累加，即得到该网格点的最后成像结果。对所有网格点执行以上操作，可得到最终的主辅天线BP成像结果，即主SAR图像S_a1(N,M)和辅SAR图像S_a2(N,M)。即有：

$S_{aq}\[i,j\]=\underset{p}{\Σ}{S}_{rq}\[a_{pq},b_{pq}\]\·\exp \{-j2{\mathrm{\πf}}_c{\mathrm{\τ}}_{pq}\}---\left(2\right)$

其中，f_c为雷达工作中心频率，q＝1，2，q＝1表示主天线，q＝2表示辅天线，i＝1，2，…N，j＝1，2，…M。

第二步，进行干涉处理，输入主辅SAR图像，输出目标的DEM。

(1)生成干涉相位

直接将由步骤一BP算法得到的成像结果的主天线SAR图像S_a1(N,M)进行共轭得到S_a1 ^*(N,M)，并乘以BP成像后得到的辅天线SAR图像S_a2(N,M)，得到干涉图S_int(N,M)的相位即为干涉相位，即

$S_{\mathrm{int}}(N,M)={S_{a1}}^{*}(N,M)*S_{a2}(N,M)---\left(3\right)$

(2)多视处理

多视处理即对干涉图S_int(N,M)进行均值滤波。首先选取一定大小C×D(如3×11)的图像窗口，则多视系数为C×D，以该窗口内所有像素的平均值作为中心像素的值，通过窗口的顺序滑动，达到平滑图像噪声的目的，最后得到的多视处理后的图像S_d[N,M]的表达式为：

$S_d\[i,j\]=\frac{1}{CD}\underset{i-(C-1)/2}{\overset{i+(C-1)/2}{\Σ}}\underset{j-(D-1)/2}{\overset{j+(D-1)/2}{\Σ}}{S}_{\mathrm{int}}\[i,j\]---\left(4\right)$

其中，i＝1，2，…N；j＝1，2，…M，对于边界点先扩展补零再处理。

(3)相位解缠

采用经典路径跟踪法的Goldstein枝切法，输入多视处理后S_d(N,M)的相位，输出解缠绕的干涉相位φ(N,M)。主要步骤为：识别残差点；生成枝切线；沿枝切线的路径积分。

(4)DEM生成

输入上一步得到解缠后的干涉相位φ(N,M)，再根据相应的目标高度反演公式(5)即可以算出每个目标点对应的高度H，将此时得到的高度记为H₀(N,M)，即输出目标的粗DEM。具体反演公式如下所示：

$H(N,M)=h+R_1(N,M)sin\[\mathrm{\α}+{\cos }^{-1}\left(\frac{\mathrm{\λ}\mathrm{\φ}(N,M)}{4\mathrm{\π}B}\right)\]---\left(5\right)$

其中，h为主天线高度，R₁(N,M)为目标到主天线斜距，α为基线水平方向的倾角，λ为雷达工作波长，φ(N,M)为解缠后的干涉相位，B为基线长度。

步骤二、将粗DEM即H₀(N,M)作为新的成像网格所在平面进行BP成像得到形变前的主辅图像，并进行干涉处理，得到新的DEM，进行多次迭代，得到比较精确的目标的DEM。

具体为：

将网格高度设为H₀(N,M)，重复步骤一，得到每个目标点对应的高度H₁(N,M)，然后将网格高度设为H₁(N,M)，重复步骤一，得到每个目标点对应的高度H₂(N,M)，以此类推，将网格高度设为H_n(N,M),重复步骤一，得到每个目标点对应的高度H_n+1(N,M)，设定迭代次数，判断n是否小于迭代次数，如果小于则继续迭代，否则，将最终迭代得到的H_n+1(N,M)设为比较精确的目标的DEM，记为H(N,M)，其中迭代次数根据精度和运算速度设定，一般取值为5～10。

步骤三、将迭代得到的目标的精确DEM即H(N,M)作为成像网格所在平面进行BP成像得到形变前后的主图，差分处理得到目标的形变信息△H(N,M)。

具体为：

第一步，BP成像，输入主天线形变前后距离压缩后的回波，输出形变前后主天线的SAR图像。

(1)网格划分

在已有主天线形变前回波的距离压缩结果S_r1(n,m)和形变后回波的距离压缩S_r1'(n,m)的情况下，对于已知的方位向分辨率ρ_a和距离向分辨率ρ_r，估计目标在方位向和地距向的范围，设置位于地平面方位向长度L和地距向长度W的BP网格(与步骤一中的L和W相等)，使网格覆盖目标范围，且网格高度为H(N,M)，网格点数设置为N×M，N为方位向网格点数，M为地距向网格点数，n为回波方位向采样点数，m为距离向采样点数。

(2)计算网格点斜距，进行相位补偿；

对于网格上的每个网格点A_ij，i＝1，2，…N，j＝1，2，…M(N为方位向网格点数，M为地距向网格点数)，计算其形变前在每个方位时间t_p到雷达主天线的斜距R_p1及对应时延τ_p1，p＝0，1，…n，并根据方位时间t_p和时延计算其形变前在主天线回波距离压缩后的S_r1(n,m)的对应位置a_p1，b_p1，其中，a_p1为形变前A_ij在t_p时刻对应S_r1(n,m)中的方位向点数，b_p1为其对应的距离向点数，则在t_p时刻网格点A_ij形变前对应的主天线成像结果S_ap1[i,j]为S_r1[a_p1,b_p1]与补偿相位exp{-j2πf_cτ_p1}的乘积，同理计算其在形变后在每个方位时间t_p到主天线的斜距R_p1'及对应的时延τ_p1'，并根据方位时间t_p和时延计算其形变后在主天线回波距离压缩后的S_r1'(n,m)的对应位置a_p1'，b_p1'，其中，a_p1'为形变后A_ij在t_p时刻对应S_r1'(n,m)中的方位向点数，b_p1'为其对应的距离向点数，则在t_p时刻网格点A_ij形变后对应的主天线成像结果S_ap1'[i,j]为S_r1'[a_p1',b_p1']与补偿相位exp{-j2πf_cτ_p1'}的乘积，具体表示如(1)式，不过，其中q＝1，1'，q＝1表示形变前，q＝1'表示形变后。

(3)相干累加

将上述步骤(2)中网格点A_ij所有方位时刻t_p，p＝0，1，…n(n为回波方位向采样点数)时，相位补偿之后的形变前成像结果S_ap1[i,j]和形变后的成像结果S_ap1'[i,j]分别进行累加，即为该网格点的最后成像结果。对所有网格点执行以上操作，可得到最终的形变前后主天线BP成像结果，即形变前的主SAR图像S_a1(N,M)和形变后的主SAR图像S_a1'(N,M)，如(2)式，不过，其中q＝1，1'，q＝1表示形变前，q＝1'表示形变后。

第二步，差分干涉处理，输入形变前后主SAR图像，输出目标的形变信息△H(N,M)。

(1)生成差分干涉相位

直接将BP算法得到的形变前的主SAR图像S_a1(N,M)进行共轭得到S_a1 ^*(N,M)，并乘以形变后的主SAR图像S_a1'(N,M)，得到干涉图S_int(N,M)的相位即为差分干涉相位，即

$S_{\mathrm{int}}(N,M)={S_{a1}}^{*}(N,M)*{S_{a1}}^{\′}(N,M)---\left(6\right)$

(2)多视处理

多视处理即对干涉图S_int(N,M)进行均值滤波。首先选取一定大小C₁×D₁(如3×3)的图像窗口，则多视系数为C₁×D₁，以该窗口内所有像素的平均值作为中心像素的值，通过窗口的顺序滑动，达到平滑图像噪声的目的，表达式为：

$S_d\[i,j\]=\frac{1}{C_1{D}_1}\underset{i-(C_1-1)/2}{\overset{i+(C_1-1)/2}{\mathrm{\Σ}}}\underset{j-(D_1-1)/2}{\overset{j+(D_1-1)/2}{\mathrm{\Σ}}}{S}_{\mathrm{int}}\[i,j\]---\left(7\right)$

其中，i＝1，2，…N；j＝1，2，…M，在一般情况下，形变导致的相位起伏不超过2π时，不需要进行相位解缠，得到的S_d(N,M)的相位即是差分干涉相位△φ(N,M)。

(3)形变反演

输入上一步得到的差分干涉相位△φ(N,M)，再根据相应的目标形变反演公式(8)即可以算出每个目标点对应的形变△H(N,M)，具体反演公式如下所示：

$\mathrm{\Δ}H(N,M)=H(N,M)-h-\sqrt{{(\sqrt{{\left(H\right(N,M)-h)}^2+{\left(r_p\right(N,M)-r)}^2}-\frac{\mathrm{\λ}\mathrm{\Δ}\mathrm{\φ}(N,M)}{4\mathrm{\π}})}^2-{\left(r_p\right(N,M)-r)}^2}---\left(8\right)$

其中，H(N,M)为目标高度，即步骤二中多次迭代得到的目标DEM，h为主天线高度，r_p为目标在地面的投影点距原点的长度，r为主天线的地面投影距原点长度，λ为雷达工作波长，△φ(N,M)为差分相位。

实施例：

本发明的一种基于曲面BP的形变反演方法，具体实施例为：

图1所示为本发明方法的整体流程图。图2所示为实施例仿真的目标地形图，坡面具体参数为方位向长度为1000m，地距向长度为200m(距原点1000m-1200m)，坡面斜度45°，最大高度为200m。

图3所示为实施例仿真中地基SAR系统的工作示意图，系统工作中心频率为16GHz，对应波长为0.0187，采用双天线同时扫描的方式以便进行干涉处理，两天线连杆长度即基线B长度为1m，基线倾角α为45°，天线绕着支点进行平行于XZ平面的扇形扫描运动，扫描角速度ω为6°/s，天线扫描范围-10°～10°，天线连杆中点与支点距离r为4m，支点高度为3m。可计算此时距离向分辨率为2.6580m，方位向分辨率为13.1300m，仿真时设地距向目标点数128，方位向目标点数256，均匀分布于坡面。

形变设置如图4所示，为在坡面上叠加的形变，即底面圆心位于(30，130，1130)，底面半径为20m，高度为0.03m的圆锥。

图5所示分别为形变前主辅图回波数据通过解线频调的方式进行距离压缩后的结果。

步骤一、将成像网格设置为地平面，进行BP成像得到形变前的主辅图像，并进行干涉处理，得到粗DEM。

具体为：

第一步、BP成像，输入主辅天线回波距离压缩后的结果，输出主辅SAR图像。

(1)网格划分

已有主天线回波的距离压缩结果S_r1(512,256)和辅天线回波的距离压缩结果S_r2(512,256)(512为回波方位向采样点数，256为距离向采样点数)，设置位于地平面方位向长度1000m和地距向长度200m(距原点1000m-1200m)的BP网格，如图6所示，网格高度均设置为零，网格点数设置为256×512，256为方位向网格点数，512为地距向网格点数。

(2)计算网格点斜距，进行相位补偿；

对于网格上的每个网格点A_ij，i＝1，2，…256，j＝1，2，…512，计算在每个方位时间t_p到雷达主天线的斜距R_p1及对应时延τ_p1，p＝0，1，…512，并根据方位时间t_p和时延计算其在主天线回波距离压缩后的结果S_r1(512,256)中的对应位置a_p1，b_p1，其中，a_p1为A_ij在t_p时刻对应S_r1(512,256)中的方位向点数，b_p1为其对应的距离向点数,则在t_p时刻网格点A_ij对应主天线的成像结果S_ap1[i,j]为S_r1[a_p1,b_p1]与补偿相位exp{-j2πf_cτ_p1}的乘积，同理，计算在每个方位时间t_p到雷达辅天线的斜距R_p2及对应的时延τ_p2，并根据方位时间t_p和时延计算其在辅天线回波距离压缩后的结果S_r2(512,256)中的对应位置a_p2，b_p2，其中，a_p2为A_ij在t_p时刻对应S_r2(512,256)中的方位向点数，b_p2为其对应的距离向点数，则在t_p时刻网格点A_ij对应辅天线的成像结果S_ap2[i,j]为S_r2[a_p2,b_p2]与补偿相位exp{-j2πf_cτ_p2}的乘积，表达式如(1)式。

(3)相干累加

将上述步骤(2)中网格点A_ij所有方位时刻t_p，p＝0，1，…512时，相位补偿之后的主天线成像结果S_ap1[i,j]和辅天线成像结果S_ap2[i,j]分别进行累加，即得到该网格点的最后成像结果。对所有网格点执行以上操作，可得到最终的主辅天线BP成像结果，即主SAR图像S_a1(256,512)和辅SAR图像S_a2(256,512)。截取照射比较完全的目标点，则S_a1[1:52,1:483]＝S_a1[103:154,15:497]和S_a2[1:52,1:483]＝S_a2[103:154,15:497]，结果如图7所示。

第二步、干涉处理，输入主辅SAR图像，输出目标的DEM。

(1)生成干涉相位

直接将BP算法得到的主图S_a1(52,483)进行共轭得到S_a1 ^*(52,483)，并乘以辅图S_a2(52,483)，得到结果S_int(52,483)的相位即为干涉相位。

(2)多视处理

多视处理即对复数干涉图S_int(52,483)进行均值滤波。选取大小3×11的图像窗口，即多视系数为3×11，以该窗口内所有像素的平均值作为中心像素的值，通过窗口的顺序滑动，达到平滑图像噪声的目的。

(3)相位解缠

采用经典路径跟踪法的Goldstein枝切法，输入多视处理后S_d(52,483)的干涉相位，输出解缠绕的干涉相位φ(52,483)。主要步骤为：识别残差点；生成枝切线；沿枝切线的路径积分。

(4)DEM生成

输入上一步得到解缠后的干涉相位φ(52,483)，再根据相应的目标高度反演公式(5)即可以算出每个目标点对应的高度H₀(52,483)，即输出目标的粗DEM。图8所示为干涉处理之后得到的DEM，解缠过程将最外一圈相位设置为零，计算时舍弃，得到此时最高点的平均误差为14.3851m，整体平均误差为6.9153m。

步骤二、将粗DEM即H₀(N,M)作为新的成像网格所在平面进行BP成像得到形变前的主辅图像，并进行干涉处理，得到新的DEM，进行多次迭代，得到比较精确的目标的DEM。具体步骤与步骤一一致，不过，在第一步的BP成像中，首先应对应于步骤一中SAR图像的截取，截取步骤一中BP成像网格，即网格点数为52×483，52为方位向网格点数，483为地距向网格点数，并且网格高度设置为H₀(52,483)。经过BP成像得到最终的主SAR图像S_a1(52,483)和辅SAR图像S_a2(52,483)，结果如图9所示。在第二步的干涉处理中，得到新的目标高度H₁(52,483)，此为一次迭代过程，得到DEM如图10所示。同样解缠过程将最外一圈相位设置为零，计算时舍弃，得到此时最高点的平均误差为9.8032m，整体平均误差为4.9572m，误差小于步骤一中得到的DEM。再将H₁(52,483)作为网格高度重复步骤一，如此，每次将新得到的目标高度作为BP网格高度按照步骤一进行10次迭代，最终可得到较高精度的目标DEM即H(52,483)，如图11所示。每次解缠过程将最外一圈相位设置为零，在最后计算时舍弃，得到此时最高点的平均误差为2.2056m，整体平均误差为2.4760m，结果表明目标DEM精度得到很大程度的提高。下表所示为每次迭代结果的误差。

表1迭代结果误差

迭代次数	整体误差均值(m)	最高点误差均值(m)	最低点误差均值(m)
				0	6.9153	14.3851	1.9437
1	4.9572	9.8032	1.9480
				2	3.8732	6.8430	1.9614
3	3.2733	5.1091	1.9727
				4	2.9377	4.0349	1.9770
5	2.7484	3.2909	1.9736
				6	2.6373	2.8269	1.9615
7	2.5719	2.6263	1.9391
				8	2.5293	2.4701	1.9035
9	2.4994	2.3261	1.8569
				10	2.4760	2.2056	1.7999

步骤三、将迭代得到的目标的精确DEM即H(52,483)作为成像网格所在平面进行BP成像得到形变前后的主图，差分处理得到目标的形变信息△H(52,483)。具体为：

第一步、BP成像，输入主天线形变前后距离压缩后的回波，输出形变前后主天线的SAR图像。

(1)网格划分

已有主天线形变前的回波的距离压缩结果S_r1(512,256)和形变后的回波的距离压缩结果S_r1'(512,256)(512为回波方位向采样点数，256为距离向采样点数)的情况下，对应于步骤一中SAR图像的截取，截取步骤一中BP成像网格，且网格高度设置为H(52,483)，即网格点数为52×483，52为方位向网格点数，483为地距向网格点数。

(2)计算网格点斜距，进行相位补偿；

对于网格上的每个网格点A_ij，i＝1，2，…52，j＝1，2，…483，计算其形变前在每个方位时间t_p到雷达主天线的斜距R_p1及对应时延τ_p1，p＝0，1，…512，并根据方位时间t_p和时延计算其形变前在主天线回波距离压缩后的S_r1(512,256)的对应位置a_p1，b_p1，其中，a_p1为形变前A_ij在t_p时刻对应S_r1(512,256)中的方位向点数，b_p1为其对应的距离向点数，则在t_p时刻网格点A_ij形变前对应的主天线成像结果S_ap1[i,j]为S_r1[a_p1,b_p1]与补偿相位exp{-j2πf_cτ_p1}的乘积，同理，计算其在形变后在每个方位时间t_p到主天线的斜距R_p1'及对应的时延τ_p1'，并根据方位时间t_p和时延计算其形变后在主天线回波距离压缩后的S_r1'(512,256)的对应位置a_p1'，b_p1'，其中，a_p1'为形变后A_ij在t_p时刻对应S_r1'(512,256)中的方位向点数，b_p1'为其对应的距离向点数，则在t_p时刻网格点A_ij形变后对应的主天线成像结果S_ap1'[i,j]为S_r1'[a_p1',b_p1']与补偿相位exp{-j2πf_cτ_p1'}的乘积，具体表示如(1)式，不过，其中q＝1，1'，q＝1表示形变前，q＝1'表示形变后。

(3)相干累加

将上述步骤(2)中网格点A_ij所有方位时刻t_p，p＝0，1，…512时，相位补偿之后的形变前成像结果S_ap1[i,j]和形变后的成像结果S_ap1'[i,j]分别进行累加，即为该网格点的最后成像结果。对所有网格点执行以上操作，可得到最终的形变前后主天线BP成像结果，即形变前的主SAR图像S_a1(52,483)和形变后的主SAR图像S_a1'(52,483)，如(2)式，不过，其中q＝1，1'，q＝1表示形变前，q＝1'表示形变后，结果如图12所示。

第二步、差分干涉处理，输入形变前后主SAR图像，输出目标的形变信息△H(52,483)。

(1)生成差分干涉相位

直接将BP算法得到的形变前的主图S_a1(52,483)进行共轭得到S_a1 ^*(52,483)，并乘以形变后的主图S_a1'(52,483)，得到结果S_int(52,483)的相位即为差分干涉相位。

(2)多视处理

多视处理即对复数干涉图S_int(52,483)进行均值滤波。首先选取一定大小3×3的图像窗口，则多视系数为3×3，以该窗口内所有像素的平均值作为中心像素的值，通过窗口的顺序滑动，达到平滑图像噪声的目的，在形变量较小时，得到S_d(52,483)的相位即是差分干涉相位△φ(52,483)。

(3)形变反演

输入上一步得到的差分干涉相位△φ(52,483)，再根据相应的目标形变反演公式(8)即可以算出每个目标点对应的形变△H(52,483)，为减小由于解缠等影响，截取最后形变△H(1:50,1:463)＝△H(2:51,20:482)，图13为形变反演结果，反演最大形变量为0.02906m，相比实际最大形变量0.02862m，相对反演误差为1.54％，图14所示为理想形变与反演得到形变的对比图，可以看到在形变位置方面，在地距向有3m的偏差。

本发明主要针对GB-DInSAR中目标精确DEM获取不易情况下形变反演的问题，首先采用常规BP算法即成像网格设置于水平面，成像后进行干涉处理得到目标DEM，再基于此采用曲面BP成像进行干涉并多次迭代获取目标较高精度的DEM，可提高最终目标形变反演的精度。通过实例分析，进一步详述了本发明方法的实施过程，实例结果显示在目标的形变反演中，反演误差可达到1.54％，验证了本发明方法的可行性。

Claims (1)

1.一种基于曲面后向投影算法的形变反演方法，包括以下几个步骤：

步骤一、将成像网格设置为地平面，进行BP成像得到形变前的主辅图像，并进行干涉处理，得到粗DEM；

具体为：

第一步，BP成像，输入主辅天线回波距离压缩后的结果，输出主辅SAR图像；

(1)网格划分

在已有的主天线回波距离压缩后的结果S_r1(n,m)，以及辅天线回波距离压缩后的结果S_r2(n,m)的情况下，对于已知的方位向分辨率ρ_a和距离向分辨率ρ_r，估计目标在方位向和地距向的范围，其中，n为回波方位向采样点数，m为距离向采样点数，设置位于地平面方位向长度L和地距向长度W的BP网格，使网格覆盖目标范围，且网格高度均设置为零，网格点数设置为N×M，其中，N为方位向网格点数，M为地距向网格点数；

(2)计算网格点斜距，进行相位补偿；

对于网格上的每个网格点A_ij，i＝1，2，…N，j＝1，2，…M，计算在每个方位时间t_p到雷达主天线的斜距R_p1及对应时延τ_p1，p＝0，1，…n，并根据方位时间t_p和时延计算其在主天线回波距离压缩后的结果S_r1(n,m)中的对应位置a_p1，b_p1，其中，a_p1为A_ij在t_p时刻对应S_r1(n,m)中的方位向点数，b_p1为其对应的距离向点数,则在t_p时刻网格点A_ij对应主天线的成像结果S_ap1[i,j]为S_r1[a_p1,b_p1]与补偿相位exp{-j2πf_cτ_p1}的乘积，同理，计算在每个方位时间t_p到雷达辅天线的斜距R_p2及对应的时延τ_p2，并根据方位时间t_p和时延计算其在辅天线回波距离压缩后的结果S_r2(n,m)中的对应位置a_p2，b_p2，其中，a_p2为A_ij在t_p时刻对应S_r2(n,m)中的方位向点数，b_p2为其对应的距离向点数，则在t_p时刻网格点A_ij对应辅天线的成像结果S_ap2[i,j]为S_r2[a_p2,b_p2]与补偿相位exp{-j2πf_cτ_p2}的乘积，具体表达式为：

S_apq[i,j]＝S_rq[a_pq,b_pq]*exp{-j2πf_cτ_pq}(1)

其中，f_c为雷达工作中心频率，q＝1，2，q＝1表示主天线，q＝2表示辅天线；

(3)相干累加

将上述步骤(2)中网格点A_ij所有方位时刻t_p时，相位补偿之后的主天线成像结果S_ap1[i,j]和辅天线成像结果S_ap2[i,j]分别进行累加，即得到该网格点的最后成像结果；对所有网格点执行以上操作，得到最终的主辅天线BP成像结果，即主SAR图像S_a1(N,M)和辅SAR图像S_a2(N,M)；即有：

$S_{aq}\[i,j\]=\underset{p}{\Σ}{S}_{rq}\[a_{pq},b_{pq}\]\·\exp \{-j2{\mathrm{\πf}}_c{\mathrm{\τ}}_{pq}\}---\left(2\right)$

其中，f_c为雷达工作中心频率，q＝1，2，q＝1表示主天线，q＝2表示辅天线，i＝1，2，…N，j＝1，2，…M；

第二步，进行干涉处理，输入主辅SAR图像，输出目标的DEM；

(1)生成干涉相位

直接将由步骤一BP算法得到的成像结果的主天线SAR图像S_a1(N,M)进行共轭得到S_a1 ^*(N,M)，并乘以BP成像后得到的辅天线SAR图像S_a2(N,M)，得到干涉图S_int(N,M)的相位即为干涉相位，即

S_int(N,M)＝S_a1 ^*(N,M)*S_a2(N,M)(3)

(2)多视处理

多视处理即对干涉图S_int(N,M)进行均值滤波，首先选取C×D的图像窗口，则多视系数为C×D，以该窗口内所有像素的平均值作为中心像素的值，得到的多视处理后的图像S_d[N,M]的表达式为：

$S_d\[i,j\]=\frac{1}{CD}\underset{i-(C-1)/2}{\overset{i+(C-1)/2}{\mathrm{\Σ}}}\underset{j-(D-1)/2}{\overset{j+(D-1)/2}{\mathrm{\Σ}}}{S}_{\mathrm{int}}\[i,j\]---\left(4\right)$

其中，i＝1，2，…N；j＝1，2，…M，对于边界点先扩展补零再处理；

(3)相位解缠

根据多视处理后S_d(N,M)的相位，得到解缠绕的干涉相位φ(N,M)；

(4)DEM生成

根据解缠后的干涉相位φ(N,M)，以及相应的目标高度反演公式(5)，得到每个目标点对应的高度H，将此时得到的高度记为H₀(N,M)，即输出目标的粗DEM；具体反演公式如下所示：

$H(N,M)=h+R_1(N,M)sin\[\mathrm{\α}+{\cos }^{-1}\left(\frac{\mathrm{\λ}\mathrm{\φ}(N,M)}{4\mathrm{\π}B}\right)\]---\left(5\right)$

其中，h为主天线高度，R₁(N,M)为目标到主天线斜距，α为基线水平方向的倾角，λ为雷达工作波长，φ(N,M)为解缠后的干涉相位，B为基线长度；

步骤二、将粗DEM即H₀(N,M)作为新的成像网格所在平面进行BP成像得到形变前的主辅图像，并进行干涉处理，得到新的DEM，进行多次迭代，得到精确的目标的DEM；

具体为：

将网格高度设为H₀(N,M)，重复步骤一，得到每个目标点对应的高度H₁(N,M)，然后将网格高度设为H₁(N,M)，重复步骤一，得到每个目标点对应的高度H₂(N,M)，以此类推，将网格高度设为H_n(N,M),重复步骤一，得到每个目标点对应的高度H_n+1(N,M)，设定迭代次数，判断n是否小于迭代次数，如果小于则继续迭代，否则，将最终迭代得到的H_n+1(N,M)设为精确的目标的DEM，记为H(N,M)；

步骤三、将迭代得到的目标的精确DEM即H(N,M)作为成像网格所在平面进行BP成像得到形变前后的主图，差分处理得到目标的形变信息△H(N,M)；

具体为：

第一步，BP成像，输入主天线形变前后距离压缩后的回波，输出形变前后主天线的SAR图像；

(1)网格划分

在已有主天线形变前回波的距离压缩结果S_r1(n,m)和形变后回波的距离压缩S_r1'(n,m)的情况下，对于已知的方位向分辨率ρ_a和距离向分辨率ρ_r，估计目标在方位向和地距向的范围，设置位于地平面方位向长度L和地距向长度W的BP网格，使网格覆盖目标范围，且网格高度为H(N,M)，网格点数设置为N×M，N为方位向网格点数，M为地距向网格点数，n为回波方位向采样点数，m为距离向采样点数；

(2)计算网格点斜距，进行相位补偿；

对于网格上的每个网格点A_ij，i＝1，2，…N，j＝1，2，…M，计算其形变前在每个方位时间t_p到雷达主天线的斜距R_p1及对应时延τ_p1，p＝0，1，…n，并根据方位时间t_p和时延计算其形变前在主天线回波距离压缩后的S_r1(n,m)的对应位置a_p1，b_p1，其中，a_p1为形变前A_ij在t_p时刻对应S_r1(n,m)中的方位向点数，b_p1为其对应的距离向点数，则在t_p时刻网格点A_ij形变前对应的主天线成像结果S_ap1[i,j]为S_r1[a_p1,b_p1]与补偿相位exp{-j2πf_cτ_p1}的乘积，同理计算其在形变后在每个方位时间t_p到主天线的斜距R_p1'及对应的时延τ_p1'，并根据方位时间t_p和时延计算其形变后在主天线回波距离压缩后的S_r1'(n,m)的对应位置a_p1'，b_p1'，其中，a_p1'为形变后A_ij在t_p时刻对应S_r1'(n,m)中的方位向点数，b_p1'为其对应的距离向点数，则在t_p时刻网格点A_ij形变后对应的主天线成像结果S_ap1'[i,j]为S_r1'[a_p1',b_p1']与补偿相位exp{-j2πf_cτ_p1'}的乘积，具体表示如(1)式，其中q＝1，1'，q＝1表示形变前，q＝1'表示形变后；

(3)相干累加

将上述步骤(2)中网格点A_ij所有方位时刻t_p时，p＝0，1，…n，相位补偿之后的形变前成像结果S_ap1[i,j]和形变后的成像结果S_ap1'[i,j]分别进行累加，即为该网格点的最后成像结果；对所有网格点执行以上操作，得到最终的形变前后主天线BP成像结果，即形变前的主SAR图像S_a1(N,M)和形变后的主SAR图像S_a1'(N,M)，如(2)式，其中q＝1，1'，q＝1表示形变前，q＝1'表示形变后；

第二步，差分干涉处理，输入形变前后主SAR图像，输出目标的形变信息△H(N,M)；

(1)生成差分干涉相位

直接将BP算法得到的形变前的主SAR图像S_a1(N,M)进行共轭得到S_a1 ^*(N,M)，并乘以形变后的主SAR图像S_a1'(N,M)，得到干涉图S_int(N,M)的相位即为差分干涉相位，即

S_int(N,M)＝S_a1 ^*(N,M)*S_a1'(N,M)(6)

(2)多视处理

多视处理即对干涉图S_int(N,M)进行均值滤波；首先选取C₁×D₁的图像窗口，则多视系数为C₁×D₁，以该窗口内所有像素的平均值作为中心像素的值，表达式为：

$S_d\[i,j\]=\frac{1}{C_1{D}_1}\underset{i-(C_1-1)/2}{\overset{i+(C_1-1)/2}{\mathrm{\Σ}}}\underset{j-(D_1-1)/2}{\overset{j+(D_1-1)/2}{\mathrm{\Σ}}}{S}_{\mathrm{int}}\[i,j\]---\left(7\right)$

其中，i＝1，2，…N；j＝1，2，…M，设S_d(N,M)的相位即是差分干涉相位△φ(N,M)；

(3)形变反演

输入上一步得到的差分干涉相位△φ(N,M)，再根据相应的目标形变反演公式(8)，得到每个目标点对应的形变△H(N,M)，具体反演公式如下所示：

$\mathrm{\Δ}H(N,M)=H(N,M)-h-\sqrt{{(\sqrt{{\left(H\right(N,M)-h)}^2+{\left(r_p\right(N,M)-r)}^2}-\frac{\mathrm{\λ}\mathrm{\Δ}\mathrm{\φ}(N,M)}{4\mathrm{\π}})}^2-{\left(r_p\right(N,M)-r)}^2}---\left(8\right)$

其中，H(N,M)为目标高度，即步骤二中多次迭代得到的目标DEM，h为主天线高度，r_p为目标在地面的投影点距原点的长度，r为主天线的地面投影距原点长度，λ为雷达工作波长，△φ(N,M)为差分相位。