专利名称：一种地形测量方法
技术领域：
本发明涉及信息获取与处理技术领域，特别是一种利用单路径全极化合成孔径雷达复图像数据的地形测量方法。
背景技术：
利用合成孔径雷达的全极化信息进行地形测量是一种新的技术，如图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日
发明者洋 李, 文 洪, 芳 曹, 王彦平, 吴一戎 申请人:中国科学院电子学研究所