專利名稱：單次全極化合成孔徑雷達圖像反演數(shù)字高程模型的方法
技術(shù)領(lǐng)域：
本發(fā)明屬于空間遙感與圖像處理技術(shù)領(lǐng)域，具體涉及一種通過合成孔徑雷達圖像獲取數(shù)字高程模型(DEM)的方法。
背景技術(shù)：
合成孔徑雷達(Synthetic Aperture Radar，SAR)是一種主動的機載或星載傳感器。它的基本原理是衛(wèi)星在軌道飛行中于不同位置定時地對同一地物發(fā)射電磁波脈沖信號，同時接收回波信號。在某種意義上，可以認為是延伸了雷達天線的長度，從而大大提高了分辨率。另外，由于合成孔徑雷達還具有全天時、全天候、不受大氣傳播和氣候影響、穿透力強等優(yōu)點，因此在民事和軍事方面的應(yīng)用非常廣泛。
數(shù)字高程模型(Digital Elevation Model)，簡稱DEM。它是用一組有序數(shù)值陣列形式表示地面高程的一種實體地面模型，是數(shù)字地形模型(Digital Terrain Model，簡稱DTM)的一個分支。一般認為，DTM是描述包括高程在內(nèi)的各種地貌因子，如坡度、坡向、坡度變化率等因子在內(nèi)的線性和非線性組合的空間分布，其中DEM是零階單純的單項數(shù)字地貌模型，其他如坡度、坡向及坡度變化率等地貌特性可在DEM的基礎(chǔ)上派生。
目前，通過SAR圖像獲取數(shù)字高程模型(DEM)的方法主要有以下四種雷達測角(Radarclinometry)、雷達立體影像測量(Radargrammetry)、雷達干涉測量(SAP Interferometry)和雷達極化測量(SAR Polarimetry)。其中，雷達立體影像測量和雷達干涉測量必須使用兩幅SAR圖像，前者存在由不同成像視角引起配準困難的問題，而后者的精度雖高，但相關(guān)技術(shù)的實現(xiàn)難度大；雷達測角只需要單幅SAR圖像，直接從地表的后向散射強度反演DEM，但在反演過程中約束條件過多，且忽略了方位向坡度，因此精度不高；而雷達極化測量是三維SAR成像技術(shù)中一種最新的發(fā)展方向，這種方法不同于雷達干涉測量的三維成像技術(shù)，不是利用兩幅圖像的相位差提取高度信息，而是利用雷達回波的極化綜合信息檢測地形高度。
Schuler等人研究表明，當(dāng)原本平坦的地表在方位向(即飛行方向)上發(fā)生傾斜時，其同極化響應(yīng)的極大值位置會發(fā)生改變，即極化方位角ψ將有一個偏移，這個現(xiàn)象可以用于估計地形起伏。由于單次全極化飛行通過極化方位角估計技術(shù)只能獲得方位向坡度，然而還需要已知距離向坡度，才能利用類似干涉SAR相位解纏的最小二乘技術(shù)構(gòu)建整景圖像的DEM。通常情況下，需要研究區(qū)域的參考DEM或者另一次全極化的正交飛行才能獲得與方位向正交的距離向坡度。由于目前極化SAR數(shù)據(jù)以及參考DEM的缺乏，上面提到的方法往往不容易實現(xiàn)。
現(xiàn)有技術(shù)中只通過單次全極化飛行獲取DEM的方法，是通過歐拉角變換得到一個有關(guān)射程角β和水平方位角γ的表達式，再結(jié)合從極化數(shù)據(jù)中估計的極化方位角ψ得到射程角β。又由于方位向坡度SA和距離向坡度SR可用射程角β和水平方位角γ表示，以此得到DEM反演。但此方法在提取水平方位角γ時需要使用復(fù)雜的數(shù)字圖像處理技術(shù)。

發(fā)明內(nèi)容
本發(fā)明的目的在于克服了上述現(xiàn)有技術(shù)的缺點，提出了一種結(jié)合雷達極化和雷達測角技術(shù)使用單次全極化的數(shù)字高程模型(DEM)的反演方法。
按照本發(fā)明的單次全極化合成孔徑雷達圖像反演數(shù)字高程模型的方法，包括如下步驟(1)使用圓極化方法，從極化合成孔徑雷達數(shù)據(jù)估計極化方位角ψ偏移，并對ψ偏移進行去奇異點處理，建立極化方位角ψ偏移與本領(lǐng)域技術(shù)人員公知的距離向坡度角α和方位向坡度角β之間的數(shù)據(jù)關(guān)系；(2)使用雷達測角中的從陰影到形狀技術(shù)，建立地表后向散射強度與距離向坡度角α和方位向坡度角β之間的數(shù)據(jù)關(guān)系；(3)從上述步驟(1)和(2)兩個數(shù)據(jù)關(guān)系，計算距離向坡度α和方位向坡度β；(4)從極化合成孔徑雷達數(shù)據(jù)計算共極化響應(yīng)信號圖，從所述共極化響應(yīng)信號圖中提取權(quán)重圖；(5)從上述步驟(3)和(4)的結(jié)果，用加權(quán)的最小二乘算法求解有關(guān)地面高程的離散泊松方程，從而反演得到合成孔徑雷達觀測區(qū)域的數(shù)字高程模型。
本發(fā)明是基于雷達極化測量方法對方位向坡度敏感，而雷達測角技術(shù)對距離向坡度敏感的基礎(chǔ)，結(jié)合這兩種三維SAR成像技術(shù)的各自特點，從單次全極化SAR圖像反演數(shù)字高程模型(DEM)的方法，從而具有容易實現(xiàn)、可滿足大范圍中等精度地形測繪需求的優(yōu)點。


圖1為奇異點檢測示意圖；圖2為面元的幾何參數(shù)模型；圖3為本發(fā)明的計算程序流程框圖；圖4為中國新疆省烏魯木齊地區(qū)的SIR-C數(shù)據(jù)(C波段Pauli分解RGB圖像)；
圖5為烏魯木齊地區(qū)方位向坡度圖；圖6為烏魯木齊地區(qū)距離向坡度圖；圖7為烏魯木齊地區(qū)的權(quán)重質(zhì)量圖；圖8為烏魯木齊地區(qū)由本發(fā)明反演出的DEM圖；圖9為烏魯木齊地區(qū)的參考DEM圖；圖10為俄羅斯色棱格河地區(qū)的SIR-C數(shù)據(jù)(L波段Pauli分解RGB圖像)；圖11為俄羅斯色棱格河地區(qū)的方位向坡度圖像；圖12為俄羅斯色棱格河地區(qū)的距離向坡度圖像；圖13為俄羅斯色棱格河地區(qū)的權(quán)重質(zhì)量圖；圖14為俄羅斯色棱格河地區(qū)由本發(fā)明反演的DEM圖；圖15為俄羅斯色棱格河地區(qū)的反演DEM的等高線圖。
具體實施例方式
下面參照附圖對本發(fā)明進行具體描述。根據(jù)本發(fā)明的方法，首先，利用圓極化方法從極化SAR圖像中估計出極化方位角ψ偏移，建立極化方位角ψ偏移與地形坡度，即距離向坡度角α和方位向坡度角β之間的數(shù)據(jù)關(guān)系。在只有單次全極化SAR圖像的情況下，假設(shè)成像區(qū)域滿足地表覆蓋的同質(zhì)性以及坡度起伏不大的條件，引入雷達測角中的從陰影到形狀(Shape from Shading)技術(shù)，建立地表的后向散射強度與地表幾何參數(shù)(包括距離向坡度角α和方位向坡度角β)之間的數(shù)據(jù)關(guān)系。這樣，由上述兩個數(shù)據(jù)關(guān)系，就可以直接計算出距離向坡度α和方位向坡度β。最后，從共極化響應(yīng)信號圖中提取權(quán)重圖，用加權(quán)的最小二乘算法求解有關(guān)地面高程的離散泊松方程，從而反演得到DEM。
圖1為奇異點檢測示意圖；圖2為面元的幾何參數(shù)模型；圖3為本發(fā)明的計算程序流程框圖。下面參照附圖1、2和3對本發(fā)明的內(nèi)容具體描述。
I.極化方位角ψ偏移的估計和優(yōu)化(1)輸入用于生成DEM的極化矩陣[C3]或[T3]數(shù)據(jù)，雷達視角φ以及圖像分辨率(包括距離向分辨率Rd和方位向分辨率Ra)。
(2)對輸入的極化矩陣[C3]或[T3]進行極化SAR濾波，可選擇Boxcar和Refined Lee極化濾波；(3)由以上步驟得到濾波后的極化矩陣[C3]或[T3]，利用極化方位角估計方法進行初步的極化方位角估計。Krogager和Czyz提出用右旋和左旋圓極化估計方位角。該方法被Lee等人進一步修改和提煉而得到圓極化方法，具體推導(dǎo)如下極化數(shù)據(jù)一般提供兩種極化類型單視復(fù)(SLC)數(shù)據(jù)和多視復(fù)(MLC)數(shù)據(jù)。由于使用單視復(fù)數(shù)據(jù)可以推導(dǎo)出多視復(fù)數(shù)據(jù)，所以在此我們以單視復(fù)數(shù)據(jù)為例進行初步的極化方位角估計。
極化SAR數(shù)據(jù)可以用散射矩陣(1)表示S=ShhShvSvhSvv---(1)]]>
化發(fā)射或接收。
由于目前SAR衛(wèi)星基本上都是單基系統(tǒng)，即發(fā)射和接收天線為同一天線，滿足天線互易條件，即Shv＝Svh。方位角的旋轉(zhuǎn)通過下式得到S~=cos(ψ)sin(ψ)-sin(ψ)cos(ψ)ShhShvSvhSvvcos(ψ)-sin(ψ)sin(ψ)cos(ψ)---(2)]]>表達式(2)式中矩陣S上面的～表示旋轉(zhuǎn)后的矩陣。
各圓極化分量可以從散射矩陣中得到Srr＝(Shh-Svv+i2Shv)/2Sll＝(Svv-Shh+i2Shv)/2 (3)Srl＝i(Shh+Svv)/2方位角的旋轉(zhuǎn)可以應(yīng)用表達式(2)式得到，整理后為S~rr=Srre-i2ψ]]>S~ll=Sllei2ψ---(4)]]>S~rl=Srl]]>定義一個圓極化向量c=Srr2SrlSll---(5)]]>圓極化基的協(xié)方差矩陣G可以從向量c中得到G=⟨cc*T⟩=⟨|Srr|2⟩2⟨(SrrSrl*)⟩⟨(SrrSll*)⟩2⟨(SrlSrr*)⟩2⟨|Srl|2⟩2⟨(SrlSll*)⟩⟨(SllSrr*)⟩2⟨(SllSrl*)⟩⟨|Sll|2⟩---(6)]]>
應(yīng)用表達式(4)，旋轉(zhuǎn)過的圓極化協(xié)方差矩陣 為G~=⟨|Srr|2⟩2⟨(SrrSrl*)e-i2ψ⟩⟨(SrrSll*)e-i4ψ⟩2⟨(SrlSrr*)ei2ψ⟩2⟨|Srl|2⟩2⟨(SrlSll*)e-i2ψ⟩⟨(SllSrr*)ei4ψ⟩2⟨(SllSrl*)ei2ψ⟩⟨|Sll|2⟩---(7)]]>顯然，對于圓極化協(xié)方差矩陣 ，方位角的變化只影響非對角線元素項的相位。右-右和左-左圓極化項，即(1，3)項，可以用于估計方位角。如果方位角的偏移主要是由方位向坡度引起的，那么⟨S~rrS~ll*⟩=⟨SrrSll*⟩e-i4ψ---(8)]]>由上式可知，只有當(dāng)<SrrSll*>為實數(shù)時，才不會破壞與相位項e-i4ψ有關(guān)的方位角。由于自然界大多數(shù)地物具有反射對稱性，而且可以證明，對于反射對稱介質(zhì)，<SrrSll*>為實數(shù)，交叉和共極化相關(guān)項為0。將式(3)代入<SrrSll*>，并設(shè)包含Shv的項(除了<Shv|2>)為0，則有⟨SrrSll*⟩=(-⟨|Shh-Svv|2⟩+4⟨|Shv|2⟩)/4---(9)]]>該項為實數(shù)，因此<SrrSll*>的幅角為0或π。所以，Shh和Svv之間的相位差不會對方位角估計產(chǎn)生偏差。
表達式(8)的因子4ψ限制了ψ角的范圍為[-π/4，π/4]。為了推導(dǎo)出一般表達式，我們用線性基表示 ⟨S~rrS~ll*⟩=14{⟨-|S~hh-S~vv|2+4|S~hv|2⟩-i4Re(⟨(S~hh-S~vv)S~hv*⟩)}---(10)]]>如果直接應(yīng)用上式將引起偏差，因為對于反射對稱介質(zhì)， 常常大于 分母為負數(shù)。因此，當(dāng)分子接近0時，反正切值接近±π，方位角將為±π/4，而不是像它應(yīng)該的那樣接近于0。為了使方位角與相應(yīng)的方位向坡度相匹配，偏差必須通過加上π來去除。最終，可以得到圓極化估計為ψ=ηη≤π/4η-π/2η>π/4---(11)]]>其中η=14[tan-1(-4Re(⟨(S~hh-S~vv)S~hv⟩)-⟨|S~hh-S~vv|2⟩+4⟨|S~hv|2⟩)+π]---(12)]]>表達式(12)中反正切的計算范圍在[-π，π]。
(4)極化方位角的估計精度取決于主要散射體的方位角偏移是否只由方位向坡度引起。不同散射機制的散射體的混合等效于散射(或相干)矩陣的疊加，當(dāng)旋轉(zhuǎn)相同的方位角，不會帶來誤差，然而，如果有一定數(shù)量的散射具有由于樹枝或其它因素引起的方位角平移，那么就可能產(chǎn)生問題。這時往往由于去極化效果明顯或峰值點不唯一，因此估計出來的極化方位角ψ偏移可能出現(xiàn)角度奇異點。特別地，本發(fā)明采用了去奇異點處理技術(shù)。其具體步驟如下a)設(shè)置一個角度變化閾值T；b)檢測像素點的奇異性。按先行后列的次序依次在整個圖像中掃描每一個像素。如圖1所示，0表示當(dāng)前像素，1，2，3，4表示與其相鄰四個像素點，如果當(dāng)前像素0與1、4中任一像素的角度差值大于給定閾值T時，則判定當(dāng)前像素點為奇異點。這樣做是為了確保該奇異點在奇異區(qū)域內(nèi)部，而不是奇異區(qū)域邊緣的非奇異點。
c)如果當(dāng)前像素點為非奇異點，則返回步驟(b)；否則，進入步驟(d)。
d)用給定窗口大小(如3*3)內(nèi)所有像素值的中值來代替當(dāng)前奇異像素點。
e)重復(fù)(b)-(d)步驟3次即可，最后得到去奇異點后的極化方位角。
II.圖像/地表關(guān)系的建立(1)從步驟I(1)輸入的極化矩陣[C3]或[T3]數(shù)據(jù)中提取HH極化強度圖，并進行SAR圖像濾波(如Lee，F(xiàn)rost，Gamma等)。但注意不要使用增強型濾波器，因為它會保留點目標，而點目標的存在(如孤立的樹木等)無法真實反映地表的實際高度，對于三維高程反演來說是不利的。
(2)計算濾波后整幅HH極化強度圖的平均值，作為平地后向散射強度I0的估計值；(3)Lee和Pottier分別獨立推導(dǎo)出極化方位角ψ與地形坡度和雷達視角具有如下的關(guān)系tanψ=tanβ-tanαcosφ+sinφ---(13)]]>其中，φ為雷達視角(已經(jīng)在步驟I(1)中輸入)，tanβ為方位向坡度SA(下標A表示方位向)，tanα為距離向坡度SR(下標R表示距離向)。由于有兩個未知量(α，β)和一個方程式(8)，要得到坡度信息，需要另一個關(guān)系式。特別地，本發(fā)明采用雷達測角方法中的從陰影到形狀(Shapefrom Shading)技術(shù)得到另一個關(guān)于圖像/地表之間的表達式。
從陰影到形狀最初用于光學(xué)圖像，是一種當(dāng)反射函數(shù)已知的情況下從反射強度圖像恢復(fù)相應(yīng)目標形狀的技術(shù)。R.L.Widley首先將從陰影到形狀技術(shù)應(yīng)用到雷達圖像中，并重新命名為雷達測角。它的原理是從SAR圖像的后向散射強度估計相應(yīng)的地表地形參數(shù)，完全使用單幅SAR圖像來揭示地表坡度信息。
從陰影到形狀技術(shù)使我們能夠由地物的后向散射強度估計相應(yīng)地表的幾何參數(shù)，需要假設(shè)成像地表覆蓋區(qū)域同質(zhì)以及坡度起伏平緩。對于地表的幾何參數(shù)，我們引入面元(Parcel)的概念，它被定義為一個雷達分辨單元所截取的地表，并成像在一個像元內(nèi)。
圖2是雷達圖像采樣點地面較為平坦時(即坡度不是很大時)對應(yīng)的地表形狀近似圖。
這時，雖然L(α)和L(β)構(gòu)成的單位面積并不和斜距垂直，但由于坡度不大，可以近似認為它就是面積A。對于每一個成像單元，包括三個幾何成像參數(shù)雷達入射電磁波的俯角θ(與雷達視角φ互余)距離向分辨率Rd，方位向分辨率Ra，因此，L(α)＝Rdsin(θ)/sin|θ-α|(14)L(β)＝Ra/cos(β) (15)A(α，β)＝L(α)L(β)(16)其中φ、Rd、Ra已經(jīng)在上述步驟I(1)中輸入。另外，每個像素強度不但依賴于成像地面的坡度，而且也取決于其地表覆蓋類型。在光學(xué)系統(tǒng)中有反射模型，在雷達系統(tǒng)中則有后向散射模型。為了建立地面幾何參數(shù)和SAR圖像強度之間的關(guān)系，我們采用簡單而有效的經(jīng)驗后向散射模型--朗伯后向散射模型。由朗伯假設(shè)可以建立起圖像/地面之間的關(guān)系Ir(α，β)＝Kσ0RdRacos2(θ-α)cos2βA(α，β)(17)其中，Ir為后向散射強度，α為距離向坡度角，β為方位向坡度角，θ為雷達入射電磁波的俯角，K為雷達定標常量，σ0為地面覆蓋物的后向散射系數(shù)，A為SAR影像上的一個像素對應(yīng)的地面面積。這個關(guān)系式也是雷達測角的基礎(chǔ)。
由于每個像素的未知量(K，σ0，α，β)眾多，為了從表達式(10)解出地形坡度，我們假設(shè)對于一個平坦地面的后向散射強度I0，地面類型相同，因此σ0在相當(dāng)大范圍內(nèi)為常數(shù)，那么計算平地的后向散射強度I0時α和β為0，即I0＝Ir(0，0)＝Kσ0RdRacos2(θ)(18)此時的平地后向散射強度I0可認為是常量。在通常情況下，我們可以把一幅影像的統(tǒng)計平均值當(dāng)作I0，以此為基礎(chǔ)，建立下面的比例關(guān)系，記為R(α，β)R(α,β)=Ir(α,β)I0=sinθcos2(θ-α)cosβcos2(θ)sin(|θ-α|)---(19)]]>通過求解上式，我們獲得了距離向坡度α的表達式α=θ-arcsin((R′2+4-R′)/2)---(20)]]>其中，R′=R(α,β)cos2θsinθcosβ---(21)]]>
因此，結(jié)合雷達極化和雷達測角的各自特點，從上面兩個數(shù)據(jù)關(guān)系(13)和(20)，就可以直接計算出地表坡度角(α，β)，從而得到方位向坡度SA和距離向坡度SR。
III.通過地表坡度S(x，y)反演DEM(1)由地形坡度反演DEM類似于干涉SAR里的相位解纏，是尋求在M×N矩形網(wǎng)格區(qū)域上具有紐曼邊界約束條件的離散泊松(Possion)方程的解，該泊松方程表示為2h(x，y)＝ρ(x，y)(22)其中，h(x，y)為地面高程，2為拉普拉斯算子2/x2+2/y2。用hij表示在(i，j)處的高程值，源函數(shù)ρij則是由輸入坡度數(shù)據(jù)Sij得到的表面曲率值，ρijR=SijR-Si-1,jR,ρijA=SijA-Si,j-1A(0≤i≤M,0≤j≤N)---(23)]]>為了得到地面高程，必須最小化hij偏導(dǎo)之間的差值，用下面求和公式χ2=Σi,j(hij-hi-1,j-SijR)2+Σi,j(hij-hi,j-1-SijA)2---(24)]]>對hij求導(dǎo)，并令其結(jié)果等于0，則得到(hi-1，j-2hij+hi-1，j)+(hi，j+1-2hij+hi，j+1)＝ρij(25)這就是要求解的泊松方程。其中，離散的源函數(shù)ρijρij=(Si+1,jR-SijR)+(Si,j+1A-SijA)---(26)]]>
現(xiàn)有技術(shù)的求解方法有交替方向隱式迭代(ADI)算法和無權(quán)重完全多重網(wǎng)格(FMG)算法。
FMG是一種經(jīng)典的在較大網(wǎng)格基礎(chǔ)上解決偏微分方程(PDE)的快速方法，在更粗、更小的網(wǎng)格上進行Gauss-Seidel松弛方法的計算。一般地，對于解決N×N網(wǎng)格上離散PDE問題，F(xiàn)MG與基于傅立葉變換的直接方法速度一樣快，而且對網(wǎng)格的大小沒有限制。其關(guān)鍵在于能夠把誤差的低頻成分轉(zhuǎn)換為高頻成份，這樣高頻成分可以很快地用Gauss-Seidel松弛法去除，這是通過把這個問題轉(zhuǎn)移到更粗的網(wǎng)格中實現(xiàn)的。
雖然無權(quán)重的最小二乘法運算速度快，但它有一個明顯的缺陷，即它反演時是穿過坡度的不一致區(qū)，而不是繞過，所以常常會造成地形表面的誤差。這個缺陷可以通過引入權(quán)重來加以彌補。特別地，本發(fā)明采用加權(quán)的FMG算法。一般來說，權(quán)重數(shù)據(jù)wi，j的取值范圍是
，由于極化SAR不存在干涉SAR相位解纏中的相干質(zhì)量圖作為權(quán)重圖，因此我們從共極化信號響應(yīng)圖中定義類似的權(quán)重質(zhì)量圖wi,j=σmax(0,ψ′)-σmin(0,ψ′)σmax(0,ψ′)---(27)]]>其中，σmax(0，ψ′)和σmin(0，ψ′)表示只考慮像素(i，j)在線性極化基(χ＝0)情況下的共極化信號的最大/最小響應(yīng)。
目標的后向散射截面(或接收功率)是入射波極化狀態(tài)的函數(shù)，即是入射波的極化橢圓率和極化方位角的函數(shù)。當(dāng)發(fā)射和接收天線的極化狀態(tài)一致時極化響應(yīng)稱為共極化響應(yīng)，當(dāng)它們極化狀態(tài)正交時的極化響應(yīng)稱為交叉極化響應(yīng)。這里，只考慮共極化響應(yīng)，那么目標的后向散射截面可表示為σ(χ,ψ)=1cos2χcos2ψ′cos2χsin2ψ′sin2χTK1cos2χcos2ψ′cos2χsin2ψ′sin2χ---(28)]]>其中，K為Stokes矩陣；χ為橢圓率角，取值范圍在 ψ′為極化方位角，取值范圍在 式(28)遍歷整個(χ，ψ′)平面，即可獲得目標的共極化響應(yīng)圖。當(dāng)目標散射接近單一的表面散射時，wi，j越大；當(dāng)目標散射接近體散射，或者說，去極化效應(yīng)增強時，會引起共極化信號響應(yīng)圖基底的增高，因此wi，j越小。
這樣，無權(quán)最小二乘問題就轉(zhuǎn)換為加權(quán)最小二乘問題?？紤]如下的加權(quán)最小二乘問題min{Σi,jwijR(hi+1,j-hi,j-Si,jR)2+Σi,jwijA(hi,j+1-hi,j-Si,jA)2}---(29)]]>梯度權(quán)重wijR和wijA定義如下wijR=min(wi+1,j2,wi,j2)wijA=min(wi,j+12,wi,j2)---(30)]]>于是，這個問題的最小二乘解hi，j定義為ρ^ij=wi,jR(hi+1,j-hi,j)-wi-1,jR(hi,j-hi-1,j)]]>+wijA(hi,j+1-hi,j)-wi,j-1A(hi,j-hi,j-1)---(31)]]>
其中， 是加權(quán)相位拉普拉斯操作數(shù)，ρ^ij=wijRSijR-wi-1,jRSi-1,jR+wijASi,jA-wi,j-1ASi,j-1A---(32)]]>經(jīng)典的Gauss-Seidel松弛法通過對方程hi,j=wijRhi+1,j+wi-1,jRhi-1,j+wijAhi,j+1+wi,j-1Ahi,j-1-ρ^ijvi,j---(33)]]>多次迭代求解hi，j。
表達式(33)中，vi,j=wijR+wi-1,jR+wijA+wi,j-1A---(34)]]>(2)使用步驟I(1)得到的距離向分辨率Rd和方位向分辨率Ra、步驟II(3)得到的方位向坡度SA和距離向坡度SR和步驟III(1)得到的權(quán)重質(zhì)量圖W，使用加權(quán)FMG算法就可以得到DEM。
加權(quán)FMG算法通過坡度數(shù)據(jù)擬合地形表面的高程，對數(shù)據(jù)“損壞”區(qū)通過設(shè)置0權(quán)重的方式，防止誤差的傳遞，然后利用多分辨率的網(wǎng)格上的Gauss-Seidel松弛來加速收斂。
實驗1新疆烏魯木齊地區(qū)采用美國加州理工學(xué)院噴氣推進實驗室(JPL)提供的1994年10月8日SIR-C對中國天山地區(qū)觀測的C波段數(shù)據(jù)，具體參數(shù)如下所示實驗地點中國天山成像時間1994/10/08 22:06:45.647
圖像中心緯度(度)43.6740608圖像中心經(jīng)度(度)88.1214600圖像中心入射角(度) 24.098雷達頻率(GHz) 5.298產(chǎn)品類型單視復(fù)(SLC)行距(米)5.0520549像素間隔(米)6.6624818每行圖像的像素數(shù)1100圖像行數(shù)9896地圖投影描述子 斜距軌道方向升軌圖4是從原圖中截取的大小為273×367像素的Pauli分解合成圖(即R為HH-VV，G為HV，B為HH+VV)。距離向分辨率為16.3177m，方位向分辨率為15.1562m(方位向做了三視處理)，其中心位置為(88.0322E，43.5610N)，位于中國新疆省烏魯木齊地區(qū)。
從圖4、圖8以及圖9可知，圖8反演的DEM與圖4中的主要地物的走向基本一致，這也在圖8的參考DEM中再次得到證實，圖8中左邊的四條深溝與圖9中的四條深溝的明暗紋理以及走向也是非常符合的。
為了更客觀地評價DEM反演精度，將圖8中反演的DEM配準到圖9中的參考DEM上，在兩者的共同區(qū)域均勻選取20個DEM控制點，得到的詳細精度分析如表1所示表1反演DEM中控制點的精度分析


根據(jù)均勻分布20個的控制點計算得到平均高程偏差為-0.91m，均方根誤差為20.29m。這個精度完全可以滿足大范圍地形測繪的需求。由表1可見，本發(fā)明所測量得到的高程與參考高程相比存在著一定的誤差，其誤差來源一個來自測量DEM與參考DEM之間的配準誤差，二是由于極化數(shù)據(jù)在C波段中采集，其較短的波長對植被覆蓋物穿透深度有限，使得測量的高程與實際地表的高程有所偏差，而使用波長更長的L波段則有望提高DEM精度。
實驗2俄羅斯色棱格河地區(qū)采用美國加州理工學(xué)院噴氣推進實驗室(JPL)提供的1994年10月7日SIR-C對俄羅斯色棱格河地區(qū)觀測的L波段數(shù)據(jù)，具體參數(shù)如下所示實驗地點俄羅斯色棱格河地區(qū)成像時間1994/10/07 22:33:30.642圖像中心緯度(度)52.1651077圖像中心經(jīng)度(度)106.6797943圖像中心入射角(度) 24.569雷達頻率(GHz) 1.249產(chǎn)品類型單視復(fù)(SLC)行距(米)5.0520549像素間隔(米)6.6624818
每行圖像的像素數(shù)1100圖像行數(shù)9030地圖投影描述子 斜距軌道方向升軌圖10是從原圖中截取的大小為562×380像素的Pauli分解合成圖(即R為HH-VV，G為HV，B為HH+VV)。距離向分辨率為16.3177m，方位向分辨率為15.1562m(方位向做了三視處理)，其中心位置為(107.0003E，52.2443N)，位于俄羅斯色棱格河地區(qū)。
由于該地區(qū)在俄羅斯地區(qū)，無法得到精度較高可與反演DEM相比較的參考DEM。但是，從圖10和圖12可以很直觀容易地看出，圖14反演得到的DEM與圖像10中地形起伏吻合得很好，圖10右下部分一些較細微的地形起伏變化，也在反演DEM圖14中有良好的體現(xiàn)。
權(quán)利要求
1.一種單次全極化合成孔徑雷達圖像反演數(shù)字高程模型的方法，該方法利用了極化方位角ψ偏移、距離向坡度角α和方位向坡度角β，其特征在于，包括如下步驟(1)使用圓極化方法，從極化合成孔徑雷達數(shù)據(jù)估計極化方位角ψ偏移，并對ψ偏移進行去奇異點處理，建立極化方位角ψ偏移與距離向坡度角α和方位向坡度角β之間的數(shù)據(jù)關(guān)系；(2)使用雷達測角中的從陰影到形狀技術(shù)，建立地表后向散射強度與距離向坡度角α和方位向坡度角β之間的數(shù)據(jù)關(guān)系；(3)從上述步驟(1)和(2)兩個數(shù)據(jù)關(guān)系，計算距離向坡度α和方位向坡度β；(4)從極化合成孔徑雷達數(shù)據(jù)計算共極化響應(yīng)信號圖，從所述共極化響應(yīng)信號圖中提取權(quán)重圖；(5)從上述步驟(3)和(4)的結(jié)果，用加權(quán)的最小二乘算法求解有關(guān)地面高程的離散泊松方程，從而反演得到合成孔徑雷達觀測區(qū)域的數(shù)字高程模型。
2.根據(jù)權(quán)利要求1所述的方法，其特征在于，對極化方位角ψ偏移的估計進行去奇異點處理，其中，極化方位角ψ偏移的估計為ψ=ηη≤π/4η-π/2η>π/4---(11)]]>其中η=14[tan-1(-4Re(⟨(S~hh-S~vv)S~hv*⟩)-⟨|S~hh-S~vv|2⟩+4⟨|S~hv|2⟩)+π]---(12)]]>表達式(12)中反正切的計算范圍為[-π，π]， 是從極化合成孔徑雷達數(shù)據(jù)得到的；奇異點處理的具體步驟如下a)設(shè)置一個角度變化閾值T；b)檢測像素點的奇異性，按先行后列的次序依次在整個圖像中掃描每一個像素，0表示當(dāng)前像素，1，2，3，4表示與其相鄰四個像素點，如果當(dāng)前像素0與1、4中任一像素的角度絕對差大于給定閾值T時，則判定當(dāng)前像素點為奇異點；c)如果當(dāng)前像素點為非奇異點，則返回步驟(b)；否則，進入步驟(d)；d)用給定窗口大小內(nèi)所有像素值的中值來代替當(dāng)前奇異像素點；e)重復(fù)(b)-(d)步驟3次，最后得到去奇異點后的極化方位角。
3.根據(jù)權(quán)利要求1所述的方法，其特征在于，采用雷達測角方法中的從陰影到形狀技術(shù)得到距離向坡度α的表達式α=θ-arcsin((R′2+4-R′)/2)---(20)]]>其中R′=R(α,β)cos2θsinθcosβ---(21)]]>雷達入射電磁波的俯角θ由輸入的SAR數(shù)據(jù)得到；又因為極化方位角ψ與地形坡度和雷達視角具有如下的關(guān)系tanψ=tanβ-tanαcosφ+sinφ---(13)]]>其中，φ為雷達視角，tanβ為方位向坡度SA，tanα為距離向坡度SR，從表達式(13)和(20)，就可以直接計算出地表坡度角(α，β)，從而得到方位向坡度SA和距離向坡度SR。
4.根據(jù)權(quán)利要求1所述的方法，其特征在于，從共極化信號響應(yīng)圖中定義權(quán)重質(zhì)量圖wi,j=σmax(0,ψ′)-σmin(0,ψ′)σmax(0,ψ′)---(27)]]>其中σmax(0，ψ′)和σmin(0，ψ′)表示像素(i，j)在線性極化基(x＝0)情況下的共極化信號的最大/最小響應(yīng)，共極化信號響應(yīng)圖的計算表達公式為σ(χ,ψ′)=1cos2χcos2ψ′cos2χsin2ψ′sin2χTK1cos2χcos2ψ′cos2χsin2ψ′sin2χ---(28)]]>其中，K為Stokes矩陣，x為橢圓率角，ψ′為極化方位角。
5.根據(jù)權(quán)利要求1所述的方法，其特征在于采用加權(quán)的最小二乘算法求解有關(guān)地面高程的離散泊松方程。
全文摘要
本發(fā)明是結(jié)合雷達極化和雷達測角技術(shù)的各自特點，從單次全極化合成孔徑雷達圖像反演數(shù)字高程模型的方法。使用圓極化方法從極化SAR數(shù)據(jù)估計極化方位角ψ偏移，建立極化方位角ψ偏移與距離向坡度角α和方位向坡度角β之間的數(shù)據(jù)關(guān)系。在只有單次全極化SAR圖像的情況下，假設(shè)成像區(qū)域滿足地表覆蓋的同質(zhì)性以及坡度起伏不大的條件，引入雷達測角中的從陰影到形狀技術(shù)，建立地表后向散射強度與地表幾何參數(shù)(包括距離向坡度角α和方位向坡度角β)之間的數(shù)據(jù)關(guān)系。這樣，由上述數(shù)據(jù)關(guān)系，計算出距離向坡度α和方位向坡度β。從共極化響應(yīng)信號圖中提取權(quán)重圖，用加權(quán)的最小二乘算法求解有關(guān)地面高程的離散泊松方程，從而反演得到DEM。
文檔編號G01S13/00GK101078769SQ20061008127
公開日2007年11月28日 申請日期2006年5月25日 優(yōu)先權(quán)日2006年5月25日
發(fā)明者張紅, 陳曦, 王超 申請人:中國科學(xué)院中國遙感衛(wèi)星地面站