[發明專利]一種基于曲面后向投影算法的形變反演方法

有效

• 申請號: 201510570376.0
• 申請日: 2015-09-09
• 公開（公告）號: CN105182337B
• 公開（公告）日: 2018-04-17
• 發明（設計）人: 孫兵;聶琦;張施雨;張潔瓊
• 申請（專利權）人: 北京航空航天大學
• 主分類號: G01S13/90
• 分類號: G01S13/90
• 代理公司: 北京永創新實專利事務所11121
• 代理人: 趙文穎

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)；即有：

Saq[i,j]=ΣpSrq[apq,bpq]·exp{-j2πfcτpq}---(2)]]>

其中，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]的表達式為：

Sd[i,j]=1CDΣi-(C-1)/2i+(C-1)/2Σj-(D-1)/2j+(D-1)/2Sint[i,j]---(4)]]>

其中，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+R1(N,M)sin[α+cos-1(λφ(N,M)4πB)]---(5)]]>

其中，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₁，以該窗口內所有像素的平均值作為中心像素的值，表達式為：

Sd[i,j]=1C1D1Σi-(C1-1)/2i+(C1-1)/2Σj-(D1-1)/2j+(D1-1)/2Sint[i,j]---(7)]]>

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

(3)形變反演

輸入上一步得到的差分干涉相位△φ(N,M)，再根據相應的目標形變反演公式(8)，得到每個目標點對應的形變△H(N,M)，具體反演公式如下所示：

ΔH(N,M)=H(N,M)-h-((H(N,M)-h)2+(rp(N,M)-r)2-λΔφ(N,M)4π)2-(rp(N,M)-r)2---(8)]]>

其中，H(N,M)為目標高度，即步驟二中多次迭代得到的目標DEM，h為主天線高度，r_p為目標在地面的投影點距原點的長度，r為主天線的地面投影距原點長度，λ為雷達工作波長，△φ(N,M)為差分相位。