一種地形測量方法

2023-05-06 00:49:01 1

專利名稱：一種地形測量方法
技術領域：
本發明涉及信息獲取與處理技術領域，特別是一種利用單路徑全極化合成孔徑雷達復圖像數據的地形測量方法。

背景技術：
利用合成孔徑雷達的全極化信息進行地形測量是一種新的技術，如圖1所示。全極化雷達通過在正交方向上發射和接收電磁波信號，可以形成四個通道的全極化復圖像數據SHH、SHV、SVH和SVV，如圖2所示。在後向散射坐標系下，當目標滿足互易性假設條件時有SHV＝SVH，故一般僅使用三個通道的數據。SHV代表利用天線的垂直極化模式(V，電場方向垂直於電磁波入射方向，且在入射平面以內)來接收水平極化模式(H，電場方向垂直於電磁波入射方向，且垂直於入射平面)發射信號的回波圖像複數據。全極化雷達不僅可以記錄回波的幅度和相位信息，同時還可以記錄回波的矢量信息。由於合成孔徑雷達圖像強度對地表的地距向坡度變化非常敏感，而回波的極化信息會受到地表的方位向坡度變化影響，因此可以通過單次飛行獲取的單路徑全極化合成孔徑雷達復圖像數據獲取一對正交方向的坡度圖方位向坡度圖和地距向坡度圖。從這兩幅坡度圖可以推導任意方向的地形坡度值以及觀測場景內的數字高程模型。
目前，國內外就利用全極化合成孔徑雷達進行地形測量的理論與方法開展了一些研究。1993年，Dale L.Schuler等人(D.L Schuler，J.S.Lee，Karl W.Hoppel.「Polarimetric SAR Image Signatures of the Ocean andGulf Stream Features」，IEEE Transactions on Geosciences and RemoteSensing，Vol.31，No.6，pp1210-1221，Nov.1993.)發現海洋表面波浪會引起全極化雷達同極化特徵圖偏移，且偏移量的大小與方位向坡度相關。2000年，Dale L.Schuler等人(D.L.Schuler，J.S.Lee.「Thomas.L.Ainsworth and M.R.Grunes.Terrain Topography Measurement UsingMultipass Polarimetric Synthetic Aperture Radar Data」，Radio Science，Vol.35，No.3，pp813-832，2000.)發現了這種關係內在的數學表達式，建立了地形—極化方位角偏移模型。但是模型方程中包含的兩個未知數方位向坡度和地距向坡度無法僅利用一個已知的極化方位角偏移值來求解。為解決這個問題，他們用兩條正交飛行路徑上的全極化合成孔徑雷達復圖像數據計算兩個方向上的極化方位角偏移值，並反演出地形信息。其缺點是需要兩條路徑上的全極化數據，增加圖像配準步驟，且數據獲取難度較大。
為了解決利用單路徑全極化合成孔徑雷達數據完成地形測量的問題，金亞秋等人(JIN Ya-Qiu and LUO Lin.「Terrain Topographic InversionUsing Single-Pass Polarimetric SAR Image Data」，SPIE Vol.4894，pp425-433，2003.)於2003年發表了利用雷達圖像灰度信息的方法。該文建立了雷達圖像中的紋理方向與地表正交坡度之間的數值關係模型，隱含限定了「沿紋理方向的高程變化率為零」的條件，並與地形—極化方位角偏移模型方程聯立來解得方位向坡度和地距向坡度。該方法依賴於精準的圖像形態學處理，對紋理的識別和提取要求較高，特別是找到符合限定條件的紋理的難度較大。
2006年，CHEN Xi等人(CHEN Xi，ZHANG Hong，WANG Chao.「ANew Method for DEM Generation Using A Single POLSAR Flight Pass」，IGRASS』2006.)引入了雷達測角技術及陰影成型技術中使用的朗伯反射模型，建立了雷達圖像亮度和正交坡度之間的數值關係。將該模型與地形—極化方位角偏移模型方程聯立來解得方位向坡度和地距向坡度。由於朗伯反射模型中的未知數較多，使用前還需要估計無地形變化時的區域後向散射係數。但後向散射係數估計問題本身就需要已知地表坡度，因此整個過程還需要迭代進行，直到終止條件得到滿足。該算法需要引入較強的假設條件，且計算效率較低。
2000年，Jong-Sen Lee等人(J.S.Lee，D.L.Schuler，T.L.Ainsworth.「Polarimetric SAR Data Compensation for Terrain Azimuth Slope Variation」，IEEE Transactions on Geoscience and Remote Sensing，Vol.38，No.5，pp2153-2163，2000.)首次較系統地提出極化補償的概念，揭示了利用極化方位角偏移值可以補償方位向坡度變化對極化合成功率的影響，從而提高地物參數反演的精度。該技術雖然沒有被用於解決地形測量問題，但卻是本發明需要使用的技術之一。
1996年，Mark D.Pritt(Mark D.Pritt.「Phase Unwrapping by Means ofMultigrid Techniques for Interferometric SAR」.IEEE Transactions onGeosciences and Remote Sensing，Vol.34，No.3，pp728-738，1996.)將完全多重網格算法用於解決相位展開中的離散泊松方程問題。本發明也將由坡度估計高度的問題抽象為解離散泊松方程問題，因此也使用了完全多重網格算法。


發明內容
為了解決現有技術的什麼問題，本發明的目的在於利用可以由單路徑的全極化合成孔徑雷達斜距復圖像數據快速地獲得中等精度方位向坡度、地距向坡度和數字高程模型，無需另外附加斜地轉換處理步驟，為此，本發明提供一種利用單路徑全極化合成孔徑雷達復圖像數據的地形測量方法。
為達成所述目的，本發明提供一種地形測量方法步驟如下 步驟A利用單路徑全極化合成孔徑雷達復圖像數據並採用圓極化算法計算，得到照射區域內每一個地表單元的極化方位角偏移中間值為 並進行值域映射為得到照射區域內地表單元的極化方位角偏移值θ(x，y)，式中SHH、SHV、SVH和SVV為四幅全極化復圖像數據，每一幅圖像為一個二維矩陣，圖像中每一個像素點在矩陣中的位置為(x，y)，表示圖像中每一個地表單元的行、列坐標，從0開始計數； 步驟B利用雷達載體平臺高度HA、近地端入射角度η0及斜距解析度RSL來計算每一個地表單元的雷達入射角度η(x，y)、方位向解析度Ra(x，y)和地距解析度Rg(x，y)； 步驟C對地表單元的極化方位角偏移值和雷達入射角使用補償朗伯算法，計算每一個地表單元的方位向坡度猜測值ωg(x，y)和地距向坡度猜測值βg(x，y)； 步驟D對地表單元的方位向解析度及地距向解析度、方位向坡度猜測值和地距向坡度猜測值使用完全多重網格算法，計算每一個地表單元的相對高程猜測值Hg(x，y)； 步驟E對地表單元的方位向解析度、地距向解析度以及相對高程猜測值做二維差分運算，計算每一個地表單元的方位向坡度修正值ω(x，y)和地距向坡度修正值β(x，y)； 步驟F對地表單元的方位向解析度及地距向解析度、方位向坡度修正值和地距向坡度修正值使用完全多重網格算法，結合一個外部地表單元參考點HT(r，c)計算每一個地表單元的絕對高程估計值H(x，y)。
本發明的有益效果本發明提供一種利用單路徑全極化合成孔徑雷達復圖像數據的地形測量方法，由單路徑的全極化合成孔徑雷達斜距復圖像數據快速地獲得中等精度方位向坡度、地距向坡度和數字高程模型，無需另外附加斜地轉換處理步驟，這種方法的特點在於對雷達數據獲取的限制較少、對外部數據的依賴較少、自動化程度較高以及整體流程計算效率較高，從而使其在地形可視化應用領域具有較大的應用價值。本發明避免了Dale L.Schuler方法中「雙路徑數據獲取難度大，後期配準處理複雜」的問題，突破了金亞秋方法中「沿紋理方向的高程變化率為零」這一條件的限制，解決了CHENXi方法中由朗伯反射模型未知變量較多帶來的「假設條件過強，計算效率較低」的問題。本發明能夠簡單、快速地獲取觀測區域的中等精度方位向坡度圖、地距向坡度圖和數字高程圖。



圖1是全極化雷達信號收發模式示意圖。
圖2是四通道全極化復圖像數據結構示意圖。
圖3a是本發明的主流程圖。
圖3b是圖3a中步驟C的子流程圖。
圖3c是圖3a中步驟D的子流程圖。
圖3d是圖3a中步驟F的子流程圖。
圖3e是圖3c中步驟DP的子流程圖。
圖3f是圖3c中步驟DS的子流程圖。
圖3g是圖3f中步驟DS8的子流程圖。
圖4是地形—極化方位角偏移模型幾何示意圖。
圖5是朗伯反射模型幾何示意圖。
圖6a是機載全極化合成孔徑雷達復圖像數據SHH的dB功率圖(全景)。
圖6b是機載全極化合成孔徑雷達復圖像數據SHV的dB功率圖。
圖6c是機載全極化合成孔徑雷達復圖像數據SVV的dB功率圖。
圖7a是本發明輸出的極化方位角偏移值圖(部分區域)。
圖7b是極化方位角偏移地表真實值圖。
圖7c是圖7a和圖7b中的剖面曲線1對照圖。
圖7d是圖7a和圖7b中的剖面曲線2對照圖。
圖8a是本發明輸出的修正後的方位向坡度圖(部分區域)。
圖8b是方位向坡度地表真實值圖。
圖8c是本發明輸出的修正後的地距向坡度圖。
圖8d是地距向坡度地表真實值圖。
圖9a是圖8a和圖8b的剖面曲線對比圖。
圖9b是圖8c和圖8d的剖面曲線對比圖。
圖10a是絕對高程地表真實值圖(全景)。
圖10b是本發明輸出的絕對高程值圖。
圖11a是絕對高程地表真實值圖(部分區域)。
圖11b是本發明輸出的絕對高程圖。
圖12a是圖11a和圖11b的剖面曲線1對比圖。
圖12b是圖11a和圖11b的剖面曲線2對比圖。
圖12c是圖11a和圖11b的剖面曲線3對比圖。
圖12d是圖11a和圖11b的剖面曲線4對比圖。

具體實施例方式 下面結合附圖詳細說明本發明技術方案中所涉及的各個細節問題。應指出的是，所描述的實施例僅旨在便於對本發明的理解，而對其不起任何限定作用。
本發明提出的利用單路徑全極化合成孔徑雷達復圖像數據的地形測量方法的流程圖如圖3所示。在地形測量進行步驟的具體闡述之前，先對發明中所用到的重要符號進行統一說明 全極化復圖像數據包含4個複數矩陣，每個矩陣的大小相同M1行，M2列。這四個矩陣分別命名為SHH、SHV、SVH和SVV。SHH(x，y)代表SHH矩陣中第x行、第y列元素(從0開始計數)。在本發明中，全極化數據處理默認在後向散射坐標系中且目標滿足互易性假設即SHV＝SVH，因此輸入複數矩陣實際為3個。θ是產生的極化方位角偏移矩陣；HA為雷達載體平臺高度，η0為近地端入射角度，RS為斜距解析度；η為雷達入射角度矩陣；Ra和Rg分別是方位向解析度和地距向解析度矩陣；ωg為方位向坡度猜測值矩陣，βg為地距向坡度猜測值矩陣；Hg為相對高程猜測值矩陣；ω為方位向坡度修正值矩陣，β為地距向坡度修正值矩陣，NC為最粗糙尺寸；每一個地表單元的絕對高程猜測值為H(x，y)＝HT(r，c)-Hg(r，c)+Hg(x，y)；(r，c)為外部參考點在高程猜測值矩陣中的行列坐標，Hg(r，c)為外部參考點處的相對高程猜測值，HT(r，c)為外部參考點處的絕對高程值。H為絕對高程估計值矩陣，HT(r，c)為外部參考點絕對高程值，其中(r，c)代表H矩陣中的行列坐標(從0開始計數)。以上所述的變量名如果代表矩陣，則它們的尺寸都與SHH矩陣相同。
圖4給出了地形—極化方位角偏移模型幾何示意圖，圖中VH代表全極化合成孔徑雷達，X軸代表雷達平臺運動方向即方位向，Y軸代表地距方向，Z軸代表天底方向，XY平面代表地表平面，YZ平面代表雷達入射平面，坐標原點處的灰色方形區域代表具有一定坡度值的散射單元，N是它的表面法向量，η代表雷達入射角。圖5給出了朗伯反射模型幾何示意圖，平行四邊形是由雷達解析度決定的一個單元截取的一小塊地表，坐標系定義與圖4相同。地表單元尺寸L(ω)和L(β)表示為地距向解析度Rg和方位向解析度Ra的函數。圖4和圖5中所使用的坡度、入射角等變量對應於前述同變量名矩陣中的一個元素，如同圖中的地表面元對應於復圖像數據矩陣中的一個元素一樣。
地形—極化方位角偏移模型為 朗伯反射模型為 其中I(x，y)代表極化合稱圖像中每一個地表單元的亮度值，ω(x，y)代表每一個地表單元的方位向坡度，β(x，y)代表每一個地表單元的地距向坡度，K為雷達系統的定標常數，σ0(x，y)為每一個地表單元的歸一化後向散射係數，Rg(x，y)為每一個地表單元的地距向解析度，Ra(x，y)為每一個地表單元的方位向解析度，η(x，y)為每一個地表單元的入射角度，θ(x，y)為每一個地表單元的極化方位角偏移值。
在所述所的完全多重網格算法中需要不斷對高程值矩陣HN執行抽取和插值運算，因此在進行多重網格算法之前需要對輸入矩陣補零，為了使矩陣尺寸(行、列數)都拓展為2的整數次冪。拓展後的行列數中較小的值定義為矩陣尺寸N。運算完畢後再將矩陣恢復為補零前的尺寸。本發明中其它帶有下標N的矩陣變量名也具有相同的約定。方位向和地距向相鄰像素點間的高程差分別為Δa和Δg；ρ為源函數項矩陣；A·為(17)式定義的矩陣元素運算符；R(·)為全加權限制運算符，P(·)為雙線性延拓運算符；←為賦值運算符；v1，v2為鬆弛運算次數。
本發明提出的利用單路徑全極化合成孔徑雷達復圖像數據的地形測量方法的流程圖如圖3a所示，下面將針對具體步驟A-F進行詳細說明 步驟A將三幅全極化復圖像數據SHH、SHV和SVV代入圓極化算法計算公式(1)並進行值域映射(2)，得到每一個像素點對應的每一個地表單元的極化方位角偏移值θ(x，y)，Re(·)為取複數實部運算，*代表復共軛，arcctan代表反正切運算。每一幅圖像為一個二維矩陣，圖像中每一個像素點在矩陣中的位置為(x，y)，表示圖像中每一個地表單元的行、列坐標，從0開始計數。
步驟B利用雷達載體平臺高度HA、近地端入射角度η0及斜距解析度RSL來計算每一個地表單元的雷達入射角度η(x，y)、方位向解析度Ra(x，y)和地距解析度Rg(x，y)； 利用已知的雷達載體平臺高度HA、近地端入射角度η0計算近地端地距 d(0，0)＝HA·tan(η0) (3) 令入射角矩陣第一行第一列元素η(0，0)＝η0，開始計算入射角矩陣第一行其餘元素，i＝1，2，…，M2-1 入射角矩陣其餘各行的元素與第一行元素相同 η(j，i)＝η(0，i)，j＝1，2，…，M1-1 (5) 方位向解析度矩陣中的每一個元素Ra(x，y)都賦值為已知的系統方位向解析度參數RA；地距向解析度矩陣中的每一個元素用(6)式計算 步驟C對地表單元的極化方位角偏移值和雷達入射角使用補償朗伯算法，計算每一個地表單元的方位向坡度猜測值ωg(x，y)和地距向坡度猜測值βg(x，y)，請參見圖3b步驟C的子流程，主要包含以下步驟 步驟CS1對四幅全極化復圖像數據每一個地表單元的SHH(x，y)、SHV(x，y)、SVH(x，y)和SVV(x，y)進行組合，得到每一個地表單元的斯託克斯矩陣元素m22(x，y)，m23(x，y)和m33(x，y) m22(x，y)＝0.25(|SHH(x，y)|2+|SVV(x，y)|2-2|SHV(x，y)|2) (7) m23(x，y)＝0.5Re(SHH(x，y)S*HV(x，y)-SVV(x，y)S*HV(x，y)) (8) m33(x，y)＝0.5(|SHV(x，y)|2+Re(SHH(x，y)S*VV(x，y))) (9) 步驟CS2對每一個地表單元的斯託克斯矩陣元素m22(x，y)，m23(x，y)和m33(x，y)進行極化補償後，獲得補償後的極化合成圖像每一個地表單元的亮度為I′(x，y) 式中補償後的每一個地表單元的極化方位角偏移值θ(x，y)為0；在地形—極化方位角偏移模型中等價於每一個地表單元的方位向坡度值ω(x，y)為0；極化補償等價於令朗伯反射模型變為 步驟CS3保持原有極化方位角偏移值不變的極化合成圖像每一個地表單元的亮度為I(x，y) I(x，y)＝m22(x，y) (11) 步驟CS4利用補償及未補償的極化合成圖像每一個地表單元的亮度進行每一個地表單元的方位向坡度猜測值ωg(x，y)的計算為 (12)； 令每一個地表單元的方位向坡度猜測值的符號(正號或負號)等於相同地表單元對應像素點處的極化方位角偏移值θ(x，y)的符號，accos為反餘弦運算子。
步驟CS5利用每一個地表單元的極化方位角偏移值及方位向坡度猜測值進行每一個地表單元的地距向坡度猜測值βg(x，y)計算 步驟D請參見圖3c是步驟D的子流程圖，對地表單元的方位向解析度及地距向解析度、方位向坡度猜測值和地距向坡度猜測值使用完全多重網格算法，計算每一個地表單元的相對高程猜測值Hg(x，y)，主要包含以下步驟 步驟DP請參見圖3e是步驟DP的子流程圖，利用每一個地表單元的方位向解析度、地距向解析度、方位向坡度和地距向坡度來計算每一個地表單元的源函數項，主要包含以下步驟 步驟DP1利用地表單元的方位向解析度及地距向解析度、方位向坡度修正值和 地距向坡度修正值計算每一個地表單元的方位向高程差 Δa(x，y)＝Ra(x，y)·tan(ωg(x，y)) (14) 和地距向高程差 Δg(x，y)＝Rg(x，y)·tan(βg(x，y))(15) 步驟DP2利用每一個地表單元的方位向高程差和地距向高程差計算每一個地表單元的源函數項 ρ(x，y)＝Δa(x+1，y)-Δa(x，y)+Δg(x，y+1)-Δg(x，y)(16) 執行邊界校正令源函數項矩陣的第一行元素為對應同一列上第二行元素的相反數；令源函數項矩陣的第一列元素為對應同一行上第二列元素的相反數。
步驟DS將相對高程猜測矩陣中每一個地表單元的相對高程猜測值Hg(x，y)初始化為0；以源函數項矩陣和相對高程猜測矩陣為輸入進行「完全多重網格計算」，計算結果將更新相對高程猜測矩陣中每一個地表單元的相對高程猜測值Hg(x，y)。若循環次數NF的典型值設定為3，則再重複「完全多重網格計算」2次。從第二次開始，輸入的相對高程猜測矩陣中每一個元素的值為前一次「完全多重網格計算」更新的每一個地表單元的相對高程猜測值。由於一次「完全多重網格計算」過程中需要改變相對高程猜測值矩陣Hg和源函數項矩陣ρ的大小，因此令帶有下標N的矩陣HN和ρN分別表示不同分辨單元大小的高程值矩陣和源函數項矩陣。例如，若Hg和ρ是大小為M1行M2列(設M1小於M2)的矩陣(這兩個矩陣的大小總是相同的)，在進行一次「完全多重網格計算」之前，首先令N＝M1，開始進行第一次「完全多重網格計算」之後，下標N將陸續變為N＝M1/2，N＝M1/4，N＝M1/8，…，N＝Nc。其中Nc為預先設定的最粗糙矩陣尺寸，在這裡令Nc＝4。

表示一個行數和列數都減小為原始矩陣

的矩陣，而矩陣中每一個元素所代表的分辨單元大小也相應變為原來的8×8(行和列都乘以8)倍。
一次「完全多重網格計算」，請參見圖3f是步驟DS的子流程圖，主要包含以下步驟 步驟DS1對高程值矩陣HN的行列數N進行判斷如果符合最粗糙尺寸定義即N≤NC，則執行步驟DS8；如果不符合最粗糙尺寸定義，則執行步驟DS2； 步驟DS2對行列數為N的源函數項矩陣中每一個地表單元的源函數項和行列數為N的高程猜測值矩陣中每一個地表單元的初始高程猜測值進行殘留誤差限制計算，得到行列數為N/2的源函數項矩陣中每一個地表單元的源函數項ρN/2←R(ρN-A·HN)，其中 A·H(x，y)＝H(x+1，y)+H(x-1，y)+H(x，y+1)+H(x，y-1)-4H(x，y) (17) 全加權限制算子R(·)為 其中f(x，y)輸入全加權限制算子的矩陣元素，c(x，y)為輸出矩陣元素。
步驟DS3對行列數為N/2的源函數項矩陣中第一行和第一列中的每一個地表單元的源函數項進行邊界校正ρN/2(0，y)＝-ρN/2(1，y)，ρN/2(x，0)＝-ρN/2(x，1)。
步驟DS4對行列數為N的高程猜測值矩陣中每一個地表單元的初始高程猜測值進行限制操作，得到行列數為N/2的高程猜測值矩陣中每一個地表單元的初始高程猜測值HN/2←R(HN)。
步驟DS5對行列數為N/2的高程猜測值矩陣中第一行和第一列中每一個地表單元的初始高程猜測值進行邊界校正HN/2(0，y)＝-HN/2(1，y)，HN/2(x，0)＝-HN/2(x，1)。
步驟DS6將行列數為N/2的高程猜測值矩陣中的每一個地表單元的初始高程猜測值以及源函數項代入步驟DS1進行遞歸計算，更新行列數為N/2的高程猜測值矩陣中每一個地表單元的高程猜測值N←N/2，HN/2←FMG(HN/2，ρN/2)。
步驟DS7對行列數為N/2的高程猜測值矩陣中的每一個地表單元的高程猜測值使用延拓算子得到行列數為N的每一個地表單元的高程延拓值，與行列數為N的每一個地表單元的初始高程猜測值相加，更新行列數為N的高程猜測值矩陣中每一個地表單元的高程猜測值HN←HN+P(HN/2)，雙線性延拓算子P(·)為 f(2x，2y)＝c(x，y) 其中c(x，y)為輸入雙線性延拓算子的矩陣元素，f(x，y)為輸出矩陣元素。
步驟DS8對行列數為N的高程猜測值矩陣中的每一個地表單元的高程猜測值以及行列數為N的源函數項矩陣中的每一個地表單元的源函數項使用V形多重網格算法，更新行列數為N的高程猜測值矩陣中每一個地表單元的高程猜測值HN←MV(HN，ρN)，請參見圖3g是步驟DS8的子流程圖，主要包含以下步驟 步驟DS81對行列數為N的高程猜測值矩陣HN進行n次鬆弛操作，n為任意自然數，在本實例中取值為9，一次鬆弛運算表示為 H(x，y)＝[(H(x+1，y)+H(x-1，y)+H(x，y+1)+H(x，y-1))-ρ(x，y)]/4 (20) 步驟DS82對矩陣尺寸N進行判斷如果符合最粗糙尺寸定義，則執行步驟DS88；如果不符合最粗糙尺寸定義，則執行步驟DS83。
步驟DS83對行列數為N的源函數項矩陣中每一個地表單元的源函數項和行列數為N的高程猜測值矩陣中每一個地表單元的初始高程猜測值進行殘留誤差限制計算，得到行列數為N/2的源函數項矩陣中每一個地表單元的源函數項ρN/2←R(ρN-A·HN)。
步驟DS84對行列數為N/2的源函數項矩陣的第一列ρN/2(x，0)和第一行ρN/2(0，y)中每一個地表單元的源函數項進行邊界校正 ρN/2(0，y)＝-ρN/2(1，y)，ρN/2(x，0)＝-ρN/2(x，1)。
步驟DS85將行列數為N/2的高程猜測值矩陣中每一個地表單元的初始高程猜測值設為0HN/2←0。
步驟DS86將行列數為N/2的高程猜測值矩陣中的每一個地表單元的初始高程猜測值及源函數項代入步驟DS81開始進行遞歸計算，更新行列數為N/2的高程猜測值矩陣中每一個地表單元的高程猜測值N←N/2，HN/2←MV(HN/2，ρN/2)。
步驟DS87對行列數為N/2的高程猜測值矩陣中的每一個地表單元的高程猜測值使用延拓算子得到行列數為N的每一個地表單元的高程延拓值，與行列數為N的每一個地表單元的初始高程猜測值相加，更新行列數為N的高程猜測值矩陣中每一個地表單元的高程猜測值HN←HN+P(HN/2)。
步驟DS88對行列數為N的高程猜測值矩陣HN進行n次鬆弛操作，n為任意自然數，在本實例中取值為9，且可以與步驟DS81中的次數不同。
步驟E對地表單元的方位向解析度、地距向解析度以及相對高程猜測值做二維差分運算，計算每一個地表單元的方位向坡度修正值ω(x，y)和地距向坡度修正值β(x，y)。
方位向坡度修正值 地距向坡度修正值 對方位向坡度矩陣ω進行邊界校正 ω(0，y)＝-ω(1，y)(23) 對地距向坡度矩陣β進行邊界校正 β(0，y)＝-β(1，y)(24) 步驟F請參見圖3d是步驟F的子流程圖，對地表單元的方位向解析度及地距向解析度、方位向坡度修正值和地距向坡度修正值使用三次完全多重網格算法，結合一個外部地表單元參考點HT(r，c)計算每一個地表單元的絕對高程估計值H(x，y)。
H(x，y)＝HT(r，c)-Hg(r，c)+Hg(x，y)(25) 本發明所述的地形測量方法中，圓極化算法利用了地形—極化方位角偏移模型，當地表場景中以裸露土地和低矮植被為主時，該模型的有效範圍限定於L波段(波長為19.3—76.9釐米)和P波段(波長為76.9—133釐米)的正側視條帶成像模式。當地表場景中以茂密的森林為主時，該模型的有效範圍限定於P波段(波長為76.9—133釐米)的正側視條帶成像模式。
機載單路徑全極化三通道復圖像數據如圖6所示，圖6a每一個像素點的灰度值GHH(x，y)是由HH通道復圖像數據經過下式處理得到的 GHH(x，y)＝10·log10(|SHH(x，y)|2)(26) 圖6b每一個像素點的灰度值GHV(x，y)是由HV通道復圖像數據經過下式處理得到的 GHV(x，y)＝10·log10(|SHV(x，y)|2)(27) 圖6c每一個像素點的灰度值GVV(x，y)是由HV通道復圖像數據經過下式處理得到的 GVV(x，y)＝10·log10(|SVV(x，y)|2)(28) 圖像尺寸為3797行，1024列，已知雷達工作波長為21釐米(L波段)，載機平臺平均高度HA為7681米，近地端入射角η0為26.2度，方位向解析度Ra為10米，斜距解析度RS為10米，最粗糙尺寸定義為NC＝4。以該數據為例，遵照本發明所述的步驟，圖7a為本發明輸出的圖6a中方框區域內的極化方位角偏移值圖，每個地表單元的極化方位角偏移值用每個像素點處的灰度值來表示；圖7b為地表真實極化方位角偏移值數據，可以檢驗本發明實例輸出結果的正確程度。將圖7a中的兩條白色線段覆蓋的像素點處的灰度值，即本發明輸出的極化方位角偏移值用「*」顯示在圖7c(剖面曲線1)和圖7d(剖面曲線2)中。相同圖像坐標位置處的地表真實極化方位角偏移值用「—」顯示在圖7c(剖面曲線1)和圖7d(剖面曲線2)中。圖7c和圖7d的縱坐標軸表示極化方位角偏移值，單位是度；圖7c的橫坐標表示圖7a中的列坐標，圖7d的橫坐標表示圖7a中的行坐標。圖7c和圖7d中的「*」線與「—」線的變化趨勢十分相似，表明本發明得到的極化方位角偏移值很接近地表實際測量結果。本發明得到的極化方位角偏移值與地表實際測量結果之間的均方根誤差為9.3度。圖8a和圖8c分別為本發明輸出的圖6a中方框區域內的方位向坡度圖和地距向坡度圖，每個地表單元的方位向和地距向坡度值用每個像素點處的灰度值來表示；圖8b和圖8d分別為地表真實方位向坡度和地距向坡度數據，可以檢驗本發明實例輸出結果的正確程度。將圖8a和圖8c中的白色線段覆蓋的像素點處的灰度值，即本發明輸出的方位向坡度值和地距向坡度值用「*」顯示在圖9a和圖9b中。相同圖像坐標位置處的地表真實方位向坡度值和地距向坡度值用「—」顯示在圖9a和圖9b中。圖9a和圖9b的縱坐標軸分別表示方位向坡度值和地距向坡度值，單位是度；圖9a的橫坐標表示圖8a中的行坐標，圖9b的橫坐標表示圖8c中的列坐標。圖9a和圖9b中的「*」線與「—」線的變化趨勢十分相似，表明本發明得到的方位向坡度值和地距向坡度值很接近地表實際測量結果。本發明得到的方位向坡度值和地距向坡度值與地表實際測量結果之間的均方根誤差分別為6.2度和8.5度。圖10b為本發明輸出的絕對高程圖，每個地表單元的絕對高程值用每個像素點處的灰度值來表示；圖10a為地表真實絕對高程數據，可以檢驗本發明實例輸出結果的正確程度。圖11b為本發明輸出的圖10a中方框區域內的絕對高程圖，每個地表單元的絕對高程值用每個像素點處的灰度值來表示；圖11a為地表真實絕對高程數據，可以檢驗本發明實例輸出結果的正確程度。將圖11a中的四條白色線段覆蓋的像素點處的灰度值，即地表真實絕對高程值用細線顯示在圖12a、b、c、d(剖面曲線1、2、3、4)中。相同圖像坐標位置處的本發明輸出的絕對高程值用粗線顯示在圖12a、b、c、d(剖面曲線1、2、3、4)中。圖12a、b、c、d的縱坐標軸表示絕對高程值，單位是米；圖12a、b的橫坐標表示圖11a中的列坐標，圖12c、d的橫坐標表示圖11a中的行坐標。圖12a、b、c、d中的粗線與細線的變化趨勢十分相似，表明本發明得到的絕對高程值很接近地表實際測量結果。本發明得到的絕對高程值與地表實際測量結果之間的均方根誤差為62米。
以上所述，僅為本發明中的具體實施方式
，但本發明的保護範圍並不局限於此，任何熟悉該技術的人在本發明所揭露的技術範圍內，可理解想到的變換或替換，都應涵蓋在本發明的包含範圍之內，因此，本發明的保護範圍應該以權利要求書的保護範圍為準。
權利要求
1、一種地形測量方法，其特徵在於
步驟A利用單路徑全極化合成孔徑雷達復圖像數據並採用圓極化算法計算，得到照射區域內每一個地表單元的極化方位角偏移中間值為
並進行值域映射為得到照射區域內地表單元的極化方位角偏移值θ(x，y)，式中SHH、SHV、SVH和SVV為四幅全極化復圖像數據，每一幅圖像為一個二維矩陣，圖像中每一個像素點在矩陣中的位置為(x，y)，表示圖像中每一個地表單元的行、列坐標，從0開始計數；
步驟B利用雷達載體平臺高度HA、近地端入射角度η0及斜距解析度RSL來計算每一個地表單元的雷達入射角度η(x，y)、方位向解析度Ra(x，y)和地距解析度Rg(x，y)；
步驟C對地表單元的極化方位角偏移值和雷達入射角使用補償朗伯算法，計算每一個地表單元的方位向坡度猜測值ωg(x，y)和地距向坡度猜測值βg(x，y)；
步驟D對地表單元的方位向解析度及地距向解析度、方位向坡度猜測值和地距向坡度猜測值使用完全多重網格算法，計算每一個地表單元的相對高程猜測值Hg(x，y)；
步驟E對地表單元的方位向解析度、地距向解析度以及相對高程猜測值做二維差分運算，計算每一個地表單元的方位向坡度修正值ω(x，y)和地距向坡度修正值β(x，y)；
步驟F對地表單元的方位向解析度及地距向解析度、方位向坡度修正值和地距向坡度修正值使用完全多重網格算法，結合一個外部地表單元參考點HT(r，c)計算每一個地表單元的絕對高程估計值H(x，y)。
2、根據權利要求1中所述的地形測量方法，其特徵在於，其中圓極化算法利用了地形—極化方位角偏移模型，該模型的有效範圍限定於L波段和P波段的正側視條帶成像模式。
3、根據權利要求1中所述的地形測量方法，其特徵在於，其中補償朗伯算法定義符合極化補償處理後的地形—極化方位角偏移模型和朗伯反射模型為
地形—極化方位角偏移模型為
朗伯反射模型為
其中I(x，y)代表極化合稱圖像中每一個地表單元的亮度值，ω(x，y)代表每一個地表單元的方位向坡度，β(x，y)代表每一個地表單元的地距向坡度，K為雷達系統的定標常數，σ0(x，y)為每一個地表單元的歸一化後向散射係數，Rg(x，y)為每一個地表單元的地距向解析度，Ra(x，y)為每一個地表單元的方位向解析度，η(x，y)為每一個地表單元的入射角度，θ(x，y)為每一個地表單元的極化方位角偏移值。
4、根據權利要求1中所述的地形測量方法，其特徵在於，其中補償朗伯算法的計算步驟為
步驟CS1對四幅全極化復圖像數據每一個地表單元的SHH(x，y)、SHV(x，y)、SVH(x，y)和SVV(x，y)進行組合，得到每一個地表單元的斯託克斯矩陣元素m22(x，y)，m23(x，y)和m33(x，y)；
步驟CS2對每一個地表單元的斯託克斯矩陣元素m22(x，y)，m23(x，y)和m33(x，y)進行極化補償後，獲得補償後的極化合成圖像每一個地表單元的亮度為
式中補償後的每一個地表單元的亮度值為I′(x，y)，補償後的每一個地表單元的極化方位角偏移值θ(x，y)為0；在地形—極化方位角偏移模型中等價於每一個地表單元的方位向坡度值ω(x，y)為0；極化補償等價於令朗伯反射模型變為
步驟CS3保持原有極化方位角偏移值不變的極化合成圖像每一個地表單元的亮度為I(x，y)＝m22(x，y)；
步驟CS4利用補償及未補償的極化合成圖像每一個地表單元的亮度進行每一個地表單元的方位向坡度猜測值ωg(x，y)的計算為令每一個地表單元的方位向坡度猜測值ωg(x，y)的符號等於相同地表單元的極化方位角偏移值θ(x，y)的符號；
步驟CS5利用每一個地表單元的極化方位角偏移值及方位向坡度猜測值進行每一個地表單元的地距向坡度猜測值βg(x，y)計算
5、根據權利要求1中所述的地形測量方法，其特徵在於，其中完全多重網格算法計算源函數項的步驟為
步驟DP1利用地表單元的方位向解析度及地距向解析度、方位向坡度修正值和地距向坡度修正值計算每一個地表單元的方位向高程差Δa(x，y)＝Ra(x，y)·tan(ωg(x，y))和地距向高程差Δg(x，y)＝Rg(x，y)·tan(βg(x，y))；
步驟DP2利用每一個地表單元的方位向高程差和地距向高程差計算每一個地表單元的源函數項
ρ(x，y)＝Δa(x+1，y)-Δa(x，y)+Δg(x，y+1)-Δg(x，y)；
執行邊界校正令源函數項矩陣的第一行元素為對應同一列上第二行元素的相反數；令源函數項矩陣的第一列元素為對應同一行上第二列元素的相反數。
6、根據權利要求1中所述的地形測量方法，其特徵在於，其中完全多重網格算法計算一個高程值估計周期的步驟為
步驟DS1對高程值矩陣HN的行列數N進行判斷如果符合最粗糙尺寸定義N≤NC，則執行步驟DS8；如果不符合最粗糙尺寸定義，則執行步驟DS2；
步驟DS2對行列數為N的源函數項矩陣中每一個地表單元的源函數項和行列數為N的高程猜測值矩陣中每一個地表單元的初始高程猜測值進行殘留誤差限制計算，得到行列數為N/2的源函數項矩陣中每一個地表單元的源函數項ρN/2←R(ρN-A·HN)；
步驟DS3對行列數為N/2的源函數項矩陣中第一行和第一列中的每一個地表單元的源函數項進行邊界校正ρN/2(0，y)＝-ρN/2(1，y)，ρN/2(x，0)＝-ρN/2(x，1)；
步驟DS4對行列數為N的高程猜測值矩陣中每一個地表單元的初始高程猜測值進行限制操作，得到行列數為N/2的高程猜測值矩陣中每一個地表單元的初始高程猜測值HN/2←R(HN)；
步驟DS5對行列數為N/2的高程猜測值矩陣中第一行和第一列中每一個地表單元的初始高程猜測值進行邊界校正HN/2(0，y)＝-HN/2(1，y)，HN/2(x，0)＝-HN，2(x，1)；
步驟DS6將行列數為N/2的高程猜測值矩陣中的每一個地表單元的初始高程猜測值以及源函數項代入步驟DS1進行遞歸計算，更新行列數為N/2的高程猜測值矩陣中每一個地表單元的高程猜測值HN/2←FMG(HN/2，ρN/2)；
步驟DS7對行列數為N/2的高程猜測值矩陣中的每一個地表單元的高程猜測值使用延拓算子得到行列數為N的每一個地表單元的高程延拓值，與行列數為N的每一個地表單元的初始高程猜測值相加，更新行列數為N的高程猜測值矩陣中每一個地表單元的高程猜測值HN←HN+P(HN/2)；
步驟DS8對行列數為N的高程猜測值矩陣中的每一個地表單元的高程猜測值以及行列數為N的源函數項矩陣中的每一個地表單元的源函數項使用V形多重網格算法，更新行列數為N的高程猜測值矩陣中每一個地表單元的高程猜測值HN←MV(HN，ρN)。
7、根據權利要求6中所述的地形測量方法，其特徵在於，其中V形多重網格算法的計算步驟為
步驟DS81對行列數為N的高程猜測值矩陣HN進行n次鬆弛操作，n為任意自然數，一次鬆弛運算表示為
HN(x，y)＝[(HN(x+1，y)+HN(x-1，y)+HN(x，y+1)+HN(x，y-1))-ρN(x，y)]/4；
步驟DS82對矩陣行列數N進行判斷如果符合最粗糙尺寸定義，則執行步驟DS88；如果不符合最粗糙尺寸定義，則執行步驟DS83；
步驟DS83對行列數為N的源函數項矩陣中每一個地表單元的源函數項和行列數為N的高程猜測值矩陣中每一個地表單元的初始高程猜測值進行殘留誤差限制計算，得到行列數為N/2的源函數項矩陣中每一個地表單元的源函數項ρN/2←R(ρN-A·HN)；
步驟DS84對行列數為N/2的源函數項矩陣的第一列ρN/2(x，0)和第一行ρN/2(0，y)中每一個地表單元的源函數項進行邊界校正為
ρN/2(0，y)＝-ρN/2(1，y)，ρN/2(x，0)＝-ρN/2(x，1)；
步驟DS85將行列數為N/2的高程猜測值矩陣中每一個地表單元的初始高程猜測值設為0HN/2←0；
步驟DS86將行列數為N/2的高程猜測值矩陣中的每一個地表單元的初始高程猜測值及源函數項代入步驟DS81開始進行遞歸計算，更新行列數為N/2的高程猜測值矩陣中每一個地表單元的高程猜測值HN/2←MV(HN/2，ρN/2)；
步驟DS87對行列數為N/2的高程猜測值矩陣中的每一個地表單元的高程猜測值使用延拓算子得到行列數為N的每一個地表單元的高程延拓值，與行列數為N的每一個地表單元的初始高程猜測值相加，更新行列數為N的高程猜測值矩陣中每一個地表單元的高程猜測值HN←HN+P(HN/2)；
步驟DS88對行列數為N的高程猜測值矩陣HN進行n次鬆弛操作，n為任意自然數，且與步驟DS81中的次數相同或不同。
8、根據權利要求1中所述的地形測量方法，其特徵在於，其中完全多重網格算法需要進行NF次循環，NF的典型值為3，最終得到每一個地表單元的相對高程猜測值Hg(x，y)。
9、根據權利要求1中所述的地形測量方法，其特徵在於，其中每一個地表單元的方位向坡度修正值為
每一個地表單元的地距向坡度修正值為
令方位向坡度矩陣ω的第一行和第一列元素等於其第二行和第二列元素的相反數；令地距向坡度矩陣β的第一行和第一列元素等於其第二行和第二列元素的相反數。
10、根據權利要求1中所述的地形測量方法，其特徵在於，其中每一個地表單元的絕對高程猜測值為H(x，y)＝HT(r，c)-Hg(r，c)+Hg(x，y)；(r，c)為外部參考點在高程猜測值矩陣中的行列坐標，Hg(r，c)為外部參考點處的相對高程猜測值，HT(r，c)為外部參考點處的絕對高程值。
全文摘要
本發明涉及一種地形測量方法，利用單路徑全極化合成孔徑雷達復圖像數據並採用圓極化算法計算得到照射區域內每一個地表單元的極化方位角偏移值；利用系統參數計算每一個地表單元的入射角度和解析度；使用補償朗伯算法計算每一個地表單元的坡度猜測值；利用完全多重網格算法計算每一個地表單元的絕對高程值。利用單路徑的全極化合成孔徑雷達斜距復圖像數據快速地獲得中等精度坡度和數字高程模型，本發明方法的特點在於對雷達數據獲取的限制較少、對外部數據的依賴較少、自動化程度較高以及整體流程計算效率較高，從而使其在地形可視化應用領域具有較大的應用價值。
文檔編號G01S13/00GK101482616SQ20081011832
公開日2009年7月15日 申請日期2008年8月13日 優先權日2008年8月13日
發明者洋 李, 文 洪, 芳 曹, 王彥平, 吳一戎 申請人:中國科學院電子學研究所