专利名称：极化干涉合成孔径雷达数据估计地形的方法及其软件的制作方法
技术领域：
本发明属于合成孔径雷达遥感信息处理技术领域，具体涉及一种极化干涉合成孔径雷达数据估计地形的方法及其软件，它利用极化分类技术优势将土地类型构成复杂区域的数据划分成若干子块处理，降低地形干涉相位估计的难度的同时，提高地形估计效率；在估计各个子块的地形干涉相位中，对单散射中心区域和多散射中心区域分别采用不同极化干涉方法估计地形干涉相位；不同子块之间由于多散射中心区域类型的差异选择不同方法提取地形干涉相位，更适合土地类型复杂区域的处理需求；从而保证了地形干涉相位、地形的估计精度。各子块地形干涉相位并行处理，缩短了地形估计时间、提高了地形估计效率。本发明方法的配套软件包含完成本发明方法中各处理步骤的代码段。
背景技术：
作为一种成熟的测量技术，干涉合成孔径雷达技术(InSAR)利用某单一种极化状态的干涉相位估计地形。除了系统参数、轨道测量准确性外，InSAR数据估计地形的性能主要取决于干涉相位质量。作为干涉相位的质量指标，干涉相干性越大表示干涉相位质量越高。对于单散射中心区域，InSAR数据的干涉相干性较高可获得准确的地形估计；但对于多散射中心区域如森林、冰雪区、城区等而言，由于除地形外其它散射中心的影响，InSAR数据的干涉相干性明显下降，造成InSAR数据的地形估计中出现明显误差(Cloude S.R.and Papathanassiou K.P.(1998)，Polarimetric SARInterferometry)。因而，InSAR数据估计多散射中心区域地形性能不佳。
极化干涉合成孔径雷达技术利用全极化数据进行干涉处理。它利用干涉相位对极化状态的依赖性，不但最大限度地保持干涉相位质量(CloudeS.R.and Papathanassiou K.P.(1998)，Polarimetric SAR Interferometry)、保证地形估计精度；还能去除地形外其它散射中心的影响，用于估计森林、冰雪、城区等多散射中心区域的地形干涉相位(Brandfass M.(2002)，Generation of bald earth digital elevation models as applied to polarimetricSAR Interferomety)，从而获得准确的地形信息。因而，极化干涉数据估计地形的方法比干涉数据估计地形的方法优势更明显，其在土地类型复杂区域的地形估计中具有巨大的优势和应用潜力。
极化干涉数据估计地形干涉相位的方法主要可划分为三个大类相干最优技术、基于相干散射模型的方法以及基于估计来波指向的技术。
相干最优技术通过改变极化状态、补偿时间去相关的影响优化干涉相干性，分辨出具有最高干涉相干性的干涉相位，并将其作为地形干涉相位(Mounira Souissi Ouarzeddine，(2002).Generation of Digital Models usingPolarimetric SAR Interferometry)。该技术适用于单散射中心区域或者地形散射中心占主导地位的土地类型。
基于相干散射模型的方法是利用复相干系数与地形干涉相位等参数间的联系，根据不同极化下相干系数的拟合直线估计地形干涉相位(Papathanassiou K.P.and Cloude S.R.(2001)，Single-Baseline PolarimetricSAR Interferometry)。该方法基于随机指向大量散射中心覆盖地面的二层简化假设，适用于森林、冰雪区等多散射中心区域地形干涉相位估计。
基于估计来波指向的技术假设存在若干相互独立的、稳定的散射中心(Guillaso S.，et al.(2003)，Urban area analysis based on ESPRIT/MUSICmethods using polarimetric interferometric SAR)，引入谱函数用于估计各个散射中心的干涉相位，地形干涉相位是对应于最大谱函数的干涉相位(Kasilingam D.，Nomula M.，Cloude S.(2002)，A technique for removingvegetation bias from polarimetric SAR interferometry)。鉴于基于估计来波指向的技术的假设，该技术适于存在稳定的、相互独立的多个散射中心的城区，估计城区的地形干涉相位。
目前，极化干涉数据估计地形的方法未充分考虑到土地类型的多样性以及不同的地形干涉相位估计技术适用范围，对所有数据统一采用某一方法估计地形干涉相位；未将单散射中心区域与多散射中心区域分别处理、未考虑在不同类型的多散射中心区域地形干涉相位估计差异，未充分发挥极化干涉数据在估计地形干涉相位的优势，影响复杂土地类型区域的地形估计准确性。
目前包含极化干涉数据处理功能的软件主要有柏林科技大学开发的RATS软件、瑞士的GAMMA遥感与咨询研究组开发的GAMMA软件、德国的CREASO公司的SARscape软件、荷兰代夫特科技大学地球观测与空间系统研究所开发的Doris软件、欧空局资助法国雷恩大学发开的POLSARPRO软件、中国林业科学研究院资源信息研究所等单位合作开发的SARINFORS软件。上述软件均包含极化干涉数据处理模块，但是它们统一采用相干最优或基于相干散射模型方法估计地形干涉相位，而未具备在不同土地类型区域选择不同极化干涉方法提取地形干涉相位的功能。
发明内容
本发明的目的在于提出一种极化干涉合成孔径雷达数据估计地形的方法及软件，以克服现有技术存在的缺陷，既降低土地类型构成复杂的区域地形估计的难度，又提高该类区域地形估计效率。。
为达到上述目的，本发明的技术解决方案是 一种基于极化干涉合成孔径雷达数据的地形估计方法，其包括步骤 A.对极化干涉数据滤波；B.利用极化分类技术区分不同的土地类型；C.将数据划分为土地类型构成简单的若干子块；D.分别估计各子块的地形干涉相位；E.并行处理各子块地形干涉相位，估计各子块地形；F.组合各子块地形，计算完整的地形结果。
所述的方法，其所述A步，对极化干涉数据滤波，滤波后的数据无需补偿平地相位，极化干涉数据滤波抑制了噪声对地形干涉相位估计的不利影响，其处理包括以下步骤 步骤S11基于成像几何关系，从主、副天线之间的平均水平基线距Bh、平均垂直基线距Bv、主天线相对参考平面高度h，近距rn、斜距采样间隔rs，雷达信号波长λ，估计方位向标号为m、斜距向标号为n的位置的平地相位φj(m，n) 其中π表示180°角对应的弧度值，符号·表示乘运算，m、n是方位向、斜距向位置标号，它们是自然数且满足1≤m≤M，1≤n≤N；符号M、N表示极化干涉数据在方位向、斜距向总的像素个数；估计平地相位是为抑制相邻像素之间平地相位差异对极化干涉数据滤波造成的不利影响； 步骤S12去除平地相位对极化干涉数据的影响 从而增强相邻像素之间极化干涉数据的连续性；其中注脚数字2表示副天线数据；p、q表示接收、发送信号的极化方式，是水平极化方式h或者垂直极化方式v；S2，pq(m，n)表示方位向标号为m、斜距向标号为n的位置的pq极化副天线数据；S′2，pq(m，n)表示去除平地相位影响后的方位向标号为m、斜距向标号为n的位置的pq极化副天线数据；e表示自然对数，符号j表示虚数单位，符号
表示自然对数e的-j·φj(m，n)次幂； 步骤S13在待滤波位置为中心的滤波窗口内，对极化干涉数据进行平滑处理，计算方位向标号为m、斜距向标号为n的位置滤波后的极化干涉数据T11(m，n)、T22(m，n)和T12(m，n) 其中上标H表示矩阵的共轭转置；括号符号<>表示平均。
表示由主天线数据构成的矢量 其中上标T表示矩阵转置，注脚数字1表示主天线数据，S1，hh(m，n)表示方位向标号为m、斜距向标号为n的位置的hh极化主天线数据；S1，hv(m，n)表示方位向标号为m、斜距向标号为n的位置的hv极化主天线数据；S1，vh(m，n)表示方位向标号为m、斜距向标号为n的位置的vh极化主天线数据；S1，vv(m，n)表示方位向标号为m、斜距向标号为n的位置的vv极化主天线数据； 符号
是由去除平地相位后的副天线数据S′2，hh(m，n)、S′2，hv(m，n)、S′2，vh(m，n)、S′2，vv(m，n)构成，其表达式为 滤波后的T12(m，n)无需补偿平地相位，直接用于后续处理。
所述的方法，其所述B步，将极化分类技术应用于主、副天线数据平均结果来区分不同的土地类型，利用极化分类是为分析数据覆盖的土地类型复杂程度，为把数据划分为若干土地类型构成简单的子块提供依据，其具体处理步骤如下 步骤S21计算主、副天线全极化数据T11(m，n)与T22(m，n)的算术平均 其中T(m，n)表示主、副天线全极化数据平均结果； 步骤S22从T(m，n)计算单次散射、二次散射以及体散射强度S(m，n)、D(m，n)、V(m，n) S(m，n)＝T(m，n)(1，1) D(n，n)＝T(m，n)(2，2)(7) V(m，n)＝T(m，n)(3，3) 其中注脚(1，1)表示矩阵T(m，n)中第1行、第1列处数据，注脚(2，2)、(3，3)分别表示矩阵T(m，n)第2行、第2列处数据和第3行、第3列处数据； 步骤S23依次将单次散射、二次散射以及体散射强度的分贝值LS(m，n)、LD(m，n)、LV(m，n) LS(m，n)＝10·log10(S(m，n)) LD(m，n)＝10·log10(D(m，n)) (8) LV(m，n)＝10·log10(V(m，n)) 均匀量化为蓝色、红色、绿色通道整数颜色值B(m，n)、R(m，n)、G(m，n)
其中log10()表示计算10为底的数据的对数值，符号
计算小于等于数据的最大整数，符号min表示指定的对应于颜色值0的数据，max表示指定的对应于颜色值255的数据，根据蓝色、红色、绿色通道的颜色值B(m，n)、R(m，n)、G(m，n)生成极化分类结果图。
所述的方法，其所述C步，将数据划分为土地类型构成简单的若干子块，该步骤降低了土地类型构成的复杂度、降低了地形干涉相位估计难度，并兼顾不同土地类型区域采用不同方法估计地形干涉相位需求；根据极化分类结果图，分析土地类型构成的复杂度，结合地面控制点的分布将其划分为土地类型构成简单且至少包含一个地面控制点的若干子块。
所述的方法，其所述D步，分别估计各子块地形干涉相位，且子块内单散射中心区域、多散射中心区域采用不同的极化干涉方法提取地形干涉相位，各子块数据计算地形干涉相位的步骤为 步骤S41标记子块内多散射中心区域和单散射中心区域，该步骤是在多散射中心区域、单散射中心区域应用不同方法提取地形干涉相位，从而满足不同土地类型地形干涉相位估计需求；多散射中心区域判定原则为 V(m，n)≥AV (10) 其中AV表示判定阈值，非多散射中心区域标记为单散射中心区域，阈值AV由对应于极化分类结果图中某一多散射中心区域的体散射强度最小值确定； 步骤S42在单散射中心区域，应用极化干涉相干最优方法通过求解具有最大干涉相干系数的干涉相位估计地形干涉相位； 步骤S43利用指定极化干涉方法估计子块内多散射中心区域的地形干涉相位，当子块内多散射中心区域为森林或者冰雪区，指定相干散射模型方法计算地形干涉相位；当子块内多散射中心区域为城区，指定CAPON算法计算地形干涉相位。
所述的方法，其所述步骤S42中的极化干涉相干最优方法，是利用T(m，n)提取散射中心信息，将具有最大干涉相干系数的干涉相为作为地形干涉相位，具体的计算步骤为 步骤S61基于特征分解，将矩阵T(m，n)分解为三个相互独立的部分 其中vi、
分别表示矩阵T(m，n)第i个特征值和对应的特征矢量(i＝1，2，3)； 步骤S62计算对应于上述三个相互独立部分的复干涉相干系数 步骤S63地形干涉相位φg(m，n)是对应于最大干涉相干系数的干涉相位 其中Arg()表示计算复数的幅角主值，||表示计算复数的模，
表示三组复干涉相干系数
之中具有最大模值的复干涉相干系数。
所述的方法，其所述步骤S43中的相干散射模型方法，是利用极化干涉相干最优方法获取的三组复干涉相干系数，拟合复干涉相干系数所在直线，地形干涉相位由该直线与单位圆交点确定，其具体的估计步骤为 步骤S71利用极化干涉相干最优方法估计三组复干涉相干系数
(i＝1，2，3)； 步骤S72采用最小二乘方法，用三组复干涉相干系数拟合复平面内直线方程y＝a·x+b，具体的计算方法为 其中x、y表示复数据的实部、虚部，a、b表示直线方程中的一次项、常数项系数，符号Re()表示复数的实部，Im()表示复数的虚部； 步骤S73计算直线y＝a·x+b与单位圆x2+y2＝1的两个交点(x1，y1)、(x2，y2)，其中x1、y1表示第一个交点的实部、虚部，x2、y2表示第二个交点的实部、虚部。将这两个交点表示为复数
其中φ1、φ2表示与第一个交点、第二各交点对应的复数的幅角主值； 步骤S74计算hh+vv极化、hv+vh极化下干涉复相干系数
(15) 步骤S75根据hv+vh极化状态的干涉相干系数距地形干涉相位最远的假设，确定地形干涉相位φg(m，n) (16) 其中Sign()表示计算数据的符号，当数据大于等于零时，Sign()的取值为1；当数据小于零，Sign()取值为-1。
所述的方法，其所述步骤S43中的CAPON算法，是利用分块矩阵以及特征分解技术，通过计算谱函数峰值求解地形干涉相位；地形干涉相位估计经过粗略峰搜索、精细谱峰搜索两阶段计算获得，具体的计算步骤为 步骤S81先利用T11(m，n)、T22(m，n)、T12(m，n)组合为6×6矩阵R6 步骤S82计算R6的逆矩阵
并将
划分为四个3×3矩阵
和
其中上标-1表示矩阵求逆； 步骤S83粗略搜索达到谱函数峰值处相位，粗略搜索是在(-π，π]均匀选取K个角度采样点 其中符号k为采样点标号，它是小于等于K的自然数；计算各个角度采样点φk谱函数P(φk) 其中Γmin()表示计算矩阵非负实数特征值中最小的特征值，比较搜索具有最大谱函数的角度采样点
，其中kmax表示对应于最大谱函数的角度采样点标号； 步骤S84精细搜索达到谱函数峰值处相位，计算地形干涉相位；在粗略搜索到的对应于最大谱函数的角度采样点
及与
最近临近前后两采样点
范围内，再进行采样数为K的均匀精细采样，各个精细采样点φ′k 计算各个精细采样点的谱函数值P(φ′k)，比较搜索具有最大谱函数的精细采样点；地形干涉相位φg(m，n)是使得谱函数P(φ′k)达到最大的精细角度采样点相位值。
所述的方法，其所述E步，并行处理各子块地形干涉相位，估计各子块的地形，各子块并行处理缩短了地形估计时间、提高地形估计效率；估计各子块的地形，按照常规干涉相位估计地形的方法，分别经过相位滤波、相位展开、平地相位补偿、地形高度估计。
所述的方法，其所述F步，组合子块地形，计算完整的地形结果，按照各子块在原始数P(φ′k)据中的位置，排列各子块的地形结果，从而组合成完整的地形结果。
一种用于所述的方法的极化干涉合成孔径雷达数据估计地形的软件，该软件是由记录在计算机刻度介质上的程序代码构成，该程序代码用于实现极化干涉合成孔径雷达数据估计地形的方法，其特征在于该程序代码包括 -用于完成步骤S11平地相位估计的代码段； -用于完成步骤S12去除平地相位对极化干涉数据影响的代码段； -用于完成步骤S13极化干涉数据平滑的代码段； -用于完成极化分类的代码段； -用于完成子块划分代码段； -用于完成地形干涉相位估计代码段； -用于完成相位滤波、相位展开、平地相位补偿、地形高度估计的代码段； -用于完成组合子块地形的代码段； 该极化干涉合成孔径雷达数据估计地形的软件可在中国科学院电子学研究所微波成像技术国家级重点实验室购买。
所述的软件，其所述平地相位估计代码段，是在A步中的步骤S11的处理过程，计算各位置(m，n)的平地相位φf(m，n)。
所述的软件，其所述去除平地相位对极化干涉数据影响的代码段，是在A步中步骤S12的处理过程，去除平地相位对指定极化状态pq的各位置(m，n)的副天线数据的影响其中pq为hh、hv、vh或者vv。
所述的软件，其所述极化干涉数据平滑的代码段，是在A步中步骤S13的处理过程，该代码段在各位置(m，n)为中心的指定大小Mf×Nf的滤波窗口内计算T11(m，n)、T22(m，n)和T12(m，n)
其中m′、n′是指定范围内的整数。
所述的软件，其所述极化分类的代码段，是在B步中步骤S21～步骤S23的处理过程，其中min、max是输入的对应于颜色值0、255的散射强度分贝值。
所述的软件，其所述子块划分代码段，是在C步的处理过程通过设置各子块的左上角的位置(m1，n1)以及各子块大小Ml×Nl，计算各子块内数据T11(m，n)、T22(m，n)和T12(m，n)。
所述的软件，其所述地形干涉相位估计代码段，是在D步的步骤S41～步骤S43处理过程其中判定阈值AV以及处理估计多散射中心区域的地形干涉相位的方法需预先设定。
所述的软件，其所述相位滤波、相位展开、平地相位补偿、地形高度估计的代码段，是在步骤E中相位滤波、相位展开、平地相位补偿、地形高度估计处理过程。
所述的软件，其所述组合子块地形的代码段，是在F步的处理过程按照各子块左上角位置(ml，nl)以及各子块大小Ml×Nl将各子块的地形高度结果排列为完整的地形高度结果。
本发明提出极化干涉合成孔径雷达数据估计地形的方法，其创新之处在于充分考虑到土地类型的多样性与地形干涉相位估计需求，利用极化分类技术区分不同的土地类型，并将其划分为土地类型构成简单的若干子块分别处理，从而降低地形估计难度、提高地形估计效率；在估计地形干涉相位过程中，根据对不同土地类型区域的特点选用不同的方法估计地形干涉相位，满足不同土地类型区域数据处理需求，保证了地形干涉相位估计、地形估计准确性。各子块地形干涉相位并行处理，缩短了地形估计时间、提高了地形估计效率。
本发明方法的配套软件是专门为本发明方法设计的，它包含完成本发明方法中各处理步骤的代码段。
图1本发明方法所述的极化干涉数据示意图。
图2本发明提出的极化干涉合成孔径雷达数据的地形估计方法的流程图。
图3本发明软件显示实例中德国ESAR系统L波段数据极化分类结果图。
图4本发明实例德国ESAR系统L波段数据划分为子块F1、F2、F3、F4、F5、F6的示意图。
图5a本发明实例德国ESAR系统L波段数据划分为子块F1应用本发明方法估计的地形干涉相位结果。
图5b本发明实例德国ESAR系统L波段数据划分为子块F2应用本发明方法估计的地形干涉相位结果。
图5c本发明实例德国ESAR系统L波段数据划分为子块F3应用本发明方法估计的地形干涉相位结果。
图5d本发明实例德国ESAR系统L波段数据划分为子块F4应用本发明方法估计的地形干涉相位结果。
图5e本发明实例德国ESAR系统L波段数据划分为子块F5应用本发明方法估计的地形干涉相位结果。
图5f本发明实例德国ESAR系统L波段数据划分为子块F6应用本发明方法估计的地形干涉相位结果。
图6a本发明实例德国ESAR系统L波段数据划分为子块F1应用本发明方法估计的地形结果。
图6b本发明实例德国ESAR系统L波段数据划分为子块F2应用本发明方法估计的地形结果。
图6c本发明实例德国ESAR系统L波段数据划分为子块F3应用本发明方法估计的地形结果。
图6d本发明实例德国ESAR系统L波段数据划分为子块F4应用本发明方法估计的地形结果。
图6e本发明实例德国ESAR系统L波段数据划分为子块F5应用本发明方法估计的地形结果。
图6f本发明实例德国ESAR系统L波段数据划分为子块F6应用本发明方法估计的地形结果。
图7a本发明实例中德国ESAR系统数据应用本发明方法计算出的完整地形高度结果。
图7b本发明实例中应用本发明方法与HH极化的干涉数据计算的地形高度结果在方位位置为910的剖面处对比图。
具体实施例方式 首先介绍一下极化干涉数据，它是由两组全极化数据构成的如图1所示的数据。全极化数据是指全极化天线系统发送水平极化(h)、垂直极化(v)发送雷达信号，并采用h、v极化接收信号获得的hh、hv、vh、vv四种极化状态下的原始回波数据，再经过成像处理获得的相同方位向、斜距向位置的hh、hv、vh、vv极化状态复数据Shh(m，n)、Shv(m，n)、Svh(m，n)、Svv(m，n)构成。极化干涉合成孔径雷达数据是由不同时间采集的或者同时在不同位置采集的两组全极化合成孔径雷达数据构成，且每组的各个极化状态下的雷达信号要分别经过成像处理后，还需要经过极化定标处理、图像配准和重采样处理，使得方位向、斜距向位置为m、n的数据对应同一目标，从而构成极化干涉数据S1，hh(m，n)、S1，hv(m，n)、S1，vh(m，n)、S1，vv(m，n)、S2，hh(m，n)、S2，hv(m，n)、S1，vh(m，n)、S1，vv(m，n)描述。其中，标号1表示图像配准中作为参考的一组全极化数据，本发明中叫做主天线数据，采集主天线数据的全极化天线系统叫做主天线；标号2表示未作为参考的一组全极化数据，本发明中叫做副天线数据，采集副天线数据的全极化天线系统叫做副天线。
如图2所示，本发明提出的极化干涉合成孔径雷达数据的估计地形高度方法主要包含六个步骤 步骤A极化干涉数据滤波，其包含以下步骤 步骤S11基于成像几何关系，从主、副天线之间的平均水平基线距Bh、平均垂直基线距Bv、主天线相对参考平面高度h，近距点斜距rn、斜距采样点间隔rs，雷达信号波长λ，估计方位向标号为m、斜距向标号为n的位置的平地相位φf(m，n) 其中m、n是方位向、斜距向位置标号，它们满足1≤m≤M，1≤n≤N。符号M、N表示极化干涉数据在方位向、斜距向总的像素个数。估计平地相位是为抑制相邻像素之间平地相位差异对极化干涉数据滤波的造成的不利影响； 步骤S12去除平地相位对极化干涉数据的影响 从而增强相邻像素之间极化干涉数据的连续性。其中注脚数字2表示副天线数据；p、q表示接收、发送信号的极化方式，它们可表示水平极化方式h，或者垂直极化方式v；S2，pq(m，n)表示方位向标号为m、斜距向标号为n的位置的pq极化副天线数据；S′2，pq(m，n)表示去除平地相位影响后的方位向标号为m、斜距向标号为n的位置的pq极化副天线数据；e表示自然对数，符号j表示虚数单位，符号
表示自然对数e的-j·φf(m，n)次幂。去除平地相位后的极化干涉数据对S1，hh(m，n)、S1，hv(m，n)、S1，vh(m，n)、S1，vv(m，n)、S′2，hh(m，n)、S′2，hv(m，n)、S′2，vh(m，n)、S′2，vv(m，n)包含的干涉相位仅与高度有关。
步骤S13在待滤波位置为中心的滤波窗口内，对极化干涉数据进行平滑处理，计算方位向标号为m、斜距向标号为n的位置滤波后的极化干涉数据T11(m，n)、T22(m，n)和T12(m，n) 其中上标H表示矩阵的共轭转置；括号符号<>表示平均。
表示由主天线数据构成的矢量 其中上标T表示矩阵转置；注脚数字1表示主天线数据；S1，hh(m，n)表示方位向标号为m、斜距向标号为n的位置的hh极化主天线数据；S1，vv(m，n)表示方位向标号为m、斜距向标号为n的位置的vv极化主天线数据；S1，hv(m，n)表示方位向标号为m、斜距向标号为n的位置的hv极化主天线数据；S1，vh(m，n)表示方位向标号为m、斜距向标号为n的位置的vh极化主天线数据。
符号
是由去除平地相位后的副天线数据S′2，hh(m，n)、S′2，hv(m，n)、S′2，vh(m，n)、S′2，vv(m，n)计算得到的，其计算式为 滤波后的T12(m，n)无需补偿平地相位，直接用于后续处理。
步骤B将极化分类技术应用于主、副天线数据平均结果，区分不同的土地类型，其具体处理步骤为 步骤S21计算主、副天线数据T11(m，n)与T22(m，n)的平均结果T(m，n) 步骤S22基于极化分类技术，从T(m，n)计算对应于单次散射、二次散射以及体散射强度S(m，n)、D(m，n)、V(m，n) S(m，n)＝T(m，n)(1，1) D(m，n)＝T(m，n)(2，2) (7) V(m，n)＝T(m，n)(3，3) 其中注脚(1，1)表示矩阵T(m，n)中第1行、第1列处数据，注脚(2，2)、(3，3)分别表示矩阵T(m，n)第2行、第2列处数据和第3行、第3列处数据。
步骤S23依次将单次散射、二次散射以及体散射强度分贝值LS(m，n)、LD(m，n)、LV(m，n) LS(m，n)＝10·log10(S(m，n)) LD(m，n)＝10·log10(D(m，n)) (8) LV(m，n)＝10·log10(V(m，n)) 均匀量化为蓝色、红色、绿色通道整数颜色值B(m，n)、R(m，n)、G(m，n)
其中符号
计算小于等于数据的最大整数，符号min表示指定的对应于颜色值0的数据，max表示指定的对应于颜色值255的数据。根据蓝色、红色、绿色通道的颜色值B(m，n)、R(m，n)、G(m，n)生成极化分类结果图。
步骤C根据极化分类的结果图，分析数据覆盖区域的土地类型复杂程度，并结合地面控制点的分布将其划分为土地类型较为简单且包含至少一个地面控制点的若干子块。
步骤D分别估计各子块地形干涉相位，且子块内单散射中心区域、多散射中心区域采用不同的极化干涉方法提取地形干涉相位。各子块计算地形干涉相位的步骤为 步骤S41标记子块内多散射中心区域和单散射中心区域，该步骤是为了在多散射中心区域、单散射中心应用不同方法提取地形干涉相位，从而满足不同土地类型地形干涉相位估计需求。多散射中心区域的判定原则为 V(m，n)≥AV (10) 其中AV表示判定阈值。非多散射中心区域标记为单散射中心区域。阈值AV由对应于极化分类结果图中某一多散射中心区域的体散射强度最小值确定。
步骤S42在单散射中心区域，应用极化干涉相干最优方法通过求解具有最大干涉相干系数的干涉相位估计地形干涉相位。其具体的计算步骤为 (1)基于特征分解，将矩阵T(m，n)分解为三个相互独立的部分 其中vi、
分别表示矩阵T(m，n)第i个特征值和对应的特征矢量(i＝1，2，3)。
(2)计算对应于上述三个相互独立部分的复干涉相干系数
(3)地形干涉相位φg(m，n)是对应于最大干涉相干系数的干涉相位 其中Arg()表示计算复数的幅角主值，||表示计算复数的模，
表示三组复干涉相干系数
之中具有最大模值的复干涉相干系数。
步骤S43利用指定极化干涉方法估计子块内多散射中心区域的地形干涉相位。当子块内多散射中心区域为森林或者冰雪区，指定相干散射模型方法计算地形干涉相位；当子块内多散射中心区域为城区，指定CAPON算法计算地形干涉相位。
基于相干散射模型方法利用极化干涉相干最优方法计算的三组复干涉相干系数，拟合复干涉相干系数所在直线，地形干涉相位可由该直线与单位圆交点确定。基于相干散射模型方法估计地形干涉相位的步骤为 (1)按照步骤S42所述相干最优方法计算三组复干涉相干系数
(i＝1，2，3)； (2)采用最小二乘方法，用三组复干涉相干系数拟合复平面内直线方程y＝a·x+b，直线方程得系数计算公式如下 其中x、y表示复数据的实部、虚部，a、b表示直线方程中的一次项、常数项系数，符号Re()表示复数的实部，Im()表示复数的虚部。
(3)计算直线方程y＝a·x+b与单位圆x2+y2＝0的两个交点(x1，y1)、(x2，y2)，其中x1、y1表示第一个交点的实部、虚部，x2、y2表示第二个交点的实部、虚部。将这两个交点表示为复数
其中φ1、φ2表示与第一个交点、第二各交点对应的复数的幅角主值； (4)计算hh+vv极化、hv+vh极化下干涉复相干系数
(15) (5)根据hv+vh极化状态的干涉相干系数距地形干涉相位最远的假设，确定地形干涉相位φg(m，n) (16) 其中Sign()表示计算数据的符号。当数据大于等于零时，Sign()的取值为1；当数据小于零，Sign()取值为-1。
CAPON算法利用分块矩阵以及特征分解技术，通过计算谱函数缝制求解地形干涉相位；地形干涉相位经过粗略峰搜索、精细谱峰搜索两阶段计算获得。具体的计算步骤为 (1)先利用T11(m，n)、T22(m，n)、T12(m，n)组合为6×6矩阵R6 (2)计算R6的逆矩阵
并将
划分为四个3×3矩阵
和
其中上标-1表示矩阵求逆。
(3)粗略搜索达到谱函数峰值处相位。粗略搜索是在(-π，π]均匀选取K个角度采样点 其中符号k为采样点标号，它是小于等于K的自然数。计算各个角度采样点φk谱函数P(φk) 其中Γmin()表示计算矩阵非负实数特征值中最小的特征值。比较搜索具有最大谱函数的角度采样点
其中kmax表示对应于最大谱函数的角度采样点标号。
(5)精细搜索达到谱函数峰值处相位，计算地形干涉相位。在粗略搜索到的对应于最大谱函数的角度采样点
及与
最近临近前后两采样点
范围内，再进行采样数为K的均匀精细采样。各个精细采样点φ′k 计算各个精细采样点的谱函数值，比较搜索具有最大谱函数的精细采样点。地形干涉相位φg(m，n)是使得谱函数P(φ′k)达到最大的精细角度采样点相位值。
步骤E由各子图块的地形干涉相位并行处理，分别按照常规干涉相位估计地形的方法，经过相位滤波、相位展开、平地相位补偿，地形高度估计，最终获得各子块的地形结果。其具体的计算步骤包含 步骤S51地形干涉相位φg(m，n)滤波 符号
表示滤波后的地形干涉相位。该步骤是为增强地形干涉相位的连续性，降低后续相位展开处理的难度。
步骤S52将滤波后的地形干涉相位进行相位展开处理。由于滤波后的地形干涉相位仍为相位主值，相位展开处理可恢复丢失的整数倍的2π信息 其中
表示相位展开后的地形干涉相位结果，2·Z·π表示待恢复的整数倍的2π信息。
步骤S53将平地相位补偿至展开后的地形干涉相位中
表示补偿了平地相位后的地形干涉相位结果，符号ni表示当前处理的第i个子块左上顶点像素的斜距标号。
步骤S54地形高度估计是利用展开后的地形干涉相位计算地形高度。地形高度估计包含如下步骤 (1)由补偿平地相位后的地形干涉相位计算主、副天线到目标点的斜距差Δr(m，n) (2)利用斜距差计算主天线视角θ(m，n)；
其中符号arccos()为反余弦函数；arctan()为反正切函数；Bh主、副天线之间的平均水平基线距；Bv是平均垂直基线距。
(3)计算地形高度htopo(m，n) htopo(m，n)＝H-(rn+(n+ni-2)rs)·cos(θ(m，n)) (28) 其中H表示主天线系统的绝对高度；rn为近距点斜距；rs是斜距采样间隔，符号ni表示当前处理的第i个子块左上顶点像素的斜距标号。
步骤F组合子图块的地形高度，按照各子块在原始数据中的位置，排列各子块的地形高度结果即获得组合成覆盖复杂土地类型复杂区域的完整的地形高度。
本发明提出的极化干涉合成孔径雷达数据的估计地形的软件提供了实现本发明方法中各处理步骤的代码段平地相位估计的代码段、去除平地相位对极化干涉数据影响的代码段、极化干涉数据平滑的代码段、极化分类的代码段、子块划分代码段、地形干涉相位估计代码段、相位滤波代码段、相位展开代码段、平地相位补偿代码段、地形高度估计代码段以及组合子块地形的代码段。本发明的极化干涉合成孔径雷达数据估计地形的软件，是中国科学院电子学研究所编写的，可在中国科学院电子学研究所微波成像技术国家级重点实验室购买。
平地相位估计代码段，完成极化干涉合成孔径雷达数据的估计地形高度方法步骤A中步骤S11的处理过程，计算各位置(m，n)的平地相位φf(m，n) 去除平地相位对极化干涉数据影响的代码段，完成极化干涉合成孔径雷达数据的估计地形高度方法步骤A中步骤S12的处理过程，去除平地相位对指定极化状态pq的各位置(m，n)的副天线数据的影响其中pq为hh、hv、vh或者vv。
极化干涉数据平滑的代码段，完完成极化干涉合成孔径雷达数据的估计地形高度方法步骤A中步骤S13的处理过程，该代码段在各位置(m，n)为中心的指定大小Mf×Nf的滤波窗口内计算T11(m，n)、T22(m，n)和T12(m，n)
其中m′、n′是指定范围内的整数。
极化分类的代码段，完成极化干涉合成孔径雷达数据的估计地形高度方法步骤B中步骤S21～步骤S23的处理过程，其中min、max是输入的对应于颜色值0、255的散射强度分贝值。
子块划分代码段，完成极化干涉合成孔径雷达数据的估计地形高度方法步骤C的处理过程通过设置各子块的左上角的位置(m1，nl)以及各子块大小Ml×Nl，计算各子块内数据T11(m，n)、T22(m，n)和T12(m，n)。
地形干涉相位估计代码段，完成极化干涉合成孔径雷达数据的估计地形高度方法步骤D中步骤S41～步骤S43处理过程其中判定阈值AV以及处理估计多散射中心区域的地形干涉相位的方法需预先设定。
相位滤波、相位展开、平地相位补偿、地形高度估计的代码段，完成极化干涉合成孔径雷达数据的估计地形高度方法步骤E中相位滤波、相位展开、平地相位补偿、地形高度估计处理过程。
组合子块地形的代码段，完成极化干涉合成孔径雷达数据的估计地形高度方法步骤F的处理过程按照各子块左上角位置(m1，nl)以及各子块大小Ml×Nl将各子块的地形高度结果排列为完整的地形高度结果。
下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。
本发明实例中使用的数据是德国ESAR系统于2003年10月16日19:30、19:11采集自Traunstein试验区L波段全极化数据构成的单基线极化干涉数据，这两组全极化数据的编号为03treesr0110x1和03treesr0109x1。图像尺寸为1500×1000个像素(方位向×斜距向)。该区域包含大片的森林，森林的平均高度约为40米，且森林和非森林区域交错分布。因而，由于该数据覆盖森林这一多散射中心区域，因而常规InSAR数据地形估计结果高于实际地形结果，出现明显的估计误差。由于森林与多散射中心区域交错分布，统一采用极化干涉相干最优方法或者基于相干散射模型的方法估计地形干涉相位不能满是实际土地类型地形干涉相位估计需求。
应用本发明所述的极化干涉合成孔径雷达数据估计地形高度的方法及其软件处理实例数据，最终获得该森林与分森林交错分布区域的地形信息。
实例中数据处理过程的六个步骤 步骤A极化干涉数据滤波，其包含以下步骤 步骤S11基于成像几何关系，从主、副天线之间的平均水平基线距Bh＝-4.3503316、平均垂直基线距Bv＝0.067152181、主天线相对参考平面高度h＝3024.33，近距点斜距rn＝3787.786309923、斜距采样点间隔rs＝1.49854，雷达信号波长λ＝0.23054400，计算平地相位φf(m，n) 其中m、n是方位向、斜距向位置标号，它们满足1≤m≤1500，1≤n≤1000。
步骤S12计算去除平地相位后的副天线数据，由下列公式
(2)得出 (30) 其中e表示自然对数，符号j表示虚数单位，符号
表示自然对数e的-j·φf(m，n)次幂；注脚数字2表示副天线数据。
步骤S13在待滤波位置为中心的7×7矩形滤波窗口内，计算方位向标号为m、斜距向标号为n的位置滤波后的极化干涉数据T11(m，n)、T22(m，n)和T12(m，n)
其中滤波窗口大小指定为Mf＝7、Nf＝7；m′、n′是指定范围内的整数，上标H表示矩阵的共轭转置。
是由主天线数据构成的矢量
其中上标T表示矩阵转置；注脚数字1表示主天线数据；S1，hh(m，n)表示方位向标号为m、斜距向标号为n的位置的hh极化主天线数据；S1，vv(m，n)表示方位向标号为m、斜距向标号为n的位置的vv极化主天线数据；S1，hv(m，n)表示方位向标号为m、斜距向标号为n的位置的hv极化主天线数据；S1，vh(m，n)表示方位向标号为m、斜距向标号为n的位置的vh极化主天线数据。
是由去除平地相位后的副天线数据S′2，hh(m，n)、S′2，hv(m，n)、S′2，vh(m，n)、S′2，vv(m，n)计算得到的，其计算式为 步骤B将极化分类技术应用于主、副天线数据平均结果，区分不同的土地类型，其具体处理步骤为 步骤S21计算主、副天线数据T11(m，n)与T22(m，n)的平均结果T(m，n) 步骤S22基于极化分类技术，从T(m，n)计算对应于单次散射、二次散射以及体散射强度S(m，n)、D(m，n)、V(m，n) S(m，n)＝T(m，n)(1，1) D(m，n)＝T(m，n)(2，2) (7) V(m，n)＝T(m，n)(3，3) 其中注脚(1，1)表示矩阵T(m，n)中第1行、第1列处数据，注脚(2，2)、(3，3)分别表示矩阵T(m，n)第2行、第2列处数据和第3行、第3列处数据。
步骤S23依次将单次散射、二次散射以及体散射强度分贝值LS(m，n)、LD(m，n)、LV(m，n) LS(m，n)＝10·log10(S(m，n)) LD(m，n)＝10·log10(D(m，n)) (8) LV(m，n)＝10·log10(V(m，n)) 均匀量化为蓝色、红色、绿色通道整数颜色值B(m，n)、R(m，n)、G(m，n)
其中量化参数指定为min＝20、max＝62，符号
计算小于等于数据的最大整数。根据蓝色、红色、绿色通道的颜色值B(m，n)、R(m，n)、G(m，n)生成图3所示的极化分类结果图。
步骤C根据极化分类的结果图，分析数据覆盖区域的土地类型复杂程度。结合地面控制点的分布，将数据划分为R(m，n)图4所示的六个子块F1、F2、F3、F4、F5、F6。F1是方位向标号在1～360，斜距向标号为1～450的数据构成的360×450的矩形子块；F2是方位向标号在361～1500，斜距向标号为1～180的数据构成的1140×180的矩形子块；F3是方位向标号在361～910，斜距向标号为181～450的数据构成的550×270的矩形子块；F4是方位向标号在911～1500，斜距向标号为181～520的数据构成的590×340的矩形子块；F5是方位向标号在911～1500，斜距向标号为521～1000的数据构成的590×480的矩形子块；F6是方位向标号在1～910，斜距向标号为451～1000的数据构成的910×550的矩形子块。
步骤D分别估计F1、F2、F3、F4、F5、F6子块地形干涉相位，且子块内单散射中心区域、多散射中心区域采用不同的极化干涉方法提取地形干涉相位。子块F1、F2、F3、F4、F5、F6计算地形干涉相位的步骤为 步骤S41标记子块内多散射中心区域和单散射中心区域，该步骤是为了在多散射中心区域、单散射中心应用不同方法提取地形干涉相位，从而满足不同土地类型地形干涉相位估计需求。多散射中心区域的判定原则为 V(m，n)≥AV (10) 其中判定阈值设置为AV＝15800；非多散射中心区域标记为单散射中心区域。
步骤S42在单散射中心区域，应用极化干涉相干最优方法通过求解具有最大干涉相干系数的干涉相位估计地形干涉相位。其具体的计算步骤为 (1)基于特征分解，将矩阵T(m，n)分解为三个相互独立的部分 其中vi、
分别表示矩阵T(m，n)第i个特征值、特征矢量(i＝1，2，3)。
(2)计算对应于上述三个相互独立部分的复干涉相干系数
(3)地形干涉相位φg(m，n)是对应于最大干涉相干系数的干涉相位 其中Arg()表示计算复数的幅角主值，||表示计算复数的模，
表示三组复干涉相干系数
之中具有最大模值的复干涉相干系数。
步骤S43由于各子块内多散射中心区域均为森林，因而选定相干散射模型方法计算多散射中心区域地形干涉相位。
基于相干散射模型方法估计地形干涉相位的步骤为 (1)按照步骤S42所述相干最优方法计算三组复干涉相干系数
(i＝1，2，3)； (2)采用最小二乘方法，用三组复干涉相干系数拟合复平面内直线方程y＝a·x+b，直线方程得系数计算公式如下 其中x、y表示复数据的实部、虚部，a、b表示直线方程中的一次项、常数项系数，符号Re()表示复数的实部，Im()表示复数的虚部。
(3)计算直线方程y＝a·x+b与单位圆x2+y2＝0的两个交点(x1，y1)、(x2，y2)，其中x1、y1表示第一个交点的实部、虚部，x2、y2表示第二个交点的实部、虚部。将这两个交点表示为复数
其中φ1、φ2表示与第一个交点、第二各交点对应的复数的幅角主值； (4)计算hh+vv极化、hv+vh极化下干涉复相干系数
(15) (5)根据hv+vh极化状态的干涉相干系数距地形干涉相位最远的假设，确定地形干涉相位φg(m，n) (16) 其中Sign()表示计算数据的符号。当数据大于等于零时，Sign()的取值为1；当数据小于零，Sign()取值为-1。
子块F1的地形干涉相位结果如图5a所示；子块F2的地形干涉相位结果如图5b所示；子块F3的地形干涉相位结果如图5c所示；子块F4的地形干涉相位结果如图5d所示；子块F5的地形干涉相位结果如图5e所示；子块F6的地形干涉相位结果如图5f所示。
步骤E子块F1、F2、F3、F4、F5、F6的地形干涉相位并行处理，分别按照常规干涉相位估计地形的方法，经过相位滤波、相位展开、平地相位补偿，地形高度估计，最终获得各子块的地形结果。子块F1、F2、F3、F4、F5、F6的地形计算步骤为 步骤S51采用7×7矩形滤波窗口，对地形干涉相位滤波。地形干涉相位滤波结果的是在7×7矩形滤波窗口内应用公式
(23)获得，具体的计算公式为 符号
表示滤波后的地形干涉相位。子图F1中方位向、斜距向标号的变化范围为1≤m≤360，1≤n≤450；子图F2中方位向、斜距向标号的变化范围1≤m≤1140，1≤n≤180；子图F3中方位向、斜距向标号的变化范围为1≤m≤550，1≤n≤270；子图F4中方位向、斜距向标号的变化范围1≤m≤590，1≤n≤340；子图F5中方位向、斜距向标号的变化范围为1≤m≤590，1≤n≤480；子图F2中方位向、斜距向标号的变化范围1≤m≤910，1≤n≤550。
步骤S52将滤波后的地形干涉相位进行相位展开处理。采用枝切线法进行在干涉相位连续区域内进行相位积分，恢复丢失的整数倍的2π信息 其中
表示相位展开后的地形干涉相位结果。子图F1中方位向、斜距向标号的变化范围为1≤m≤360，1≤n≤450；子图F2中方位向、斜距向标号的变化范围1≤m≤1140，1≤n≤180；子图F3中方位向、斜距向标号的变化范围为1≤m≤550，1≤n≤270；子图F4中方位向、斜距向标号的变化范围1≤m≤590，1≤n≤340；子图F5中方位向、斜距向标号的变化范围为1≤m≤590，1≤n≤480；子图F2中方位向、斜距向标号的变化范围1≤m≤910，1≤n≤550。
步骤S53将平地相位补偿至展开后的地形干涉相位中。
将各子块左上顶点像素的斜距标号代入公式
(25)，计算各子块补偿平地相位后的地形干涉相位结果。
由于子块F1平地相位补偿结果为 其中子块F1方位向、斜距向标号的变化范围为1≤m≤360，1≤n≤450； 由于子块F2平地相位补偿结果为 其中子块F2方位向、斜距向标号的变化范围为1≤m≤1140，1≤n≤180； 由于子块F3平地相位补偿结果为 其中子块F3方位向、斜距向标号的变化范围为1≤m≤550，1≤n≤270； 由于子块F4平地相位补偿结果为 其中子块F4方位向、斜距向标号的变化范围为1≤m≤590，1≤n≤340； 由于子块F5平地相位补偿结果为 其中子块F5方位向、斜距向标号的变化范围为1≤m≤590，1≤n≤480； 由于子块F6平地相位补偿结果为 其中子块F6方位向、斜距向标号的变化范围为1≤m≤910，1≤n≤550； 步骤S54地形高度估计是利用展开后的地形干涉相位计算地形高度。地形高度估计包含如下步骤 (1)由补偿平地相位后的地形干涉相位、雷达信号波长λ＝0.23054400计算主、副天线到目标点的斜距差Δr(m，n) (2)将利用斜距差、平均水平基线距Bh＝-4.3503316、平均垂直基线距Bv＝0.067152181计算主天线视角θ(m，n)；
其中符号arccos()为反余弦函数；arctan()为反正切函数；Bh主、副天线之间的平均水平基线距；Bv是平均垂直基线距。
(3)由主天线系统的绝对高度H＝3674.33、近距点斜距rn＝3787.786309923、斜距采样点间隔rs＝1.49854代入公式htopo(m，n)＝H-(rn+(n+ni-2)·rs)·cos(θ(m，n))(28)计算各子块地形高度htopo(m，n)。
由于子块F1地形高度计算公式为 htopo(m，n)＝H-(rn+(n-1)·rs)·cos(θ(m，n)) (36) 其中子块F1方位向、斜距向标号的变化范围为1≤m≤360，1≤n≤450； 由于子块F2地形高度计算公式为 htopo(m，n)＝H-(rn+(n-1)·rs)·cos(θ(m，n)) (37) 其中子块F2方位向、斜距向标号的变化范围为1≤m≤1140，1≤n≤180； 由于子块F3地形高度计算公式为 htopo(m，n)＝H-(rn+(n+179)·rs)·cos(θ(m，n)) (38) 其中子块F3方位向、斜距向标号的变化范围为1≤m≤550，1≤n≤270； 由于子块F4地形高度计算公式为 htopo(m，n)＝H-(rn+(n+179)·rs)·cos(θ(m，n))(39) 其中子块F4方位向、斜距向标号的变化范围为1≤m≤590，1≤n≤340； 由于子块F5地形高度计算公式为 htopo(m，n)＝H-(rn+(n+519)·rs)·cos(θ(m，n))(40) 其中子块F5方位向、斜距向标号的变化范围为1≤m≤590，1≤n≤480； 由于子块F6地形高度计算公式为 htopo(m，n)＝H-(rn+(n+449)·rs)·cos(θ(m，n))(41) 其中子块F6方位向、斜距向标号的变化范围为1≤m≤910，1≤n≤550； 子块F1的地形高度结果如图6a所示；子块F2的地形高度结果如图6b所示；子块F3的地形高度结果如图6c所示；子块F4的地形高度结果如图6d所示；子块F5的地形高度结果如图6e所示；子块F6的地形高度结果如图6f所示。
步骤F组合子图块的地形高度，按照如图3所示的各子块在原始数据中的位置，排列各子块的地形高度结果即获得组合成覆盖复杂土地类型复杂区域的完整的地形高度。完整的地形高度结果如图7a所示。图7b是HH极化的干涉数据估计的地形与本发明的方法估计的地形在指定方位标号为910的剖面上对比结果。对比结果表明，本发明所提出的极化干涉合成孔径雷达数据估计地形的方法多散射中心区的结果在低于HH极化的干涉数据估计结果，即说明本发明方法在多散射中心的地形高度更接近于实际地形。对比结果还表明，在单散射中心区域本发明方法地形估计与HH极化的干涉数据估计结果十分接近，说明本发明方法在单散射中心区域克获得准确的地形估计。实例数据处理结果证实本发明方法即能在多散射中心区域又能在单散射中心区域获得可靠的地形信息，因而该发明方法于复在复杂土地类型区域的地形估计中有巨大的应用潜力。
本发明方法充分利用极化干涉数据的优势，引入极化分类技术将复杂土地类型区域划分为若干土地类型构成简单的子块分别处理。本发明方法不但可降低地形估计的难度，还可根据不同土地类型选择适合的极化干涉技术估计地形干涉相位从而获取可靠的地形高度信息。
参考文献Cloude S.R.and Papathanassiou K.P.(1998)，Polarimetric SARInterferometry[J]，IEEE Trans.on Geoscience and Remote Sensing，36(5)1551-1565.Brandfass M.(2002)，Generation of bald earth digital elevation modelsas applied to polarimetric SAR Interferomety，IGARSS 2002，21014-1016.Mounira Souissi Ouarzeddine，(2002).Generation of Digital Modelsusing Polarimetric SAR Interferometry，International Institute forGeo-information Science and Earth Observation，Master Degree.Kasilingam D.，Nomula M.，Cloude S.(2002)，A technique forremoving vegetation bias from polarimetric SAR interferometry，IGARSS 2002，21017-1019.Papathanassiou K.P.and Cloude S.R.(2001)，Single-BaselinePolarimetric SAR Interferometry，IEEE Trans.on Geoscience andRemote Sensing，39(11)2352-2363.Guillaso S.，et al.(2003)，Urban area analysis based onESPRIT/MUSIC methods using polarimetric interferometric SAR，2ndGRSS/ISPRS JOINT WORKSHOP ON REMOTE SENSING ANDDATA FUSION OVER UBRAN AREAS，77-81.
权利要求
1.一种基于极化干涉合成孔径雷达数据的地形估计方法，其特征在于，包括步骤
A.对极化干涉数据滤波；B.利用极化分类技术区分不同的土地类型；C.将数据划分为土地类型构成简单的若干子块；D.分别估计各子块的地形干涉相位；E.并行处理各子块地形干涉相位，估计各子块地形；F.组合各子块地形，计算完整的地形结果。
2.如权利要求1所述的方法，其特征在于所述A步，对极化干涉数据滤波，滤波后的数据无需补偿平地相位，极化干涉数据滤波抑制了噪声对地形干涉相位估计的不利影响，其处理包括以下步骤
步骤S11基于成像几何关系，从主、副天线之间的平均水平基线距Bh、平均垂直基线距Bv、主天线相对参考平面高度h，近距rn、斜距采样间隔rs，雷达信号波长λ，估计方位向标号为m、斜距向标号为n的位置的平地相位φf(m，n)
其中π表示180°角对应的弧度值，符号·表示乘运算，m、n是方位向、斜距向位置标号，它们是自然数且满足1≤m≤M，1≤n≤N；符号M、N表示极化干涉数据在方位向、斜距向总的像素个数；估计平地相位是为抑制相邻像素之间平地相位差异对极化干涉数据滤波造成的不利影响；
步骤S12去除平地相位对极化干涉数据的影响
从而增强相邻像素之间极化干涉数据的连续性；其中注脚数字2表示副天线数据；p、q表示接收、发送信号的极化方式，是水平极化方式h或者垂直极化方式V；S2，pq(m，n)表示方位向标号为m、斜距向标号为n的位置的pq极化副天线数据；S′2，pq(m，n)表示去除平地相位影响后的方位向标号为m、斜距向标号为n的位置的pq极化副天线数据；e表示自然对数，符号j表示虚数单位，符号
表示自然对数e的-j·φf(m，n)次幂；
步骤S13在待滤波位置为中心的滤波窗口内，对极化干涉数据进行平滑处理，计算方位向标号为m、斜距向标号为n的位置滤波后的极化干涉数据T11(m，n)、T22(m，n)和T12(m，n)
其中上标H表示矩阵的共轭转置；括号符号< >表示平均。
表示由主天线数据构成的矢量
其中上标T表示矩阵转置；注脚数字1表示主天线数据；S1，hh(m，n)表示方位向标号为m、斜距向标号为n的位置的hh极化主天线数据；S1，hv(m，n)表示方位向标号为m、斜距向标号为n的位置的hv极化主天线数据；S1，vh(m，n)表示方位向标号为m、斜距向标号为n的位置的vh极化主天线数据；S1，vv(m，n)表示方位向标号为m、斜距向标号为n的位置的vv极化主天线数据；
符号
是由去除平地相位后的副天线数据S′2，hh(m，n)、S′2，hv(m，n)、S′2，vh(m，n)、S′2，vv(m，n)构成，其表达式为
滤波后的T12(m，n)无需补偿平地相位，直接用于后续处理。
3.如权利要求1所述的方法，其特征在于所述B步，将极化分类技术应用于主、副天线数据平均结果来区分不同的土地类型，利用极化分类是为分析数据覆盖的土地类型复杂程度，为把数据划分为若干土地类型构成简单的子块提供依据，其具体处理步骤如下
步骤S21计算主、副天线全极化数据T11(m，n)与T22(m，n)的算术平均
其中T(m，n)表示主、副天线全极化数据平均结果；
步骤S22从T(m，n)计算单次散射、二次散射以及体散射强度S(m，n)、D(m，n)、V(m，n)
S(m，n)＝T(m，n)(1，1)
D(m，n)＝T(m，n)(2，2)(7)
V(m，n)＝T(m，n)(3，3)
其中注脚(1，1)表示矩阵T(m，n)中第1行、第1列处数据，注脚(2，2)、(3，3)分别表示矩阵T(m，n)第2行、第2列处数据和第3行、第3列处数据；
步骤S23依次将单次散射、二次散射以及体散射强度的分贝值LS(m，n)、LD(m，n)、LV(m，n)
LS(m，n)＝10·log10(S(m，n))
LD(m，n)＝10·log10(D(m，n))(8)
LV(m，n)＝10·log10(V(m，n))
均匀量化为蓝色、红色、绿色通道整数颜色值B(m，n)、R(m，n)、G(m，n)
其中log10( )表示计算10为底的数据的对数值，符号
计算小于等于数据的最大整数，符号min表示指定的对应于颜色值0的数据，max表示指定的对应于颜色值255的数据，根据蓝色、红色、绿色通道的颜色值B(m，n)、R(m，n)、G(m，n)生成极化分类结果图。
4.如权利要求1所述的方法，其特征在于所述C步，将数据划分为土地类型构成简单的若干子块，该步骤降低了土地类型构成的复杂度、降低了地形干涉相位估计难度，并兼顾不同土地类型区域采用不同方法估计地形干涉相位需求；根据极化分类结果图，分析土地类型构成的复杂度，结合地面控制点的分布将其划分为土地类型构成简单且至少包含一个地面控制点的若干子块。
5.如权利要求1所述的方法，其特征在于所述D步，分别估计各子块地形干涉相位，且子块内单散射中心区域、多散射中心区域采用不同的极化干涉方法提取地形干涉相位，各子块数据计算地形干涉相位的步骤为
步骤S41标记子块内多散射中心区域和单散射中心区域，该步骤是在多散射中心区域、单散射中心区域应用不同方法提取地形干涉相位，从而满足不同土地类型地形干涉相位估计需求；多散射中心区域判定原则为
V(m，n)≥AV(10)
其中AV表示判定阈值，非多散射中心区域标记为单散射中心区域，阈值AV由对应于极化分类结果图中某一多散射中心区域的体散射强度最小值确定；
步骤S42在单散射中心区域，应用极化干涉相干最优方法通过求解具有最大干涉相干系数的干涉相位估计地形干涉相位；
步骤S43利用指定极化干涉方法估计子块内多散射中心区域的地形干涉相位，当子块内多散射中心区域为森林或者冰雪区，指定相干散射模型方法计算地形干涉相位；当子块内多散射中心区域为城区，指定CAPON算法计算地形干涉相位。
6.如权利要求5所述的方法，其特征在于所述步骤S42中的极化干涉相干最优方法，是利用T(m，n)提取散射中心信息，将具有最大干涉相干系数的干涉相为作为地形干涉相位，具体的计算步骤为
步骤S61基于特征分解，将矩阵T(m，n)分解为三个相互独立的部分
其中vi、
分别表示矩阵T(m，n)第i个特征值和对应的特征矢量(i＝1，2，3)；
步骤S62计算对应于上述三个相互独立部分的复干涉相干系数
步骤S63地形干涉相位φg(m，n)是对应于最大干涉相干系数的干涉相位
其中Arg( )表示计算复数的幅角主值，| |表示计算复数的模，
表示三组复干涉相干系数
之中具有最大模值的复干涉相干系数。
7.如权利要求5或6所述的方法，其特征在于所述步骤S43中的相干散射模型方法，是利用极化干涉相干最优方法获取的三组复干涉相干系数，拟合复干涉相干系数所在直线，地形干涉相位由该直线与单位圆交点确定，其具体的估计步骤为
步骤S71利用极化干涉相干最优方法估计三组复干涉相干系数
(i＝1，2，3)；
步骤S72采用最小二乘方法，用三组复干涉相干系数拟合复平面内直线方程y＝a·x+b，具体的计算方法为
其中x、y表示复数据的实部、虚部，a、b表示直线方程中的一次项、常数项系数，符号Re( )表示复数的实部，Im( )表示复数的虚部；
步骤S73计算直线y＝a·x+b与单位圆x2+y2＝1的两个交点(x1，y1)、(x2，y2)，其中x1、y1表示第一个交点的实部、虚部，x2、y2表示第二个交点的实部、虚部。将这两个交点表示为复数
其中φ1、φ2表示与第一个交点、第二各交点对应的复数的幅角主值；
步骤S74计算hh+vv极化、hv+vh极化下干涉复相干系数
步骤S75根据hv+vh极化状态的干涉相干系数距地形干涉相位最远的假设，确定地形干涉相位φg(m，n)
其中Sign( )表示计算数据的符号，当数据大于等于零时，Sign( )的取值为1；当数据小于零，Sign( )取值为-1。
8.如权利要求5所述的方法，其特征在于所述步骤S43中的CAPON算法，是利用分块矩阵以及特征分解技术，通过计算谱函数峰值求解地形干涉相位；地形干涉相位估计经过粗略峰搜索、精细谱峰搜索两阶段计算获得，具体的计算步骤为
步骤S81先利用T11(m，n)、T22(m，n)、T12(m，n)组合为6×6矩阵R6
步骤S82计算R6的逆矩阵
，并将
划分为四个3×3矩阵
和
其中上标-1表示矩阵求逆；
步骤S83粗略搜索达到谱函数峰值处相位，粗略搜索是在(-π，π]均匀选取K个角度采样点
其中符号k为采样点标号，它是小于等于K的自然数；计算各个角度采样点φk谱函数P(φk)
其中Γmin( )表示计算矩阵非负实数特征值中最小的特征值，比较搜索具有最大谱函数的角度采样点
，其中kmax表示对应于最大谱函数的角度采样点标号；
步骤S84精细搜索达到谱函数峰值处相位，计算地形干涉相位；在粗略搜索到的对应于最大谱函数的角度采样点
及与
最近临近前后两采样点
范围内，再进行采样数为K的均匀精细采样，各个精细采样点φ′k
计算各个精细采样点的谱函数值P(φ′k)，比较搜索具有最大谱函数的精细采样点；地形干涉相位φg(m，n)是使得谱函数P(φ′k)达到最大的精细角度采样点相位值。
9.如权利要求1所述的方法，其特征在于所述E步，并行处理各子块地形干涉相位，估计各子块的地形，各子块并行处理缩短了地形估计时间、提高地形估计效率；估计各子块的地形，按照常规干涉相位估计地形的方法，分别经过相位滤波、相位展开、平地相位补偿、地形高度估计。
10.如权利要求1所述的方法，其特征在于所述F步，组合子块地形，计算完整的地形结果，按照各子块在原始数P(φ′k)据中的位置，排列各子块的地形结果，从而组合成完整的地形结果。
11.一种用于如权利要求1所述的方法的极化干涉合成孔径雷达数据估计地形的软件，该软件是由记录在计算机刻度介质上的程序代码构成，该程序代码用于实现极化干涉合成孔径雷达数据估计地形的方法，其特征在于该程序代码包括
-用于完成步骤S11平地相位估计的代码段；
-用于完成步骤S12去除平地相位对极化干涉数据影响的代码段；
-用于完成步骤S13极化干涉数据平滑的代码段；
-用于完成极化分类的代码段；
-用于完成子块划分代码段；
-用于完成地形干涉相位估计代码段；
-用于完成相位滤波、相位展开、平地相位补偿、地形高度估计的代码段；
-用于完成组合子块地形的代码段；
该极化干涉合成孔径雷达数据估计地形的软件可在中国科学院电子学研究所微波成像技术国家级重点实验室购买。
12.如权利要求11所述的软件，其特征在于所述平地相位估计代码段，是在A步中的步骤S11的处理过程，计算各位置(m，n)的平地相位φf(m，n)。
13.如权利要求11所述的软件，其特征在于所述去除平地相位对极化干涉数据影响的代码段，是在A步中步骤S12的处理过程，去除平地相位对指定极化状态pq的各位置(m，n)的副天线数据的影响其中pq为hh、hv、vh或者vv。
14.如权利要求11所述的软件，其特征在于所述极化干涉数据平滑的代码段，是在A步中步骤S13的处理过程，该代码段在各位置(m，n)为中心的指定大小Mf×Nf的滤波窗口内计算T11(m，n)、T22(m，n)和T12(m，n)
其中m′、n′是指定范围内的整数。
15.如权利要求11所述的软件，其特征在于所述极化分类的代码段，是在B步中步骤S21～步骤S23的处理过程，其中min、max是输入的对应于颜色值0、255的散射强度分贝值。
16.如权利要求11所述的软件，其特征在于所述子块划分代码段，是在C步的处理过程通过设置各子块的左上角的位置(ml，nl)以及各子块大小Ml×Nl，计算各子块内数据T11(m，n)、T22(m，n)和T12(m，n)。
17.如权利要求11所述的软件，其特征在于所述地形干涉相位估计代码段，是在D步的步骤S41～步骤S43处理过程其中判定阈值AV以及处理估计多散射中心区域的地形干涉相位的方法需预先设定。
18.如权利要求11所述的软件，其特征在于所述相位滤波、相位展开、平地相位补偿、地形高度估计的代码段，是在步骤E中相位滤波、相位展开、平地相位补偿、地形高度估计处理过程。
19.如权利要求11所述的软件，其特征在于所述组合子块地形的代码段，是在F步的处理过程按照各子块左上角位置(ml，nl)以及各子块大小Ml×Nl将各子块的地形高度结果排列为完整的地形高度结果。
全文摘要
本发明公开了一种极化干涉合成孔径雷达数据估计地形的方法及其软件，涉及合成孔径雷达，包括步骤A.对极化干涉数据滤波；B.利用极化分类技术区分不同的土地类型；C.将数据划分为土地类型构成简单的若干子块；D.分别估计各子块的地形干涉相位；E.并行处理各子块地形干涉相位，估计各子块地形；F.组合各子块地形，计算完整的地形结果。本发明方法及软件引入极化分类技术，将数据划分成土地类型构成简单的若干子块处理，根据子块包含的土地类型差异选用不同的极化干涉方法估计地形干涉相位，再估计各子块的地形。本发明方法既降低土地类型构成复杂的区域地形估计的难度，又提高该类区域地形估计效率。
文档编号G01S13/90GK101770026SQ20091007656
公开日2010年7月7日 申请日期2009年1月7日 优先权日2009年1月7日
发明者洪文, 白璐, 曹芳 申请人:中国科学院电子学研究所