专利名称：单次飞行全极化合成孔径雷达图像反演地面数字高程的方法
技术领域：
本发明属于空间遥感与通信技术领域，具体涉及一种星载合成孔径雷达反演地面数字高程的方法。
但是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日
发明者金亚秋, 罗霖 申请人:复旦大学