利用非線性矢量曲面構(gòu)建InSAR相位圖像模型的方法
【專利摘要】本發(fā)明涉及雷達、光學干涉圖像處理，尤其涉及合成孔徑雷達干涉測量(InSAR)相位圖像處理。首先對離散InSAR相位圖像進行多層次小波分解，將相位表示為低頻部分和高頻部分和的形式；利用估計窗口內(nèi)相位點高頻分量的最大值確定主矢量的方向，根據(jù)信號和殘差點的模極值變化不同，對估計窗口內(nèi)相位進行平滑；在估計窗口內(nèi)對各分解方向求取相位矢量和，構(gòu)建非線性相位矢量曲面模型；利用非線性相位矢量曲面顧及測區(qū)地形因素起伏變化的非線性部分，得到的相位模型能預估測區(qū)相位大小變化趨勢，幫助相位濾波和解纏，提高處理結(jié)果的精度。
【專利說明】
利用非線性矢量曲面構(gòu)建I nSAR相位圖像模型的方法
技術(shù)領(lǐng)域
[0001]本發(fā)明涉及雷達、光學干涉圖像處理，尤其涉及合成孔徑雷達干涉測量（InSAR)相 位圖像處理。
【背景技術(shù)】
[0002] InSAR技術(shù)利用側(cè)視雷達回波信號所攜帶的相位信息來計算地面的三維數(shù)據(jù)。從 最初的雷達數(shù)據(jù)格式SLC(單視復圖像)得到InSAR相位圖像，再得到最終高程數(shù)據(jù)產(chǎn)品，需 要經(jīng)過兩幅圖像配準、復共輒相乘、相位濾波、相位解纏和基線估計等數(shù)據(jù)處理步驟才能得 到數(shù)字高程模型。
[0003] InSAR相位圖像是SAR數(shù)據(jù)處理的中間產(chǎn)品，也是圖像濾波和相位解纏的基礎(chǔ)數(shù) 據(jù)，由于后續(xù)的步驟需要轉(zhuǎn)換成數(shù)學模型進行處理，構(gòu)建精確、有效的InSAR圖像的相位模 型是重點要解決的技術(shù)難題之一。
[0004] InSAR干涉圖像可表示為二維離散相位信號形式，干涉相位Φ可寫成線性項 Φ linear和非線性項Φ_-linear和的形式，如公式（1-1):
(1.-1)
[0006] 其中，Dlinear、D_-linear分別表示數(shù)據(jù)處理系統(tǒng)線性和非線性組成部分;v表示線性 地表起伏變化的速率，A h表示地表起伏的高度變化，T表示獲取干涉對的雷達數(shù)據(jù)時間基 線。公式（1-1)中等號右邊的第一項主要是垂直基線、時間基線構(gòu)成的線性關(guān)系；第二項主 要是由數(shù)據(jù)處理系統(tǒng)和地形欠采樣引起相位的非線性變化以及相位殘余構(gòu)成，非線性部分 造成傳統(tǒng)枝切線回路無法閉合以及濾波和解纏步驟誤差傳遞等現(xiàn)象，直接用線性的方法無 法獲取高精度相位，因此利用相位模型必須要顧及非線性變化的相位，即構(gòu)建非線性相位 模型。
[0007] 目前，InSAR相位模型可通過線性相位估值的方法獲取，設(shè)干涉圖的線性相位模型 為U，那么：
(1-2)
[0009]其中，尤，/;為通過二維Chirp-Z變換得到的二維空間頻率分量，X?對應于距離 向頻譜偏移量Α ω
?其中Fs為系統(tǒng)的距離向采樣率，C(x，y)為被截取局部干涉圖。 詳見文南犬：Bamler R，Adaln N，et al .Noise-Induced Slope Distortion in 2_D Phase Unwrapping by Linear Estimators with Application to SAR Interferometry，1998， 中公開的技術(shù)方案。從（1-2)式可以看出線性相位模型可以用一階近似相位平面和一個復 常數(shù)部分構(gòu)成。算法求解過程中采用二維離散傅里葉變換(FFT)方法，通過搜索窗口內(nèi)二維 信號頻譜的峰值位置來估計干涉條紋頻率，取估計窗口的中心相位頻率作為信號的相位頻 率估計，通過最大化估計窗口相位頻率的似然函數(shù)，可以估計出窗口中心(m, η)處的真實局 部空間頻率fx和fy。經(jīng)過二維離散傅立葉變換(DFT)后，得到頻譜峰值處的頻率能夠近似局 部頻率，詳見文獻:朱岱寅.利用線性相位模型提高干涉SAR圖像相干性及改進相干性估計， 2005,中公開的技術(shù)方案。
[0010]這是相位模型的線性構(gòu)造方法，模型假設(shè)估計窗口內(nèi)二維相位滿足線性變化的特 征，然而在實際數(shù)據(jù)處理中，線性相位模型仍然存在以下缺陷：
[0011] (1)相位殘差點反映在頻率域中是遠離估計頻率的干擾信號，此時使用估計窗口 的方法估計的相位不再滿足線性模型假設(shè)，容易導致線性相位估計模型出現(xiàn)偏差。
[0012] (2)線性模型平滑了相位中的非線性變化部分，使用頻譜峰值的最大化似然函數(shù) 使得窗口內(nèi)相位樣本局部統(tǒng)一，模糊了相位變化的細節(jié)。
[0013] (3)線性模型的特點導致估計窗口小且固定，計算時間長，不能有效抑制殘差點導 致低相干區(qū)域出現(xiàn)較大的模型誤差。
[0014] 目前，對于非線性相位使用二維離散傅里葉變換(FFT)方法，相位估計窗口可以根 據(jù)均值的大小自適應變化，反映了部分地形變化的輪廓，詳見文獻:Wang Q S，Huang H F， An X Y,et al. Improving phase Unwrapping Techniques by the Use of Nonlinear Phase Model，2011，中公開的技術(shù)方案。但FFT方法仍不能顧及地形和相位變化的細節(jié)，無 法依據(jù)地形變化的視線向、距離向和方位向進行區(qū)分，而且使用估計窗口的均值估計相位 不準確，也容易造成誤差的傳遞。

【發(fā)明內(nèi)容】

[0015] 本發(fā)明對離散InSAR相位圖像進行多層次小波分解，將相位表示為低頻部分和高 頻部分和的形式;利用估計窗口內(nèi)相位點高頻分量的最大值確定主矢量的方向，根據(jù)信號 和殘差點的模極值變化不同，對估計窗口內(nèi)相位進行平滑;在估計窗口內(nèi)對各分解方向求 取相位矢量和，構(gòu)建非線性相位矢量曲面模型；利用非線性相位矢量曲面顧及測區(qū)地形因 素起伏變化的非線性部分，得到的相位模型能預估測區(qū)相位大小變化趨勢，幫助相位濾波 和解纏，提高處理結(jié)果的精度。
[0016]具體方案流程為：
[0017] 1.1將相位干涉圖以復數(shù)的浮點形式讀取到Matlab中，標記矩陣行、列大小為M、N， 單位為像素。將干涉相位圖劃分大小為kXl的初始估計窗口，并計算值得信任的第一個估 計窗口，將窗口內(nèi)相位點(m，n)單頻信號形式轉(zhuǎn)換為三角函數(shù)弦分量的復數(shù)相位形式：Φ (m，n)=angle(exp[ j(fxm+fyn)]) = cos(fxm+fyn)+j sin(fxm+fyn)，其中，fx和fy分別為相位 點行向量和列向量的頻率；
[0018] 1.2在估計窗口內(nèi)，利用離散小波變換對干涉圖像正弦分量和余弦分量進行多層 分解，將相位表示為低頻部分和高頻部分和的形式；
[0019] 1.3在估計窗口內(nèi)，利用相位點高頻分量的最大值確定主矢量的方向，也即是下一 步搜索窗口的方向；
[0020] 1.4在估計窗口內(nèi)，計算相位主矢量，利用高頻部分主矢量平滑窗口內(nèi)的相位頻 率；
[0021 ] 1.5在估計窗口內(nèi)，根據(jù)各層分解分量的高頻矢量和構(gòu)建顧及相位細節(jié)變化的非 線性相位曲面模型，將相位表示成非線性相位和線性相位和的形式Φ = Φ line3ar + J _ _ _ Φ ii_r，對非線性部分利用(爪，》) = arg(Z< )= m，g(l + 6十表;）對j層分解的三 Η 個高頻方向求矢量和，其中j = 1，2···J，^.、分別表示對j層高頻的水平、垂直和視 線方向分量求矢量和，矢量求和保證了非線性曲面的特性；
[0022] 1.6對小波分解相位進行逆變換得到相位估計模型，在估計窗口 k X 1內(nèi)，對相位點 也,1 (m，η)，利用(電，《) = a取(紅_1 = arg (I + ΣI j可得到估計窗口內(nèi)非線 性相位估計模型，其中，WTl ·)為小波逆變換，(^為小波分解低頻部分，ΣΙ為小波分解 高頻部分和。
[0023] 本發(fā)明技術(shù)方案帶來的有益效果
[0024] (1)利用構(gòu)造的非線性相位模型顧及了非線性相位的偏差，實現(xiàn)了相位最優(yōu)估計；
[0025] (2)相位濾波是去除干涉圖中的噪聲，而應區(qū)分地形因素引起的相位相干斑噪聲； 相位解纏是將纏繞相位（以2π為模)展開成實際相位的過程，但由于受到地形因素的影響使 得解纏路徑有環(huán)路，不能正確解纏，若使用本發(fā)明估計相位模型來顧及測區(qū)地形因素，能預 估測區(qū)相位變化趨勢，幫助相位濾波和解纏，提高處理結(jié)果的精度。
[0026] (3)非線性相位模型不再以局部線性相位頻率估計代替局部相位，而是以反映地 形因素的非線性曲面來近似估計窗口相位。
【附圖說明】
[0027] 圖1為以相位垂直方向為主方向的相位矢量求和圖；
[0028] 圖2為以相位水平方向為主方向的相位矢量求和圖；
[0029] 圖3為以相位視線方向為主方向的相位矢量求和圖；
[0030] 圖4為非線性矢量曲面構(gòu)建InSAR相位圖像流程圖。
[0031 ]圖例說明，1-側(cè)視雷達視線向，2-水平方向，3-垂直方向，V -側(cè)視雷達視線向反方 向，2'-水平方向反方向，3'-垂直方向反方向。
【具體實施方式】
[0032]下面結(jié)合附圖對本發(fā)明技術(shù)方案進一步說明。
[0033]由于星載側(cè)視雷達以相同頻率發(fā)射窄帶脈沖，在斜距相同的條件下，不同地形條 件反射的雷達信號不同，反映在干涉圖中就是或密集或稀疏的干涉條紋，那么相位條紋的 疏密程度也反映了干涉圖中地形的特征，可以根據(jù)相位條紋頻率建立相位模型。
[0034]為了方便本發(fā)明中技術(shù)描述，首先進行了以下定義：
[0035] 定義1:采樣點地形坡度的相位差導數(shù)表示。
[0036] 設(shè)參數(shù)τ為地面采樣點的地形坡度，τ與相位信號系統(tǒng)各參數(shù)之間的關(guān)系為：
P)
[0038]其中，c為雷達波傳播速度;B丄為基線在垂直于視線方向的投影，即有效基線長度； λ為波長;R為斜距；Θ為入射角;?·φ為干涉圖相位差頻率。根據(jù)公式(2)可推出地面傾角與入 射角之間的關(guān)系公式(3):
(3)
[0040]由上式可知，可以根據(jù)距離向和方位向上的局部頻率Af來估計干涉圖相位模型。 設(shè)Δ fx、△ fy分別為距離向、方位向相位差頻率，可得下式(4):
[0042] 聯(lián)立相位差的公式： 可得：
u
[0046]本方法利用離散小波分層分解，獲取非線性相位曲面模型來逼近估計窗口內(nèi)非線 性相位，并對曲面模型進行逆變換得到非線性相位模型。具體步驟如下：
[0047]步驟1:準備基礎(chǔ)數(shù)據(jù)的工作。
[0048] 獲取兩幅SAR圖像51與52配準、濾波后的干涉圖像，包括附加的系統(tǒng)參數(shù)信息（波 長、基線、下視角、斜距和如果可以獲取的地形DEM信息），計算干涉圖像相干系數(shù)γ和離散 干涉相位Φ。在實際圖像處理過程中，
，式中1表示多視 數(shù)，取值為1:5或2:10，表示SAR圖像共輒相乘。
[0049] 步驟2:讀取干涉圖得到相位矩陣并轉(zhuǎn)換成相位頻率的形式。
[0050] 將相位干涉圖以復數(shù)的浮點形式讀取到Matlab中，設(shè)矩陣行、列大小為Μ、Ν，單位 為像素。將干涉相位圖劃分大小為ΚXL的初始估計窗口，并將窗口內(nèi)相位點(m，η)的單頻信 號形式轉(zhuǎn)換為三角函數(shù)弦分量的復數(shù)相位形式： i,kw,n、=-"ngle(CKp j{ f'iii + f'n] \
[0051 ] / 1 L , J/ ?? =cos ( /> + /;") + / sin (./:/" + />)
[0052] 其中，f4Pfy分別為相位點行向量和列向量的頻率。
[0053] 另外，使用相位導數(shù)也就是相位差的形式得到相位頻率，用下面的公式進行計算，
[0054] L = 〇Φ/^ = φη1ιυι + ~2φηιη
[0055] f, = 9φ/δν = 4,"+ι + ~ 2Φη,."
[0056] 上式中φω,η表示干涉圖像ΜΧΝ中相位點(m，n)的相位值;χ和y為矩陣行向量和列 向量。
[0059] 在實際離散相位圖像算法中，根據(jù)第1步的相位相干系數(shù)搜索干涉圖像，計算值得 信任的第一個估計窗口，可使用下面的方法:首先搜索干涉圖中局部估計窗口的平均最大 相干系數(shù)若這樣的窗口數(shù)量不止一個，設(shè)數(shù)量為N，則根據(jù)地形坡度角與相位頻率 V kxk J 公式(3 )，求取窗口內(nèi)最小地形坡度窗口設(shè)數(shù)量為Μ，確定為初始估計窗口；否則，若f > 1， 則將窗口下標值最小的窗口確定為第一個估計窗口，從第一個估計窗口開始進行相位模型 的計算。
[0060] 步驟3:利用離散小波變換的方法對相位頻率進行分解。
[0061]在估計窗口內(nèi)，利用離散小波變換對干涉圖像正弦分量和余弦分量進行多層分 解，將相位表示為低頻部分和高頻部分和的形式。對估計相位Φ在矢量域求和，分別用不同 頻率的矢量形式表示(根據(jù)頻率變化的特點，使用小波變換后的不同分辨率表示，低頻部分 表示非線性相位;高頻部分為相位的細節(jié)信息，地形的起伏變化反映在頻率上就是頻率域 中的高頻部分）。
[0062]相位Φ k, 1 (m，η)的小波頻譜變換為：
[0064] 其中，Φ 〇為窗口中心相位點，成=訓沙(為(L/)exp(-./2；r(./;/c +//))); Δ φ邊示 Φ 〇的微小變化區(qū)域;J為相位分解的層次J = 1，2，…，n;Cj和Dh，l，r分別表亦小波分解中的低 頻和高頻部分，其中Dh，l，r表示各個分辨率層的相位細節(jié)信息，H，L，R表示窗口相位搜索的三 個方向（水平、垂直、視線)或其反方向，詳細的相位搜索方向可參見附圖1、2。
[0065]在實際離散相位圖像算法中，對于相位窗口的第t次分解可以由第t-Ι次分解得 到，t-Ι次的分解窗口結(jié)構(gòu)為：
[0067] X表示估計窗口中的相位頻率。
[0068] 步驟4:利用分層分解結(jié)果獲取高頻部分主矢量的方向，確定估計窗口滑動的方 向。
[0069] 根據(jù)公式(8)對估計窗口內(nèi)小波分解的結(jié)果，得到窗口內(nèi)相位點的頻率表示形式：
[0071] 其中，Φ?ρ為小波分解低頻部分K為小波分解高頻部分;Pi為取值范圍為(〇，1) 的加權(quán)系數(shù);在矢量加權(quán)求和中，系數(shù)Pi的大小決定了第i個信號分量在估計窗口中的權(quán) 重。
[0072] 在估計窗口內(nèi)，利用相位點高頻分量的最大值確定主矢量的方向，也即是下一步 搜索窗口的方向。設(shè)當前窗口 j層的主矢量方向為M，標記為^^，是根據(jù)j-Ι層的主矢量方向 確定的：
[0075] 上式中表示H、L、R三個方向頻率對主矢量的角度偏差，也就是三分量在 主矢量方向的矢量分解。
[0076] 在實際離散相位圖像算法中，對矢量分解的角度偏差通過如下公式進行計算：
[0080]在上式中，若視線向為主矢量方向，水平、垂直方向矢量的角度偏差用A 01/31獲 取;若水平、垂直為主矢量方向，視線向矢量的角度偏差用Α Θ2Λτ、△ θ3/3?獲取。
[0081 ]步驟5:在估計窗口內(nèi)利用高頻部分主矢量平滑窗口內(nèi)的相位頻率。
[0082]由于窗口內(nèi)的相位噪聲表現(xiàn)為信號無關(guān)性，若當前估計窗口判斷為相位噪聲的相 位點頻率值大于窗口內(nèi)相位頻率值的和，則認定為噪聲信號，其非線性相位為無效估計，利 用已經(jīng)構(gòu)建相位模型的相鄰窗口頻率均值，與當前窗口高頻部分的數(shù)值之和作為當前窗口 .A'xL 的相位模型，若物?，《)>，貝1J: ?=1
[0084] 根據(jù)干涉相位頻率提取相位矢量，可以設(shè)置閾值(窗口頻率變換主瓣值的均值）， 當估計窗口內(nèi)提取得到的頻率值大于該門限時，由于窗口中頻譜變化最大值表示了相位變 化的主方向，發(fā)明中使用窗口內(nèi)信號頻譜最大值作為矢量閾值：Ize
[0085] 對于高頻部分的處理為：= Σ U.A. +/d- + /"A) 'H，L，R表示三個方向（水 平、垂直、視線）的正方向高頻gz和低頻hz，在方向搜索可時同時顧及這三個方向的反方向。 若此時根據(jù)要控制的方向個數(shù)(2個、4個、6個或8個)對各方向求矢量和。
[0086] 在實際離散相位圖像算法中，相位高頻部分通過如下公式進行計算：
[0087] Σ (盡--'?--十 W)= Σ (2_, {^ !D,, (in,n)-D, \ym-[k-l)j2,n-(k-l)l2^Cj (>?,?)] j = zeH,L,R z.eH,L,R + (2-7 (2-人1?_ (辦，π)-為(m -(左-1)/2，'卜(免" 1)/2)) (67)(??？7)) '2~i(2-iDz(nKn)-Dz(m-(k~l)/2^i-(k-l)/2)f Η- |2-?λ (?；, η)- (m-(k-1)/2 .η-(k-I )/2))
[0088] 上式中(m,n)為估計窗 口中相位點，02(111-(1^-1)/2，11-(1^-1)/2)表示(111，11)鄰域的 相位點。
[0089] (k-Ι )/2保證相位變換點在估計窗口范圍內(nèi)。
[0090] 步驟6:在估計窗口內(nèi)利用高頻部分矢量和構(gòu)造非線性相位曲面模型。
[0091] 由公式（2) Φ = Φ η_γ+Φ non-linear ? 根據(jù)各層分解分量高頻矢量和構(gòu)造顧及相位 細節(jié)變化的構(gòu)建非線性相位曲面模型。
[0092] 設(shè)式（8)中
為信號Φ:分解的i層主矢量， 那么主矢量和的相位反映非線性相位巾_得到公式(10): J
[0093] 為―·(衍，= ) (1〇)
[0094] 同時對局部窗口相位高頻矢量求和保證了非線性曲面的特性，如公式(11):
[0096] 其中，上式中右邊g、慰、g表示對j層高頻的水平、垂直和視線方向分量求矢 量和。砹1表示視線向為主矢量方向，水平、垂直方向矢量和用(/? 1獲取;若水 平、垂直為主矢量方向，視線向矢量的矢量和分別用、GfA爲/;τ獲取。
[0097] 步驟7:對小波分解相位進行逆變換，可求取非線性相位估計模型。
[0098] 在估計窗口（k，l)內(nèi)可以得到非線性相位： = arg^T'1 {fm + ./>?})
[0099] f r Λ =ill'S Σ+ AnM," {m-fl)i-'xP(-/2^(./> + A11)) Π2) v k,i ' ' J
[0100] 其中，胃!^1!： ·)為小波逆變換，<Kp為小波分解低頻部分為小波分解高頻部 分在實際離散相位圖像算法中，對于估計窗口的第t次逆變換由第t+i次分解得到，t+i次的 逆變換窗口結(jié)構(gòu)為：
[0102] X-1表示Cj或Dh,l,r的逆變換。
[0103] 步驟8:判斷是否搜索完整個圖像，是結(jié)束流程;否則，窗口向相位極值方向滑動。
[0104] 以上借助具體實施例對本發(fā)明做了進一步描述，但是應該理解的是，這里具體的 描述，不應理解為對本發(fā)明的實質(zhì)和范圍的限定，本領(lǐng)域內(nèi)的普通技術(shù)人員在閱讀本說明 書后對上述實施例做出的各種修改，都屬于本發(fā)明所保護的范圍。
【主權(quán)項】
1.利用非線性矢量曲面構(gòu)建InSAR相位圖像模型的方法，其特征在于，包括w下步驟： 步驟一:將相位干設(shè)圖W復數(shù)的浮點形式讀取到Matlab中，標記矩陣行、列大小為M、N， 單位為像素;將干設(shè)相位圖劃分大小為kXl的初始估計窗口，并計算值得信任的第一個估 計窗口，將窗口內(nèi)相位點(m，n)單頻信號形式轉(zhuǎn)換為Ξ角函數(shù)弦分量的復數(shù)相位形式：Φ (m,n) =angle(e邱[j(fxm+fyn)]) =cos(fxm+fyn)+jsin(fxm+fyn)，其中，fx和fy分別為相位 點行向量和列向量的頻率； 步驟二:在估計窗口內(nèi)，利用離散小波變換對干設(shè)圖像正弦分量和余弦分量進行多層 分解，將相位表示為低頻部分和高頻部分和的形式； 步驟Ξ:在估計窗口內(nèi)，利用相位點高頻分量的最大值確定主矢量的方向，也即是下一 步捜索窗口的方向； 步驟四：在估計窗口內(nèi)，計算相位主矢量，利用高頻部分主矢量平滑窗口內(nèi)的相位頻 率. 步驟五:在估計窗口內(nèi)，根據(jù)各層分解分量的高頻矢量和構(gòu)建顧及相位細節(jié)變化的非 線性相位曲面模型，將相位表示成非線性相位和線性相位和的形式Φ = Φ linear + Φ ηαη-1 inear，對非線性部分利月對j層分解的Ξ 個高頻方向求矢量和，其中j = l，2...J，ft、敵、ft分別表示對j層高頻的水平、垂直和視 線方向分量求矢量和，矢量求和保證了非線性曲面的特性； 步驟六:對小波分解相位進行逆變換得到相位估計模型，在估計窗口 k X 1內(nèi)，對相位點 4k,i(m，n)，利用可得到估計窗口內(nèi)非線 性相位估計模型，其中，wri(.)為小波逆變換，(iHp為小波分解低頻部分，Σ拓,,為小波分解 高頻部分和。
【文檔編號】G06T9/00GK106097404SQ201610363407
【公開日】2016年11月9日
【申請日】2016年5月27日
【發(fā)明人】不公告發(fā)明人
【申請人】山東科技大學