[發明專利]一種迭代曲面預測InSAR成像及高度估計方法

有效

• 申請號: 201410442793.2
• 申請日: 2014-09-02
• 公開（公告）號: CN104698457B
• 公開（公告）日: 2017-02-15
• 發明（設計）人: 韋順軍;張曉玲;鄒光浩;熊海進
• 申請（專利權）人: 電子科技大學
• 主分類號: G01S13/90
• 分類號: G01S13/90
• 代理公司: 電子科技大學專利中心51203
• 代理人: 曾磊

1.一種迭代曲面預測InSAR成像及高度估計方法，其特征是它包括以下幾個步驟：

步驟1、初始化InSAR成像系統參數：

初始化InSAR成像系統參數，包括：雷達載波波長，記為λ，雷達平臺主天線發射信號帶寬，記為B，雷達發射脈沖時寬，記為T_r，雷達采樣頻率，記為F_s，主天線雷達入射角，記為θ，兩組天線之間的基線長，記為B_l，基線與水平方向夾角，記為α，雷達脈沖重復頻率，記為PRF，平臺運動速度矢量，記為V_r，雷達系統距離向采樣點數，記為N_r，雷達系統方位向采樣點數，記為N_a，雷達平臺主天線初始位置，記為P_M(0)，雷達平臺副天線初始位置，記為P_S(0)；上述參數中，雷達系統工作的信號波長λ、雷達平臺主天線發射的信號帶寬B、雷達平臺主天線發射的脈沖時寬T_r、雷達平臺接收系統的采樣頻率F_s、主天線雷達入射角θ、兩組天線之間的基線長B_l、基線與水平方向夾角α，雷達脈沖重復頻率PRF，在InSAR雷達系統設計過程中已經確定；平臺運動速度矢量V_r，雷達系統距離向采樣點數N_r、雷達系統方位向采樣點數N_a、雷達平臺主天線初始位置P_M(0)、雷達平臺副天線初始位置P_S(0)在InSAR雷達成像觀測方案設計中已經確定；

步驟2、獲取InSAR原始回波數據：

InSAR主天線在距離向第t個快時刻方位向第k個慢時刻的原始回波數據，記為E_m(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a，其中t和k為自然數，t表示距離向第t個快時刻，k表示方位向第k個慢時刻，N_r為步驟1初始化得到的雷達系統距離向采樣點數，N_a為步驟1初始化得到的雷達系統方位向采樣點數；InSAR副天線在距離向第t個快時刻方位向第k個慢時刻的原始回波數據，記為E_s(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a；在InSAR實際成像中，InSAR主天線在距離向第t個快時刻方位向第k個慢時刻的原始回波數據E_m(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a，和InSAR副天線在距離向第t個快時刻方位向第k個慢時刻的原始回波數據E_s(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a，由InSAR系統數據接收機提供；

步驟3、對InSAR原始回波數據進行距離壓縮：

采用傳統標準合成孔徑雷達標準距離壓縮方法對步驟2中得到的InSAR主天線在距離向第t個快時刻方位向第k個慢時刻的原始回波數據E_m(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a，進行距離壓縮得到InSAR主天線在距離向第t個快時刻方位向第k個慢時刻的距離壓縮后數據，記為S_m(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a，其中t表示距離向第t個快時刻，k表示方位向第k個慢時刻；采用傳統標準合成孔徑雷達標準距離壓縮方法對步驟2中得到的InSAR副天線在距離向第t個快時刻方位向第k個慢時刻的原始回波數據E_s(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a，進行距離壓縮得到InSAR副天線在距離向第t個快時刻方位向第k個慢時刻的距離壓縮后數據，記為S_s(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a；

步驟4、初始化雷達投影成像曲面參數：

InSAR雷達投影成像曲面選擇為地表水平面坐標系，該坐標系水平橫軸記為X軸，該坐標系水平縱軸記為Y軸，雷達投影成像曲面的中心坐標位于[0,0,0]，雷達投影成像曲面的X軸向分辨單元數，記為N_x，雷達投影成像曲面的Y軸向分辨單元數，記為N_y，雷達投影成像曲面的X軸向成像范圍，記為W_x，雷達投影成像曲面的Y軸向成像范圍，記為W_y，將投影成像曲面均勻等間隔進行劃分，投影成像曲面中第a個Y軸向和第r個X軸向分辨單元的位置矢量，記為P_T(a,r)＝[x(a,r),y(a,r),z(a,r)]，a＝1,…,N_x，r＝1,…,N_y，其中a和r均為自然數，a表示投影成像曲面的第a個X軸向分辨單元，r表示投影成像曲面的第r個Y軸向分辨單元，x(a,r)、y(a,r)和z(a,r)分別表示投影成像曲面中第a個Y軸向第r個X軸向分辨單元的X軸向位置、Y軸向位置和高度；

步驟5、初始化迭代曲面預測InSAR成像及高度估計方法所需的參數：

初始化迭代曲面預測InSAR成像及高度估計方法的參數，包括：成像處理迭代的最大次數，記為MI，成像處理迭代閾值，記為ε，第i次迭代過程中雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元的高度，記為Z_i(a,r)，a＝1,…,N_x，r＝1,…,N_y，i＝0,…,MI，其中i為自然數，表示成像算法的第i次迭代，N_x為步驟4初始化得到的雷達投影成像曲面的X軸向分辨單元數，N_y為步驟4中初始化得到的雷達投影成像曲面的Y軸向分辨單元數；Z₀(a,r)為第0次迭代過程中雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元的高度，即該步驟中初始化得到雷達投影成像曲面高度，初始化雷達投影成像曲面高度Z₀(a,r)的值為Z₀(a,r)＝0，a＝1,…,N_x，r＝1,…,N_y；

步驟6、計算雷達投影成像曲面中分辨單元到平臺的斜距：

在成像處理第i次迭代中，若i＝1，則更新投影成像曲面中第a個Y軸向第r個X軸向分辨單元的高度為z(a,r)＝Z₀(a,r)，更新得到投影成像曲面中第a個Y軸向和第r個X軸向分辨單元的位置矢量P_T(a,r)＝[x(a,r),y(a,r),Z0(a,r)]；若i＞1，則更新投影成像曲面中第a個Y軸向第r個X軸向分辨單元的高度為z(a,r)＝Z_i-1(a,r)，得到P_T(a,r)＝[x(a,r),y(a,r),Zi(a,r)]，其中Z₀(a,r)為步驟5中初始化得到的第0次迭代過程中雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元的高度Z₀(a,r)＝0，a＝1,…,N_x，r＝1,…,N_y，而Z_i-1(a,r)為成像處理第i-1次迭代過程中雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元的高度值；

采用公式P_M(k)＝P_M(0)+V_r·k/PRF，k＝1,2,…,N_a，計算得到雷達平臺主天線第k個慢時刻的位置，記為P_M(k)，k＝1,2,…,N_a，其中P_M(0)為步驟1初始化得到的雷達平臺主天線初始位置，V_r為步驟1初始化得到的平臺運動速度矢量，PRF為步驟1初始化得到的雷達脈沖重復頻率，k表示方位向第k個慢時刻，N_a為步驟1初始化得到的雷達系統方位向采樣點數；采用公式P_S(k)＝P_S(0)+V_r·k/PRF，k＝1,2,…,N_a，計算得到雷達平臺副天線第k個慢時刻的位置，記為P_S(k)，k＝1,2,…,N_a，其中P_S(0)為步驟1初始化得到的雷達平臺副天線初始位置；

采用公式R_m(k；a,r)＝||P_M(k)-P_T(a,r)||₂，a＝1,…,N_x，r＝1,…,N_y，k＝1,2,…,N_a，計算得到雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元到InSAR主天線的斜距，記為R_m(k；a,r)，a＝1,…,N_x，r＝1,…,N_y，k＝1,2,…,N_a，其中||·||₂表示向量2范數運算符號，P_T(a,r)為步驟4得到的投影成像曲面中第a個Y軸向和第r個X軸向分辨單元的位置矢量；采用公式R_s(k；a,r)＝||P_S(k)-P_T(a,r)||₂，a＝1,…,N_x，r＝1,…,N_y，k＝1,2,…,N_a，計算得到雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元到InSAR副天線的斜距，記為R_s(k；a,r)，a＝1,…,N_x，r＝1,…,N_y，k＝1,2,…,N_a；

步驟7、采用標準后向投影成像算法進行成像：

對步驟3得到的InSAR主天線在距離向第t個快時刻方位向第k個慢時刻的距離壓縮后數據S_m(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a，以及步驟6得到的雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元到InSAR主天線的斜距R_m(k；a,r)，a＝1,…,N_x，r＝1,…,N_y，采用傳統標準后向投影成像算法進行成像，得到雷達平臺主天線數據在投影成像曲面的成像結果，記為B_m(a,r)，a＝1,…,N_x，r＝1,…,N_y；對步驟3得到的InSAR副天線在距離向第t個快時刻方位向第k個慢時刻的距離壓縮后數據S_s(t,k)，t＝1,2,…,N_r，k＝1,2,…,N_a，以及步驟6得到的雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元到InSAR副天線的斜距R_s(k；a,r)，a＝1,…,N_x，r＝1,…,N_y，采用傳統標準后向投影成像算法進行成像，得到雷達平臺副天線數據在投影成像曲面的成像結果，記為B_s(a,r)，a＝1,…,N_x，r＝1,…,N_y；

步驟8、InSAR主副天線干涉相位提?。?!-- SIPO -->

采用公式φ(a,r)＝∠((B_m(a,r))^*·B_s(a,r))，a＝1,…,N_x，r＝1,…,N_y，計算得到投影成像曲面下InSAR主副天線成像結果的干涉相位，記為φ(a,r)，a＝1,…,N_x，r＝1,…,N_y，其中B_m(a,r)和B_s(a,r)為步驟7中得到的雷達平臺主天線與副天線數據在投影成像曲面的成像結果，∠(·)表示求解相位運算符號，上標^*表示求解共軛運算符號；

步驟9、InSAR干涉相位解纏：

采用傳統標準干涉相位解纏方法對步驟8中得到的投影成像曲面下InSAR主副天線成像結果的干涉相位φ(a,r)，a＝1,…,N_x，r＝1,…,N_y，進行干涉相位解纏處理，得到解纏后的InSAR干涉相位，記為，a＝1,…,N_x，r＝1,…,N_y，其中a表示投影成像曲面第a個X軸向分辨單元，r表示投影成像曲面第r個Y軸向分辨單元；

步驟10、計算雷達投影成像曲面條件下的地形高度：

采用公式計算得到投影成像曲面第a個Y軸向第r個X軸向分辨單元P(a,r)的地形高度，記為H(a,r)，a＝1,…,N_x，r＝1,…,N_y，其中a表示投影成像曲面第a個X軸向分辨單元，r表示投影成像曲面第r個Y軸向分辨單元，λ為步驟1中初始化得到的雷達載波波長，為步驟6中時得到的雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元到InSAR主天線的斜距，ceil(·)表示四舍五入運算符號，sin(·)表示正弦函數運算符號，θ為步驟1中初始化得到的主天線雷達入射角，B_l為步驟1中初始化得到的兩組天線之間的基線長，cos(·)表示余弦函數運算符號，α為步驟1中初始化得到的基線與水平方向夾角，為步驟7得到的解纏后的InSAR干涉相位，π為圓周率；

步驟11、判定算法迭代條件：

采用公式計算得到雷達投影成像曲面條件下地形高度絕對值的平均值，記為H_mean；

采用公式Z_i(a,r)＝H(a,r)+Z_i-1(a,r)，a＝1,…,N_x，r＝1,…,N_y，將地形高度H(a,r)，a＝1,…,N_x，r＝1,…,N_y，補償到成像處理中第i次迭代過程中雷達投影成像曲面高度Z_i(a,r)，a＝1,…,N_x，r＝1,…，N_y；

如果平均值H_mean＞ε且i≤MI，將成像處理中迭代次數i加1得到i←i+1，再重復執行步驟6至步驟11；否則終止成像處理迭代，Z_i(a,r)即為雷達投影成像曲面的最終地形高程，其中i表示成像處理第i次迭代，ML為步驟4中初始化得到算法迭代的最大次數，Σ(·)表示向量求和運算符號，|·|表示求絕對值運算符號，←表示賦值運算符號，ε為步驟4中初始化得到的成像處理迭代閾值，H(a,r)，a＝1,…,N_x，r＝1,…,N_y，為步驟10中得到的投影成像曲面第a個Y軸向第r個X軸向分辨單元P(a,r)的地形高度，Z_i-1(a,r)為成像處理第i-1次迭代過程中雷達投影成像曲面中第a個Y軸向第r個X軸向分辨單元的高度值。