單次飛行全極化合成孔徑雷達圖像反演地面數字高程的方法

專利名稱：單次飛行全極化合成孔徑雷達圖像反演地面數字高程的方法
技術領域：
本發明屬于空間遙感與通信技術領域，具體涉及一種星載合成孔徑雷達反演地面數字高程的方法。
但是INSAR各相關技術難度高，并且不是每一處地形都能容易得到符合要求的相干的SAR圖像。
用全極化SAR圖像對地表分類已經得到了廣泛的研究。從全極化SAR圖像數據，可以得到2×2維的復散射振幅函數Fpq(p，q＝v，h)和4×4維的實Mueller矩陣Mij(i，j＝1，...，4)([1]Jin Y.Q.，Electromagnetic scattering modelling forquantitative remote sensing，1994，(SingaporeWorld Scientific))。同極化和交叉極化后向散射強度是入射波極化橢圓角χ與取向角ψ的函數。當地面平坦時，后向散射信號一般在入射極化取向角ψ＝0處達到極大值。但是，當地面有傾斜坡度時，同極化或交叉極化后向散射信號的極大值會從ψ＝0處發生遷移。這種ψ的遷移可用來推算地表面坡度([2]Lee J.S.，Jansen，R.W.，Schuler，D.L.，Polarimetric analysis andmodeling of multi-frequency SAR signatures from Gulf stream fronts，IEEE.Journalof Oceanic Engineering，1998，23322-333)。用實同極化和零交叉極化的散射振幅函數的假設條件推導了用散射振幅函數表示的ψ的遷移。由于要同時確定水平方位角和射程角，所以需要兩次或多次相近飛行的SAR圖像數據來確定表面坡度。近年來，用兩次相近飛行的SAR和INSAR數據來反演地面數字高程(DEM)的方法已有很好的實例(見文獻[3])。
本發明提出的用單次飛行的極化SAR數據圖像實現DEM反演的方法，是用全極化散射Mueller矩陣解，推導ψ的遷移，將其表示為三個散射Stokes參數Ivs，Ihs，Us的函數，而Ivs，Ihs，Us是由SAR測量得到的；利用主坐標系和像素元的局部坐標系之間的歐拉(Euler)角變換，把入射波極化取向角ψ與射程角β、水平方位角γ以及SAR觀測的幾何構造直接聯系起來；利用傾斜地面的水平方位排列在SAR圖像中表現出來的線性紋理，用圖像形態學細化算法確定圖像中各像素元的地形水平方位角γ，再用Ivs，Ihs，Us計算ψ遷移來確定射程角β；最后用β和γ來確定各像素元地形水平方位坡度與射程方位坡度；最后用完整多重網格算法數值求解地面高程阿松(Poisson)方程，完成SAR觀測區域地表DEM的反演。
本發明的內容進一步描述如下1、用散射Stokes參數表達取向角ψ遷移考慮一極化電磁波以角度(π-θi，i)入射到地面上，其極化可用橢圓角χ和取向角ψ來表征(見文獻[1])。散射場可寫為 其中2×2維的復散射振幅函數Fpq(p，q＝v，h)可由全極化散射測量技術得到，它包含了散射層的體散射和面散射。由式(1)，可推算得到Stokes矢量的實Mueller矩陣解(見文獻[1])， 其中歸一化的入射Stokes矢量可寫成I1(χ，ψ)＝[Ivi，Ihi，Ui，Vi]T＝
T(3)Mueller矩陣元素Mij(i，j＝1，Λ，4)由 (p，q，s，t＝v，h)的實部或虛部的系綜平均表示，表達式可見文獻[1]。應注意到，同極化項 總是比交叉極化項 (s≠t)要大得多。
同極化的后向散射系數σc(θi，π+i；π-θi，i)可表示為(見文獻[1])，σc＝4πcosiPn(4)其中Pn＝0.5[Ivs(1-cos2χcos2ψ)+Ihs(1+cos2χcos2ψ)+Uscos2χsin2ψ+Vssin2χ](5)當地面平坦時，同極化后向散射σc(χ，ψ)的極大值一般總在ψ＝0處(見文獻[2])。然而，研究表明，當地面傾斜時，σc的極大值會從ψ＝0處遷移([3]Schuler，D.L，Lee，J.S.，Answorth，T.L.，and Grunes，M.R.，Terrain topography measurements using multipasspolarimetric synthetic aperture radar data，Radio Science，2000，35813-832；[4]Lee，J.S.，Schuler，D.L.and Answorth，T.L.，Polarimetric SAR data compressionfor terrain azimuthal slope variation，IEEE Transaction on Geoscience RemoteSensing，2000，38(5)2153-2163；[5]Schuler，D.L.，Answorth，T.L.，and Lee J.S.，Topographic mapping using polarimetric SAR data，Int.J.Remote Sensing，1998，19(1)141-160)。
在Pn/ψ＝0處，σc取極大值。假設χ＝0的對稱情形，由(5)式可得到0＝-(Ihs-Ivs)sin2ψ+Uscosψ+0.5(Ihs+Ivs)′+0.5U′ssin 2ψ+0.5(Ihs-Ivs)′cos2ψ(6)其中上標′表示/ψ。由(2)式的 比較(6)式右端各項，有0.5(Ihs+Ivs)'~(M13+M23)cos2ψ~(Re<FvvFvh*>+Rc<FhvFhh*>)cos2ψ---(7a)]]>0.5Us'~M33cos2ψ~Re<FvvFhh*+FvhFhv*>cos2ψ---(7b)]]>(Ihs-Ivs)～0.5(M11+M22)cos 2ψ～0.5(<|Fvv|2>+<|Fhh|2>)cos 2ψ (7c)(Ihs-Ivs)'~0.5(M13-M23)cos2ψ~0.5(Re<FvvFvh*>-Re<FhvFhh*>)cos2ψ---(8a)]]>Us~0.5(M31+M32)~Re<FvvFhv*>+Re<FvhFhh*---(8b)]]>可看到，(7a)與(7b)均遠小于(7c)，而(8a)亦遠小于(8b)。因此(6)式右端的最后三項可略去。于是，本發明中得到σc的極大值的ψ遷移取為tan2ψ=UsIhs-Ivs---(9)]]>由此可以看出，第3個Stokes參數Us≠0引起了ψ的遷移。
順便說，如果Us和Ihs-Ivs都趨向于0，比如，當散射來自于均勻取向散射體或各向同性散射介質，如厚的植被頂層時，式(9)就不能很好地確定ψ的遷移。
進一步地，從式(2)可得到Us~0.5(M31+M32)=Re(<FvvFhv*>+<FvhFhh*>)---(10a)]]>Ihs-Ivs～0.5(M22-M11)＝0.5(<|Fhh|2>-<|Fvv|2>) (10b)2、用Euler角變換得到射程角和水平方位角入射波的極化矢量 在主坐標系中定義為(見文獻[1])，h^i=z^×k^i|z^×k^i|=-sinφix^+cosφiy^,v^i=h^i×k^i---(12)]]>其中入射波矢量為k^i=sinθicosφix^+sinθisinφiy^-cosθiz^---(13)]]>當該像素地表面有如

圖1所示的局部坡度時，由局部法向量 定義的局部坐標下的極化矢量可表示為h^b=z^b×k^ib|z^b×k^ib|,v^b=h^b×k^ib---(14)]]>通過兩個坐標系 和 之間的Euler角(β，γ)變換(見文獻(1))，有x^=cosγcosβx^b+sinγy^b-cosγβz^b---(15a)]]>y^=-sinγcosβx^b+cosγy^b+sinγsinβz^b---(15b)]]>z^=sinβx^b+cosβz^b---(15c)]]>將以上三式代入式(12)，(13)中，可得h^i=(-cosγcosβsinφi-sinγcosβcosφi)x^b+(-sinγsinφi+cosγcosφi)y^b---(16)]]>+(cosγsinβsinφi+sinγsinβcosφi)z^b]]>k^ib=(cosγcosβsinθicosφi-sinγcosβsinθisinφi-sinβcosθi)x^b]]>+(sinγsinθicosφi+cosγsinθisinφi)y^b---(17)]]>+(-cosγsinβsinθicosφi+sinγsinβsinθisinφi-cosβcosθi)z^b]]>將式(17)代入式(14)中，該像素地表面的極化矢量 可表示為h^b=ax^b+by^ba2+b2---(18)]]>其中a＝-(sinγsinθicosφi+cosγsinθisinφi)(19)b＝(cosγcosβsinθicosφi-sinγcosβsinθisinφi-sinβcosθi)于是，由式(16)，(18)可以得到入射波極化的取向角ψ為cosψ=h^b·h^i=cosβsinθi-cos(γ+φi)sinβcosθia2+b2(20a)]]>tanψ=tanβsin(γ+φi)sinθi-cosθitanβcos(γ+φi)---(20b)]]>通常令入射φi＝0， 是射程方向， 是水平方位方向，注意在文獻[2]中，水平方位角定義為 和 之間的夾角。于是，該像素地表面的射程方位坡度SR和水平方位坡度SA可定義為(見文獻[2])SR＝tanβcosγ(21)SA＝tanβsinγ因為在式(20b)中，用一個ψ不能同時確定β和γ兩個角度，所以通常需要兩次或多次相近飛行的SAR圖像數據，用式(20b)建立兩個方程來求解β和γ，這在文獻[2]中已有討論。
3、確定SAR圖像中各像素的水平方位角如果只有單次飛行的SAR圖像數據，除了用ψ與Ivs，Ihs，Us的關系式(9)外，還要找到另一個已知條件，用來確定兩個未知的β和γ。在SAR圖像中，傾斜地表面水平方位上的排列表現出的排列紋理實際上是傾斜地表排列水平取向方位的標志。本發明采用自適應閾值化方法對SAR圖像作二值化和形態學處理([6]Castleman，K.R.，Digital imageprocessing，1996(Printice Hall))，獲得每個像素的水平方位角γ，具體步驟如下(a)濾波對圖像進行斑點濾波，去除孤立的噪聲斑點，初步提高圖像信噪比。
(b)二值化用自適應閾值化的方法([7]Francis H.Y.Chan，F.K.Lam and Hui Zhu，Adaptive Thresholding by Variational Method，IEEE Transaction on Geoscience RemoteSensing，1998，7(3)468-473)，對斑點濾波后的SAR圖像進行二值化處理。由于SAR圖像背景不均勻，不適宜采用全局固定的閾值。
(c)對得到的二值圖像作形態學處理，消除孤立的像素，填補微小的孔隙，我們把二值圖像中為“1”的部分稱為前景，為“0”的部分稱為背景。
(d)對前景部分，提取邊緣以邊緣上每個像素為中心像素，開一個方形窗口，得到通過中心像素的曲線段。然后用最小二乘法作曲線的多項式擬合，確定邊緣上像素處的切線斜率，從而得到該像素的水平方位角。對已計算的像素作上標記，下一次不再計算。在計算中，我們取方形窗口的邊長為21像素。
(e)從前景部分去掉(d)中得到的邊緣，得到新的前景。重復(d)，直到獲得前景部分的每個像素的水平方位角。
(f)將步驟(c)得到的二值圖像取反，原二值圖像的背景成為前景，重復(d)和(e)步驟，直到獲得背景部分每個像素的水平方位角。
在只有單次飛行的SAR數據的情況下，以上的方法為確定水平方位角γ提供一種可靠的手段。
在得到各個像素元的水平方位角γ后，用(9)式確定的ψ和(20b)式得到各個像素元的β角為tanβ=tanψsinθisinγ+tanψcosθicosγ(22)]]>其中θi由SAR觀測的幾何結構給出。由β和γ得到式(21)中的SA和SR，在此基礎上采用文獻[11-12]的方法來反演地形DEM。
4、計算由坡度函數S(x，y)獲得的曲率值ρ(x，y)。生成DEM是尋求在M×N的矩形網格區域上的離散Poisson方程的解，類似于INSAR的相位展開。([8]H.Takajo and T.Takahashi，Least-squares phase estimation from phase difference，J.Opt.Soc.Amer.A，1998，5416-425；[9]D.C.Ghiglia and L.A.Romero，Robusttwo-dimensional weighted and unweighted phase unwrapping that uses fast transforms anditerative methods，J.Opt.Soc.Amer.A，1994，11(1)107-117；[10]M.D.Pritt and J.S.Shipman，Least-squares two-dimensional phase unwrapping using FFT’s，IEEE Transaction on GeoscienceRemote Sensing，1994，32(3)706-708)。Poisson方程為2φ(x，y)＝ρ(x，y) (23)其中2是Laplace算子2/x2+2/y2，φ(x，y)表示在(x，y)處的高度值，并令初始值為零。方程右端的源函數ρ(x，y)是由輸入坡度函數S(x，y)得到的表面曲率值。我們假定x方向為射程方位，而y方向為水平方位，則dSR(x)/dx和dSA(y)/dy分別包含了射程方位和水平方位上的表面曲率(上標A和R分別表示水平方位和射程方位)。
水平方位和射程方位的源函數的離散值由下式給出，ρijR=SijR-Si-1,jR,ρijA=SijA-Si,j-1A(0≤i≤M,0≤j≤N)---(24)]]>假定輸入坡度的誤差滿足高斯(Gauss)分布，地面的高程可看成坡度的積分，那么地面高程的誤差也滿足Gauss分布。因此，一個服從χ2分布的統計量可看成期望解與實際值相符合的量度標準(見文獻[5])。這也是干涉SAR中相位展開算法的標準假設。χ2的值由下式給出x2=Σi,j(φij-φi-1,j-SijR)2+Σi,j(φij-φi,j-1-SijA)2---(25)]]>求解φij的偏微商值與式(24)所定義的偏微商值的差必須為最小。對χ2求導取零得到(φi+1，j-2φij+φi-1，j)+(φi，j+1-2φij+φi，j-1)＝ρij(26)式(26)就是需要求解的離散Poisson方程。
離散的源函數ρij由下式給出ρij=(Si+1,jR-SijR)+(Si,j+1A-SijA)---(27)]]>文獻[2，5]中用于生成全極化SAR的DEM的算法是交替方向隱式迭代(AlternatingDirections Implicit，ADI)算法，本發明使用完整多重網格(Full Multigrid，FMG)算法。此算法收斂快，魯棒性強，計算量小。([11]M.D.Pritt，Phase unwrapping bymeans of multi-grid techniques for interferometric SAR，IEEE Transaction onGeoscience Remote Sensing，1996，34728-738；[12]J.Demmel，Applied NumericalLinear Algebra，1997，Society for Industrial and Applied Mathematics)對于M×N的矩形網格，其計算量為(MN)ln(MN)。FMG算法的具體步驟可參見[11]。計算框圖見圖13所示。
圖2為中心像素的水平方位角。
圖3為中國廣東省惠州的SIR-CSAR圖像(L波段vv，hh，vh總功率)圖4為惠州方位坡度圖，其中圖4(a)為水平方位坡度，圖4(b)為射程方位坡度。
圖5為惠州由SAR圖像反演出的地形DEM。
圖6為惠州反演出的地形高度偽彩色圖。
圖7為惠州的等高線圖。
圖8為臺灣阿里山地區的STR-CSAR圖像(L波段vv，hh，vh總功率)。
圖9為阿里山地區的方位坡度圖，其中圖9(a)為水平方位坡度，圖9(b)為射程方位坡度。
圖10為阿里山地區由SAR圖像反演出的地形DEM。
圖11為阿里山地區反演出的地形高度偽彩色圖。
圖12為阿里山地區的等高線圖。
圖13為本發明計算程序框圖。
采用美國加州理工學院噴氣推進實驗室(JPL)提供的1994年10月4日SAR SIR-C對中國廣東省觀測的L波段數據，圖像中心位置約在(114.561°E，22.815°N)，每個像素的大小是12.5×12.5m2。SAR在圖像中心的入射角是32.85°。圖3是從原圖中截取的大小為512×512像素的后向散射vv，hh，hv總功率圖像，其中心位置為(114.7826°E，23.1359°N)，位于中國廣東省的惠州地區。
從圖3可以看出南北走向的丘陵起伏。為了突出地形的主要走向，我們略去了一些更為細微的局部起伏變化。注意，這里定義的水平方位角γ是坡度走向 與方位 的夾角(如圖1(a)所示)。
SA和SR的計算結果如圖4(a)和4(b)所示，圖5是由SAR圖像反演出的DEM，圖6是反映地形高度的偽彩色圖，圖7是地形的等高線圖。
實施例2，臺灣阿里山地區采用美國加州理工學院噴氣推進實驗室(JPL)提供的1994年10月4日SAR SIR-C對中國臺灣觀測的L波段數據，圖像中心位置約在(120.964°E，23.487°N)，每個像素的大小是12.5×12.5m2。SAR在圖像中心的入射角是32.85°。圖8是從原圖中截取的大小為512×512像素的后向散射vv，hh，hv總功率圖像，其中心位置為(120.8626°E，23.6659°N)，位于中國臺灣的阿里山地區。
反演的DEM與原SAR圖像所反映的直覺的地形結構以及十分簡略的地面高程圖[13]是基本一致的，其主要地形的走向與SAR圖像的明暗紋理也是吻合的。
權利要求
1.一種實現地面數字高程反演方法，其特征在于利用全極化散射Mueller矩陣解，推導ψ的遷移，將其表示為三個散射Stokes參數Ivs，Ihs，Us的函數，而Ivs，Ihs，Us是由SAR測量得到的；利用主坐標系和像素元的局部坐標系之間的歐拉角變換，把入射波極化取向角ψ與射程角β、水平方位角γ以及SAR觀測的幾何構造直接聯系起來；利用傾斜地面的水平方位排列在SAR圖像中表現出來的線性紋理，用圖像形態學細化算法確定圖像中各像素元的地形水平方位角γ，再用Ivs，Ihs，Us計算ψ遷移來確定射程角β；最后用β和γ來確定各像素元地形水平方位坡度與射程方位坡度；最后用完整多重網格算法數值求解地面高程泊松方程，完成SAR觀測區域地表DEM的反演。
2.根據權利要求1所述的反演方法，其特征在于同極化的后向散射系數σc(θi，π+i；π-θi，i)的極大值的ψ遷移取為tan2ψ=UsIhs-Ivs---(9)]]>其中Ivs，Ihs，Us為三個散射Stokes參數Us＝-cos2xsin2ψ，Ihs＝0.5(1+cos2xcos2ψ)，Ivs＝0.5(1-cos2xcos2ψ)；經過Euler變換，得到入射波極化的取向角ψ為tanψ=tanβsin(γ+φi)sinθi-cosθitanβcos(γ+φi)---(20b)]]>
3.根據權利要求1所述的反演方法，其特征在于采用自適應閾值化方法對SAR圖像作二值化和形態學處理，獲得每個像素的水平方位角γ；然后由(9)式和(20b)式得到各個像素元的β角為tanβ=tanψsinθisinγ+tanψcosθicosγ---(22)]]>其中θi由SAR觀察的幾何結構給出，于是得到SR＝tanβcosγSA＝tanβsinγ然后反演地形DEM。
4.根據權利要求3所述的反演方法，其特征在于對SAR圖像作二值化和形態學處理的具體步驟如下(a)濾波對圖像進行斑點濾波，去除孤立的噪聲斑點，初步提高圖像信噪比；(b)二值化用自適應閾值化的方法，對斑點濾波后的SAR圖像進行二值化處理；(c)對得到的二值圖像作形態學處理，消除孤立的像素，填補微小的孔隙，我們把二值圖像中為“1”的部分稱為前景，為“0”的部分稱為背景；(d)對前景部分，提取邊緣以邊緣上每個像素為中心像素，開一個方形窗口，得到通過中心像素的曲線段，然后用最小二乘法作曲線的多項式擬合，確定邊緣上像素處的切線斜率，從而得到該像素的水平方位角；(e)從前景部分去掉(d)中得到的邊緣，得到新的前景。重復(d)，直到獲得前景部分的每個像素的水平方位角；(f)將步驟(c)得到的二值圖像取反，原二值圖像的背景成為前景，重復(d)和(e)步驟，直到獲得背景部分每個像素的水平方位角。
5.根據權利要求3所述的反演方法，其特征在于采用完整多重網格算法求解如下的離散Poisson方程，得到地面數字高程(φi+1，j-2φij+φi-1，j)+(φi+1-2φij+φi，j-1)＝ρij(26)其中φij(x，y)表示在(x，y)處的高度值，pij(x，y)是由輸入坡度函數s(x，y)得到的表面曲率值。
全文摘要
本發明是單次飛行全極化合成孔徑雷達圖像反演地面數字高程的方法。利用全極化散射Mueller矩陣解，推導取向角ψ的遷移，借助Euler角的變換，建立了由3個Stokes參數確定的ψ遷移與傾斜表面像素的射程角β和水平方位角γ之間的聯系。在只有單次飛行SAR數據的條件下，用自適應閾值的二值化和圖像形態學細化算法來確定γ以及β，以及地表面射程方位坡度和水平方位坡度。由此用完整多重網格算法數值求解DEM的Poisson方程，得到DEM反演。
文檔編號G01S13/90GK1415975SQ0213773
公開日2003年5月7日 申請日期2002年10月31日 優先權日2002年10月31日
發明者金亞秋, 羅霖 申請人:復旦大學