单次飞行全极化合成孔径雷达图像反演地面数字高程的方法

摘要

本发明是单次飞行全极化合成孔径雷达图像反演地面数字高程的方法。利用全极化散射Mueller矩阵解，推导取向角ψ的迁移，借助Euler角的变换，建立了由3个Stokes参数确定的ψ迁移与倾斜表面像素的射程角β和水平方位角γ之间的联系。在只有单次飞行SAR数据的条件下，用自适应阈值的二值化和图像形态学细化算法来确定γ以及β，以及地表面射程方位坡度和水平方位坡度。由此用完整多重网格算法数值求解DEM的Poisson方程，得到DEM反演。

说明书

发明领域

本发明属于空间遥感与通信技术领域，具体涉及一种星载合成孔径雷达反演地面数字高程的方法。

背景技术

星载合成孔径雷达(SAR)极化测量与成像技术是90年代空间对地观测与遥感的最重要的进展之一，它能全天候高分辨率观测地面与海洋的散射回波，并能成像。近年来在SAR的基础上又发展了相干飞行或双天线飞行的干涉SAR(INSAR)技术。美国海军实验室Lee等人在2000年用近似的散射振幅函数证明了地面起伏产生回波的极值处于入射极化角ψ≠0。这种ψ≠0的迁移可用来推算地面起伏高程。地面起伏高程是对地观测与地形测量的重要内容。Lee等人用两次相干飞行的INSAR技术实现了数字地面高程(DEM)的反演。

但是INSAR各相关技术难度高，并且不是每一处地形都能容易得到符合要求的相干的SAR图像。

发明内容

本发明的目的在于提出一种用单次飞行的极化SAR数据图像来实现DEM反演的方法。

用全极化SAR图像对地表分类已经得到了广泛的研究。从全极化SAR图像数据，可以得到2×2维的复散射振幅函数F_pq(p，q＝v，h)和4×4维的实Mueller矩阵M_ij(i，j＝1，...，4)([1]Jin Y.Q.，Electromagnetic scattering modelling forquantitative remote sensing，1994，(Singapore：World 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，23：322-333)。用实同极化和零交叉极化的散射振幅函数的假设条件推导了用散射振幅函数表示的ψ的迁移。由于要同时确定水平方位角和射程角，所以需要两次或多次相近飞行的SAR图像数据来确定表面坡度。近年来，用两次相近飞行的SAR和INSAR数据来反演地面数字高程(DEM)的方法已有很好的实例(见文献[3])。

本发明提出的用单次飞行的极化SAR数据图像实现DEM反演的方法，是用全极化散射Mueller矩阵解，推导ψ的迁移，将其表示为三个散射Stokes参数I_vs，I_hs，U_s的函数，而I_vs，I_hs，U_s是由SAR测量得到的；利用主坐标系和像素元的局部坐标系之间的欧拉(Euler)角变换，把入射波极化取向角ψ与射程角β、水平方位角γ以及SAR观测的几何构造直接联系起来；利用倾斜地面的水平方位排列在SAR图像中表现出来的线性纹理，用图像形态学细化算法确定图像中各像素元的地形水平方位角γ，再用I_vs，I_hs，U_s计算ψ迁移来确定射程角β；最后用β和γ来确定各像素元地形水平方位坡度与射程方位坡度；最后用完整多重网格算法数值求解地面高程阿松(Poisson)方程，完成SAR观测区域地表DEM的反演。

本发明的内容进一步描述如下：

1、用散射Stokes参数表达取向角ψ迁移

考虑一极化电磁波以角度(π-θ_i，_i)入射到地面上，其极化可用椭圆角χ和取向角ψ来表征(见文献[1])。散射场可写为其中2×2维的复散射振幅函数F_pq(p，q＝v，h)可由全极化散射测量技术得到，它包含了散射层的体散射和面散射。由式(1)，可推算得到Stokes矢量的实Mueller矩阵解(见文献[1])，其中归一化的入射Stokes矢量可写成

 I₁(χ，ψ)＝[I_vi，I_hi，U_i，V_i]^T

 ＝[0.5(1-cos 2χcos 2ψ)，0.5(1+cos 2χsin 2ψ)，-cos 2χsin 2ψ，sin 2χ]^T      (3)Mueller矩阵元素M_ij(i，j＝1，Λ，4)由(p，q，s，t＝v，h)的实部或虚部的系综平均表示，表达式可见文献[1]。应注意到，同极化项总是比交叉极化项(s≠t)要大得多。

同极化的后向散射系数σ_c(θ_i，π+_i；π-θ_i，_i)可表示为(见文献[1])，

                     σ_c＝4πcos_iP_n                               (4)其中P_n＝0.5[I_vs(1-cos2χcos2ψ)+I_hs(1+cos2χcos2ψ)+U_scos2χsin2ψ+V_ssin2χ]

                                                                  (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，35：813-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)。

在P_n/ψ＝0处，σ_c取极大值。假设χ＝0的对称情形，由(5)式可得到

            0＝-(I_hs-I_vs)sin2ψ+U_scosψ

               +0.5(I_hs+I_vs)′+0.5U′_ssin 2ψ+0.5(I_hs-I_vs)′cos2ψ            (6)其中上标′表示/ψ。由(2)式的比较(6)式右端各项，有 _hs-I_vs)～0.5(M₁₁+M₂₂)cos 2ψ～0.5(<|F_vv|²>+<|F_hh|²>)cos 2ψ                     (7c)

可看到，(7a)与(7b)均远小于(7c)，而(8a)亦远小于(8b)。因此(6)式右端的最后三项可略去。于是，本发明中得到σ_c的极大值的ψ迁移取为： $>>tan>2>\psi >=>>>U>s>>>>I>hs>>->>I>vs>>>>->->->>(>9>)>>>s>由此可以看出，第3个Stokes参数U$_s≠0引起了ψ的迁移。

顺便说，如果U_s和I_hs-I_vs都趋向于0，比如，当散射来自于均匀取向散射体或各向同性散射介质，如厚的植被顶层时，式(9)就不能很好地确定ψ的迁移。

进一步地，从式(2)可得到

        I_hs-I_vs～0.5(M₂₂-M₁₁)＝0.5(<|F_hh|²>-<|F_vv|²>)              (10b)2、用Euler角变换得到射程角和水平方位角入射波的极化矢量在主坐标系中定义为(见文献[1])， $>>ver>>h>^>>i>>=>>ver>>z>^>>\times >ver>>k>^>>i>>>>|ver>>z>^>>\times >ver>>k>^>>i>>|>>>=>->sin>>\phi >i>ver>>x>^>>+>cos>>\phi >i>ver>>y>^>>,>ver>>v>^>>i>>=>ver>>h>^>>i>>\times >ver>>k>^>>i>>->->->>(>12>)>>>s>其中入射波矢量为$ >>ver>>k>^>>i>>=>sin>>\theta >i>>cos>>\phi >i>ver>>x>^>>+>sin>>\theta >i>>sin>>\phi >i>ver>>y>^>>->cos>>\theta >i>ver>>z>^>>->->->>(>13>)>>>s>当该像素地表面有如图1所示的局部坡度时，由局部法向量定义的局部坐标下的极化矢量可表示为$ >>ver>>h>^>>b>>=>>>ver>>z>^>>b>>\times >ver>>k>^>>ib>>>>|>ver>>z>^>>b>>\times >ver>>k>^>>ib>>|>>>,>ver>>v>^>>b>>=>ver>>h>^>>b>>\times >ver>>k>^>>ib>>->->->>(>14>)>>>s>通过两个坐标系和之间的Euler角(\beta ，\gamma )变换(见文献(1))，有$ >ver>>x>^>>=>cos>>\gamma >cos>\beta >ver>>x>^>>b>>+>sin>\gamma >ver>>y>^>>b>>->cos>\gamma \beta >ver>>z>^>>b>>->->->>(>15>a>)>>>s>$ >ver>>y>^>>=>->sin>>\gamma >cos>\beta >ver>>x>^>>b>>+>cos>\gamma >ver>>y>^>>b>>+>sin>>\gamma >sin>\beta >ver>>z>^>>b>>->->->>(>15>b>)>>>s>$ >ver>>z>^>>=>sin>\beta >ver>>x>^>>b>>+>cos>\beta >ver>>z>^>>b>>->->->>(>15>c>)>>>s>将以上三式代入式(12)，(13)中，可得$ >>ver>>h>^>>i>>=>>(>->cos>>\gamma >cos>>\beta >sin>>\phi >i>>->sin>>\gamma >cos>>\beta >cos>>\phi >i>>)>>ver>>x>^>>b>>+>>(>->sin>>\gamma >sin>>\phi >i>>+>cos>>\gamma >cos>>\phi >i>>)>>ver>>y>^>>b>>->->->>(>16>)>>>s>$ >>+>>(>cos>>\gamma >sin>>\beta >sin>>\phi >i>>+>sin>>\gamma >sin>>\beta >cos>>\phi >i>>)>>ver>>z>^>>b>>>s>$ >>ver>>k>^>>ib>>=>>(>cos>>\gamma >cos>>\beta >sin>>>\theta >i>>cos>>\phi >i>>->sin>>\gamma >cos>>\beta >sin>>>\theta >i>>sin>>\phi >i>>->sin>>\beta >cos>>\theta >i>>)>>ver>>x>^>>b>>>s>$ >>+>>(>sin>>\gamma >sin>>\theta >i>>cos>>\phi >i>>+>cos>>\gamma >sin>>>\theta >i>>sin>>\phi >i>>)>>ver>>y>^>>>b>->->->>(>17>)>>>>>s>$ >>+>>(>->cos>>\gamma >sin>>\beta >sin>>>\theta >i>>cos>>\phi >i>>+>sin>>\gamma >sin>>\beta >sin>>>\theta >i>>sin>>\phi >i>>->cos>>\beta >cos>>\theta >i>>)>>ver>>z>^>>b>>>s>将式(17)代入式(14)中，该像素地表面的极化矢量可表示为$ >>ver>>h>^>>b>>=>>>a>ver>>x>^>>b>>+>b>ver>>y>^>>b>>>>>a>2>>+>>b>2>\; >>->->->>(>18>)>>>s>其中a＝-(sin\gamma sin\theta $$$$$$$$$$$$_icosφ_i+cosγsinθ_isinφ_i)

                                                       (19)b＝(cosγcosβsinθ_icosφ_i-sinγcosβsinθ_isinφ_i-sinβcosθ_i)于是，由式(16)，(18)可以得到入射波极化的取向角ψ为 $>>cos>\psi >=>ver>>h>^>>b>>·>ver>>h>^>>i>>=>>>cos>>\beta >sin>>\theta >i>>->cos>>(>\gamma >+>>\phi >i>>)>>sin>\beta >cos>>\theta >i>>>>>a>2>>+>>b>2>\; >>>(>20>a>)>>>s>$ >>tan>\psi >=>>>tan>>\beta >sin>>(>\gamma >>>+>\phi >>i>>)>>>>sin>>\theta >i>>->cos>>\theta >i>>tan>>\beta >cos>>(>\gamma >+>>\phi >i>>)>>>>->->->>(>20>b>)>>>s>通常令入射\phi $$_i＝0，是射程方向，是水平方位方向，注意在文献[2]中，水平方位角定义为和之间的夹角。于是，该像素地表面的射程方位坡度S_R和水平方位坡度S_A可定义为(见文献[2])：S_R＝tanβcosγ

                        (21)S_A＝tanβsinγ因为在式(20b)中，用一个ψ不能同时确定β和γ两个角度，所以通常需要两次或多次相近飞行的SAR图像数据，用式(20b)建立两个方程来求解β和γ，这在文献[2]中已有讨论。

3、确定SAR图像中各像素的水平方位角

如果只有单次飞行的SAR图像数据，除了用ψ与I_vs，I_hs，U_s的关系式(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>\beta >=>>>tan>>\psi >sin>>\theta >i>>>>sin>\gamma >+>tan>>\psi >cos>>>\theta >i>>cos>\gamma >>>>(>22>)>>>s>其中\theta$_i由SAR观测的几何结构给出。由β和γ得到式(21)中的S_A和S_R，在此基础上采用文献[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，5：416-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方程为  

                         ²φ(x，y)＝ρ(x，y)                           (23)其中²是Laplace算子²/x²+²/y²，φ(x，y)表示在(x，y)处的高度值，并令初始值为零。方程右端的源函数ρ(x，y)是由输入坡度函数S(x，y)得到的表面曲率值。我们假定x方向为射程方位，而y方向为水平方位，则dS^R(x)/dx和dS^A(y)/dy分别包含了射程方位和水平方位上的表面曲率(上标A和R分别表示水平方位和射程方位)。

水平方位和射程方位的源函数的离散值由下式给出， $>sup>>\rho >ij>Rsup>>=sup>>S>ij>Rsup>>-sup>>S>>i>->1>,>j>>Rsup>>,sup>>\rho >ij>Asup>>=sup>>S>ij>Asup>>-sup>>S>>i>,>j>->1>>Asup>>>(>0>\le >i>\le >M>,>0>\le >j>\le >N>)>>->->->>(>24>)>>>s>$

假定输入坡度的误差满足高斯(Gauss)分布，地面的高程可看成坡度的积分，那么地面高程的误差也满足Gauss分布。因此，一个服从χ²分布的统计量可看成期望解与实际值相符合的量度标准(见文献[5])。这也是干涉SAR中相位展开算法的标准假设。χ²的值由下式给出 $>>>x>2>>=>>\Sigma >>i>,>j>\; >>>(>>\phi >ij>>->>\phi >>i>->1>,>j>>>-sup>>S>ij>Rsup>>)>>2>>+>>\Sigma >>i>,>j>\; >>>(>>\phi >ij>>->>\phi >>i>,>j>->1>>>-sup>>S>ij>Asup>>)>>2>>->->->>(>25>)>>>s>$

求解φ_ij的偏微商值与式(24)所定义的偏微商值的差必须为最小。对χ²求导取零得到(φ_i+1，j-2φ_ij+φ_i-1，j)+(φ_i，j+1-2φ_ij+φ_i，j-1)＝ρ_ij              (26)式(26)就是需要求解的离散Poisson方程。

离散的源函数ρ_ij由下式给出 $>>>\rho >ij>>=>>(sup>>S>>i>+>1>,>j>>Rsup>>-sup>>S>ij>Rsup>>)>>+>>(sup>>S>>i>,>j>+>1>>Asup>>-sup>>S>ij>Asup>>)>>->->->>(>27>)>>>s>$

文献[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，34：728-738；[12]J.Demmel，Applied NumericalLinear Algebra，1997，Society for Industrial and Applied Mathematics)对于M×N的矩形网格，其计算量为(MN)ln(MN)。FMG算法的具体步骤可参见[11]。计算框图见图13所示。

附图说明

图1为模型的几何结构，其中图1(a)为地形的几何形态图状，图1(b)为地形的等高线图。

图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为本发明计算程序框图。

具体实施方式

实施例1，广东惠州地区。

采用美国加州理工学院喷气推进实验室(JPL)提供的1994年10月4日SAR SIR-C对中国广东省观测的L波段数据，图像中心位置约在(114.561°E，22.815°N)，每个像素的大小是12.5×12.5m²。SAR在图像中心的入射角是32.85°。图3是从原图中截取的大小为512×512像素的后向散射vv，hh，hv总功率图像，其中心位置为(114.7826°E，23.1359°N)，位于中国广东省的惠州地区。

从图3可以看出南北走向的丘陵起伏。为了突出地形的主要走向，我们略去了一些更为细微的局部起伏变化。注意，这里定义的水平方位角γ是坡度走向与方位的夹角(如图1(a)所示)。

S_A和S_R的计算结果如图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.5m²。SAR在图像中心的入射角是32.85°。图8是从原图中截取的大小为512×512像素的后向散射vv，hh，hv总功率图像，其中心位置为(120.8626°E，23.6659°N)，位于中国台湾的阿里山地区。

反演的DEM与原SAR图像所反映的直觉的地形结构以及十分简略的地面高程图^[13]是基本一致的，其主要地形的走向与SAR图像的明暗纹理也是吻合的。