CN105005047B - 后向散射优化的森林复杂地形校正及树高反演方法、系统 - Google Patents

Abstract

本发明提供后向散射优化的森林复杂地形校正及树高反演方法、系统，地形校正包括针对森林场景的极化干涉SAR图像，根据森林场景的两层散射特点，构造两层结构的相干模型，并对地形引起的参数变化进行校正；结合地形仿真和植被仿真，建立植被与地形之间的定量描述关系；最后结合距离多普勒几何校正方式和归一化辐射校正方式进行校正；树高反演包括针对不同树种之间的结构和分布的差异，提出树种多样性因子，对不同的树种的植被高度分别进行估计然后获得整个森林场景的树高反演结果。采用本发明可以有效的实现森林复杂场景的地形效应分析与补偿；此外，结合树种多样性因子，提高植被高度的反演精度。

Description

后向散射优化的森林复杂地形校正及树高反演方法、系统

技术领域

本发明属于图像处理技术领域，特别涉及一种将后向散射优化用于极化干涉SAR森林复杂地表效应分析及树高反演的技术方案。

背景技术

森林在全球水文、生态、碳循环及气候变化中起着重要作用，森林的类别、林分结构、高度及生物量等参数是林业资源信息调查中的主要参数。随着遥感技术的发展，具有极化和干涉技术优点的极化干涉SAR技术，以其独特的全天候、低成本的优势，逐渐成为森林资源调查中一种不可替代的技术，在森林制图、特别是森林参数的定量反演中，发挥着越来越重要的作用。基于极化干涉SAR的森林分类、森林高度反演模型、生物量估计模型的研究是极化干涉SAR森林遥感研究的核心问题，对其开展研究有着重要的意义。

在森林垂直结构参数提取中，由于坡度倾斜，随机起伏以及局部地形破碎等复杂地形，除了在SAR图像上产生透视收缩、叠掩和阴影以外，地形还引起了局部入射角和后向散射系数的变化，从而使得森林结构模型的精度受到影响，从而大大的限制了参数反演的精度。

目前常用的校正方法包括几何与辐射两类，在几何方面，校正方法根据校正变换模型的不同又可以分为多项式校正法、共线方程校正法以及基于SAR成像原理的距离多普勒模型校正法三大类：1)多项式校正法，该方法是一种比较传统的校正方法，在SAR图像校正理论发展的初级阶段，被应用到SAR图像的几何校正中。这种方法回避成像的空间几何过程，而直接对图像变形的本身进行数学模拟，原理比较直观，并且计算较为简单，适合于地形相对平坦的情况。多项式的系数一般利用已知控制点的坐标数值按最小二乘法原理求解，多项式校正法的精度与地面控制点的精度、分布和数量及校正范围有关。采用多项式校正法的优点是能保证政府图像变换后总误差最小，但不能保证各局部的精度完全一致。共线共线方程法是基于传感器的成像方程来建立模型的一种方法，比多项式校正法在理论上更严密，因为它是建立在恢复实际成像条件的基础之上的。2)共线方程法，又被称为数字微分法，它是建立在对传感器成像时的位置和姿态进行模拟和结算的基础上，即成像瞬间的像点与其相应地物点应位于通过传感器中心的一根直线上，是对成像空间的几何描述。因为在共线方程中可考虑高程的影响，在地形起伏影响比较大的情况下，这种方法比多项式校正法精度要高，但是必须首先生成该影响范围内的DEM，否则在缺乏数字地面高程的情况下，影像的校正将变得复杂而困难。3)基于SAR本身成像原理的校正方法，该方法最初是F.Leberal针对机载SAR提出了斜距方程和零多普勒方程，但是它在针对星载SAR图像处理的时候会有一些不足，如星载SAR的多普勒一般不为零、距离变化较大等，所以需要另外建立适合星载SAR的成像模型。

基于DME模拟SAR影像进行校正的方法需要利用SAR成像参数和DME数据建立模拟图像，然后将模拟图像与真实(待校正)SAR图像进行配准，从而建立真实SAR成像到DEM坐标的变换公式。Gunnter Schreier于1993年出版的专著中介绍了“椭球表面地理编码”(Geocoded Ellipsoid Corrected，GEC)的概念。该方法将轨道模型采用多项式进行参数化描述，距离向方程、多普勒频率方程也都采用多项式进行描述，将地球表面看成一个具有固定高程的椭球面，而不考虑地形的起伏变化。模型参数都从雷达处理参数中提取，利用这些参数建立RD模型，采用间接法进行重采样定位处理。此外，基于“地形校正地理编码”(Geocoded Terrain corrected，GTC)方法，不需要对轨道方程进行参数化，而是采用插值方法计算得到任意时刻的卫星位置和速度矢量。由地面坐标计算影像坐标的方法采用基于多普勒方程的迭代方法，这和GEC方法有所不同。GTC利用DEM描述地球表面，因此地面坐标是精确已知的。这里假设卫星轨道和其它成像参数是正确的，不需要进行控制点对模型参数进行优化处理。

近年来，随着电磁波辐射传输模型的发展，基于辐射归一化的校正方法得到了广泛的关注，该方法通过建立后向散射系数、局部入射角和有效散射面积之间的关系模型，通过寻找与斜距方向垂直的散射面的稳定后向散射估计结果，可以实现地形参数的准确估计。此外，通过森林两层相干模型，求解最佳极化状态和相干相位，也可以有效的解决地形效应的影响。

然而，随着遥感技术的发展，当前基于几何插值的方法由于精度有限，难以达到定量遥感的反演需求，辐射校正中，地形因素往往是较平坦的或者仅考虑了地形倾斜的影响，而对于景观破碎和起伏多变的环境研究较少，同时，由于现有的地形多变，未能结合仿真数据对真实森林场景进行精确的模拟分析。此外，在实际的树高参数估计中，由于不同树种在分布和结构上存在一定的差异，其对入射角的变化存在差异，不仅直接影响了辐射校正的精度，而且采用统一的模型进行树高反演会存在一定的误差。本领域尚未有能够解决问题的相关技术出现。

发明内容

针对现有技术缺陷，本发明提供一种基于后向散射优化的极化干涉SAR森林复杂地表地形校正及树高反演技术方案。

本发明技术方案提供一种后向散射优化的森林复杂地形校正方法，包括以下步骤：

步骤1，针对森林场景的极化干涉SAR图像，根据森林场景的两层散射特点，构造两层结构的相干模型如下，

其中，ω表示某一给定的极化状态，γ(ω)为相应的干涉相干系数，表示极化方位角，i是虚数单位；

m(ω)为地体辐射比，定义为

其中，m_G(ω)和m_V(ω)分别表示地表层的表面散射和树冠层的体散射的幅度，

γ_V为纯体散射，表达式如下，

其中，σ表示后向散射系数，θ表示入射角，k_z表示有效垂直波数，z表示垂直方向高度，h_V表示植被高度；

步骤2，在步骤1所得相干模型基础上，对地形引起的参数变化进行校正，实现如下，首先进行散射机制分解，

对于表面散射，采用Freeman-Durden方法进行散射机制分类；

对于体散射，选取满足下式的投影向量w₁和w₂，实现体散射机制分离，

其中，T_v为体散射单元，m_v表示每个散射体中体散射的幅度，η表示粒子的平均形状，C_v为系数，定义为

其中，θ₀表示局部入射角，e为自然常数，z'表示地面参考高度，j是虚数单位；

然后，针对地形引起的后向散射相位的变化，转换成拉格朗日优化的形式，通过求解矩阵特征值得到地形校正的结果，

F^-1Aw＝λw

其中，F、A为正半定Hermitian矩阵，λ为拉格朗日因子，w表示后向散射相应的归一化的权值向量；

步骤3，在步骤2所得散射机制分解结果的基础上，结合地形仿真和植被仿真，建立植被与地形之间的定量描述关系；

步骤4，根据步骤3所得结果，结合距离多普勒几何校正方式和归一化辐射校正方式进行校正，

所述归一化辐射校正方式实现如下，

cos(ψ)＝sin(θ₀)·cos(u)+cos(θ₀)·sin(u)·sin(v)＝A_area ^-1

其中，ψ为投影角，u和v分别表示传感器方向的地形坡度和方位角，A_area表示散射面积。

本发明还基于上述后向散射优化的森林复杂地形校正方法提供树高反演方法，针对不同树种之间的结构和分布的差异，提出树种多样性因子，对不同的树种的植被高度分别进行估计然后获得整个森林场景的树高反演结果，

树种多样性参数如下式所示，

其中，S表述树的种类数，N表示每一类树的数目，α表示多样性参数，则树高估计的表达式为，

其中，为体散射估计值，k_z表示有效垂直波数，φ₀表示干涉相位角。

本发明相应提供一种后向散射优化的森林复杂地形校正系统，包括以下模块，相干模型构建模块，用于针对森林场景的极化干涉SAR图像，根据森林场景的两层散射特点，构造两层结构的相干模型如下，

其中，ω表示某一给定的极化状态，γ(ω)为相应的干涉相干系数，表示极化方位角，i是虚数单位；

m(ω)为地体辐射比，定义为

其中，m_G(ω)和m_V(ω)分别表示地表层的表面散射和树冠层的体散射的幅度，

γ_V为纯体散射，表达式如下，

其中，σ表示后向散射系数，θ表示入射角，k_z表示有效垂直波数，z表示垂直方向高度，h_V表示植被高度；

地形变化校正模块，用于在相干模型基础上，对地形引起的参数变化进行校正，实现如下，首先进行散射机制分解，

对于表面散射，采用Freeman-Durden方法进行散射机制分类；

对于体散射，选取满足下式的投影向量w₁和w₂，实现体散射机制分离，

其中，T_v为体散射单元，m_v表示每个散射体中体散射的幅度，η表示粒子的平均形状，C_v为系数，定义为

其中，θ₀表示局部入射角，e为自然常数，z'表示地面参考高度，j是虚数单位；

然后，针对地形引起的后向散射相位的变化，转换成拉格朗日优化的形式，通过求解矩阵特征值得到地形校正的结果，

F^-1Aw＝λw

其中，F、A为正半定Hermitian矩阵，λ为拉格朗日因子，w表示后向散射相应的归一化的权值向量；

仿真关系模块，用于在散射机制分解结果的基础上，结合地形仿真和植被仿真，建立植被与地形之间的定量描述关系；

校正输出模块，用于根据仿真关系模块所得结果，结合距离多普勒几何校正方式和归一化辐射校正方式进行校正，

所述归一化辐射校正方式实现如下，

cos(ψ)＝sin(θ₀)·cos(u)+cos(θ₀)·sin(u)·sin(v)＝A_area ^-1

其中，ψ为投影角，u和v分别表示传感器方向的地形坡度和方位角，A_area表示散射面积。

本发明还提供基于上述后向散射优化的森林复杂地形校正系统实现的树高反演系统，设置树高反演模块，用于针对不同树种之间的结构和分布的差异，提出树种多样性因子，对不同的树种的植被高度分别进行估计然后获得整个森林场景的树高反演结果，

树种多样性参数如下式所示，

其中，S表述树的种类数，N表示每一类树的数目，α表示多样性参数，则树高估计的表达式为，

其中，为体散射估计值，k_z表示有效垂直波数，φ₀表示干涉相位角。

本发明提出了一种基于后向散射优化的干涉SAR地形校正技术方案，并将其用于森林高度等垂直结构参数估计。针对森林垂直结构参数估计时存在景观破碎地形起伏多变的问题，首先建立了极化干涉SAR的森林相干模型，在对森林场景典型地物的主导散射机制分析的基础上，采用迭代的后向散射优化方法获取最佳后向散射估计值，然后采用基于归一化的辐射校正方法对地形的影响进行补偿，用于消除地形的影响；此外在树高等森林垂直结构参数的估计中，考虑不同的树种的分布以及结构的多样性，在垂直参数提取之时，首先对树种的分布进行建模和分类，用于区分入射角变化对不同地物的辐射影响，然后在分类的基础上采用基于模型的方法对垂直结构参数进行估计，从而提高了树高参数的估计精度。通过对真实的极化干涉SAR数据的实验并与实地测量结果进行比较发现，地形效应的影响得到了显著地改善，同时采用本文提出的方法可以显著的提高森林垂直结构参数的估计精度。

附图说明

图1是本发明实施例的流程图；

图2是本发明实施例中SAR仿真的流程图；

图3是本发明实施例中地形校正的流程图；

图4是本发明实施例方法和传统方法的地形校正后HH极化通道后向散射系数对比图。

图5是本发明实施例方法和传统方法的地形校正后HV极化通道后向散射系数对比图。

图6是本发明实施例方法和传统方法的地形校正后VV极化通道后向散射系数对比图。

具体实施方式

以下结合附图和实施例详细说明本发明技术方案。

本发明考虑结合仿真数据，对森林的极化散射特性进行定性分析，并集合几何与辐射的方法实现地形校正。如图1，针对森林复杂场景的地形效应分析，本发明首先分析了PolInSAR数据的森林散射模型的构造方法，特别的针对其中与地形相关的极化方位角、后向散射系数、局部入射角等变量。然后，重点分析了这些参数的统计纹理特征、极化散射特征、上下文特征的描述方法，并结合仿真数据分析，形成相应的参数校正算法，结合基于距离多普勒变换的几何校正和辐射归一化的方法实现地形校正。最后，结合样地数据获得森林复杂地表的散射模型，用于垂直结构参数反演。

本发明实施例可采用计算机软件技术实现自动流程运行，以下分步骤详细说明本发明实施例流程：

步骤1，针对森林场景的极化干涉SAR图像，根据森林场景的两层散射特点(地表层的表面散射和树冠层的体散射)，构建两层相干模型，建立相干系数、后向散射系数、入射角之间的定量关系。

首先根据树冠层的体散射与地表层的表面散射特性，构造两层结构的相干模型，即建立干涉相干系数的表达式：

其中，ω表示某一给定的极化状态，γ(ω)为相应的干涉相干系数，表示极化方位角，i是虚数单位；

m(ω)为地体辐射比，定义为

其中，m_G(ω)和m_V(ω)分别表示地表层的表面散射和树冠层的体散射的幅度，γ_V为纯体散射，表达式为：

其中，σ表示后向散射系数，即消光系数；θ表示入射角，可采用平均入射角，即图像中某个变化不大的区域的入射角的均值；k_z表示有效垂直波数，z表示垂直方向高度，h_V表示植被高度。该表达式建立了相干系数γ(ω)和入射角θ、后向散射系数σ之间的定量描述关系。

步骤2，在步骤1的相干模型的基础上，为了实现相干系数的精度估计，需要对地形引起的参数变化进行校正，具体实现步骤如下：

1)散射机制分解

在森林的两层散射结构中，可能存在的主要的散射形式包括：森林树冠体散射、地表后向散射、树干的直接后向散射、树冠和地表之间的后向散射、树干和地表之间的散射等形式。为了准确描述复杂地形的散射特性，采用了散射机制分类的方法进行地形效应分析。对于体散射分量，具体实施时可以通过选取两组合适的投影向量进行分离。对于表面散射而言，根据植被单次散射和二面角散射的极化相位在极化干涉协方差矩阵中符号不同而干涉相位符号相同的原理，采用Freeman-Durden的方法进行散射机制分类。

在干涉复相干矩阵中，体散射单元T_v可以表示为：

其中，m_v表示每个散射体中体散射的幅度，η表示粒子的平均形状。图像中每个点对应的都有一个入射角，可记为局部入射角θ₀。因此结合极化相干矩阵，选取合适的投影向量w₁和w₂，满足下式的条件，则可以实现体散射机制分离。

其中，系数(其中θ₀表示局部入射角，e为自然常数，z'表示地面参考高度，j是虚数单位)，对于体散射的提取，可以表示成：

其中，是投影向量w₁中的元素的共轭形式，w₂₁、w₂₂、w₂₃是投影向量w₂中的元素。化简后为：当提取出体散射单元后，对干涉复相干系数的优化问题可看成是非线性规划问题，可采用迭代优化过程实现，具体优化过程如下：

maximizieρ(w₁,w₂)

ρ(w₁,w₂)表示复相干系数。

对于表面散射机制的分离，本发明基于相干矩阵分解理论，采用Freeman-Durden分解的方法进行表面散射提取。在Freeman-Durden分解理论中，极化干涸矩阵有体散射、二次散射和表面散射组成，对于表面散射而言，在假设一阶Bragg模型的基础上，地面的表面散射可以表示为：

其中，T₁₁(i,j)表示相干矩阵T₁₁中第i行，第j列的元素，η表示粒子的平均形状。其中，i＝1,2,3，j＝1,2,3。

2)后向散射系数估计：特别的，针对地形引起的后向散射相位的变化，研究过程中将其转换成拉格朗日优化的形式，以便通过求解矩阵特征值得到地形校正的结果。

为了使得后向散射系数的估计结果满足最小方差准则，本发明实施例采用了基于迭代的后向散射系数迭代方法。假设极化干涉数据的先验信息为先验测量信息，以及对应的极化协方差矩阵为P是对应的是极化通道个数，例如：通道HH、HV、VV等。P可以取值为4,3,2，因此协方差矩阵就对应有C2矩阵、C3矩阵和C4矩阵。

则其任意的由发射到接收的后向散射可表示为(w表示后向散射相应的归一化的权值向量(w'w＝1)，

其中，表示w的伪逆。

对于全极化数据而言，w可以通过电磁波极化方向以及椭球模型的曲率角计算得到，其中ψ和χ分别表示发射接收的极化方向角和椭圆曲率角：

其中，ψ_r,χ_r,ψ_t,χ_t表示电磁波的发射和接收极化方式的极化变量(就是前面提到的方位角和椭圆曲率角)，w(ψ_r,χ_r,ψ_t,χ_t)即相应的w。例如：H方式发射，V方式接收，对应的就是HV通道。下标r和t分别表示发射和接收极化方式(t:tranmsmit,r:receive)。

那么对于一阶和二阶的极化向量p_t/r,1/2而言，满足如下的关系：

其中，为极化向量。

因而在这种情况下，可看成是散射类型、各向同性和各向异性、奇次或者偶次散射、水平极化和垂直极化之间的组合。因此，散射机制的决定了分辨率单元内的极化敏感性，而采用归一化变化矩阵得到协方差矩阵的估计结果C'_i：

C'_i＝UC_iU⁺

其中，'+'表示伪逆矩阵，此时，后向散射系数的估计可以看成如下的过程：a)最大化预测后向散射系数σ和先验信息之间的相关系数R(w)；b)最大化判定相关系数；c)最小化后向散射估计的均方误差。假设后向散射系数和观测数据满足线性关系，则先验信息和后向散射系数之间的相关系数R(w)可以表示为：

和分别为σ_i、y_i对应的均值，那么最佳极化状态权系数可表示为：

其中，其中和A是Hermitian正半定矩阵，B_i是Hermitian矩阵，n表示所有的观测数据个数。

为了获取后向散射系数的估计结果，本发明基于协方差矩阵C_i的连续，选取最佳的极化方式，从而获取最优的极化权系数，优化准则为：

在这种情况下，最佳权向量可表示为：

假设F为正半定Hermitian矩阵，上式的最优解可以转化成拉格朗日多项式的形式：

λ为拉格朗日因子，通过最大化上式的分子并保持分母不变可获得最佳估计结果。如果该计划状态和权向量负相关，则将导致最优相关系数为复数，则此时需要对矩阵A乘以-1然后再求解最佳极化状态量。对于上式的L以及对应的该优化问题可转化为矩阵特征值求解问题(F^-1表示矩阵求逆)：

F^-1Aw＝λw

上式的最初结果最为初始化近似值，为了进一步的提高估计精度，通过协方差矩阵变化方法，估计最佳极化状态的权向量，经过多次迭代后获取后向散射的估计值。具体实施时，可根据F^-1Aw＝λw得到w的估计结果，并判断是不是该条件下的最优解，若不是则返回继续计算，直到判断为是。

步骤3，在步骤2散射机制分类的基础上，结合仿真数据，模拟不同地形，不同植被覆盖条件下局部入射角以及后向散射系数的变化规律；建立植被覆盖(植被高度h_V)与地形参数(入射角θ)之间的定量描述关系，该关系是结合具体的数据得到的模拟结果，主要是根据步骤1所得纯体散射γ_V表达式分析植被覆盖(植被高度h_V)与地形参数(入射角θ)和纯体散射γ_V的关系。

参见图2，为了定性的描述不同地形和不同制备覆盖下的地物特征，具体实施时可以预先进行基于SAR模拟数据的仿真，通过对不同地形条件(例如坡度大小、破碎情况等)和不同植被覆盖(结构差异、树种差异、分布差异，例如稀疏林和茂密林、纯净林和混交林、灌木和乔木)的输入SAR数据模拟，基于地物特征表达定性的分析他们在统计纹理特征，极化散射特征和上下文特征上的差异，为进一步的建立森林复杂场景补偿模型提供数据支持。具体实施时，本领域技术人员可以设定复杂地形模拟，例如地形倾斜、随机起伏、景观破碎等，植被差异模拟可以考虑结构差异、树种差异和分布差异等。

具体实施时，本领域技术人员可预先进行仿真获取数据。仿真过程可分为两类，一类是对地形的仿真，一类是对植被的仿真。就地形而言，根据地形倾斜、随机起伏和破碎的情况分别进行仿真分析。例如，对于地形倾斜的情况，只需要给定一个倾斜坡度参数，就可以描述地形变化和后向散射之间的关系；对于随机起伏的情况，使得坡度变化比较剧烈，需要给定局部入射角，用于刻画其和后向散射的关系；而对于地形破碎情况，也主要体现在局部入射角的变化上，只是入射角属于不连续的变化。而对于植被的模拟，主要针对纯净林和混交林进行仿真分析，对于纯净林，分别对针叶林和和阔叶林进行仿真分析，分析不同植被覆盖下的后向散射特性，然后在针对混交林(针叶林和阔叶林混合的情况)，给定植被多样性参数，分析其后向散射的变化情况。最后将地形仿真和植被仿真相结合，分析坡度、入射角、不同植被覆盖及植被多样性参数和森林后向散射之间的关系。

步骤4，结合仿真实验建立的森林地形和后向散射系数以及局部入射角等参数的定性描述的基础上，本发明结合几何与辐射的方法对地形的影响进行校正。

参见图3，实施例的校正实现如下：

1)几何校正：一般包括坐标变换、轨道校正、RD变换

对图像像素进行精确地地理校正需要考虑SAR数据的辐射地形校正。对于地形编码来讲，需要高分辨率的DEM和传感器轨道平面的附加信息。通常来讲，图像像素的编码分类两类：前向编码和后向编码。在前向编码中，地面的每一个图像像素采用距离多普勒(RD)公式单独计算，后向散射则相反。参考现有技术，主要的处理过程包括：1)将坐标变化到一般的参考框架下；2)轨道一体化处理；3)通过距离多普勒迭代过程寻找最佳的图像像素和DEM对；4)在DEM上绘制图像像素图，并计算局部图像几何；5)斜距地面重构。为了精确地估计地理编码的精度，通过计算每一个GCPs(ground Control Points)位置的残差去表征。在成像几何(E,N)和图像几何(r,a)中有很多方法去估计散射面积。E和N分别表示东经和北纬的经纬度，r和a表示斜距范围的距离和方位向坐标。公式中的辐射归一化可以在成像几何和图像几何中表示，可采用地理编码的方法实现斜距-地矩的投影。

2)辐射校正：包括斜距投影、散射面积估计、后向散射估计

当SAR图像的局部图像几何进行描述和重构之后，这些信息可用于精确地的辐射计校正：

σ⁰＝β⁰/A_area

其中，A_area表示散射面积，σ⁰为图像空间的后向散射系数(地距上的)，β⁰为成像空间的后向散射系数(斜距上的)。可以看出，散射面积的精确估计不仅依赖于局部地形信息，同时也依赖DEM的分辨率。通过散射面积的空间依赖性的不同，辐射归一化方法分为成像空间估计法和图像空间估计法。通过在成像空间估计散射面积的方法可以评价辐射归一化的性能，但需要注意的是在给定的成像元素中所有的图像像素都需要选定，用于保留图像数据中能量。

由于地形扭曲带来的散射面积A的变化可通过投影角ψ进行参数化表示，投影角ψ为地面法向量和图像平面余角的最小值，可以通过SAR系统的观测向量和局部地形的角度和方位进行推导：

cos(ψ)＝sin(θ₀)·cos(u)+cos(θ₀)·sin(u)·sin(v)＝A_area ^-1

其中，θ₀表示局部入射角，u和v分别表示传感器方向的地形坡度和方位角。

经几何与辐射校正后，利用DEM提供的高程和地理位置信息得到最终的地形校正结果。

参见图4、5、6，相应于观测信息，对比本发明基于后向散射优化的方法和传统的距离多普勒方法、辐射归一化方法所得后向散射系数，明显本发明技术方案的效果更好。

在步骤4地形校正的基础上，本发明在进一步的树高反演研究中，考虑到不同树种之间的结构和分布的差异，提出了树种多样性因子，对不同的树种的植被高度分别进行估计然后获得整个森林场景的树高反演结果。

实施例提供的引入树种多样性因子的三阶段反演算法实现如下：

1)最小平方直线拟合

根据不同极化通道复相干系数，采用最小方差准则，在干涉单位圆上获取最佳直线匹配结果，若仅采用两个极化通道，则计算通过单位圆上两点的直线。

2)地形相位估计

由于最佳地形相位值位于最佳直线和单位圆的交点之间，因此分别计算具有最大体散射和最大表面散射的相对位置进行估计，并采用地形校正的后向散射面积估计最佳地形相位。

3)植被高度估计

由于不同树种的分布和结构上的差异，因此，为了植被提高估计精度，本发明引入Fisher Alpha-Diversity方法，计算树种多样性参数α：

其中，S表示类别数，n表示植被总数，表示多样性参数α。则此时植被高度估计h_v如下所示：

其中，为体散射估计值，k_z表示有效垂直波数，φ₀表示干涉相位角。为验证本发明技术效果，可以根据本发明提出的复杂地表地形校正算法对校正前后的后向散射估计结果进行验证。此外，进一步的结合实地调研数据，比较地形校正前后的植被高度反演结果的精度，也可以对本发明技术进行验证。

具体实施时，也可以采用模块方式提供相应系统。本发明实施例相应提供一种后向散射优化的森林复杂地形校正系统，包括以下模块：

相干模型构建模块，用于针对森林场景的极化干涉SAR图像，根据森林场景的两层散射特点，构造两层结构的相干模型如下，

其中，ω表示某一给定的极化状态，γ(ω)为相应的干涉相干系数，表示极化方位角，i是虚数单位；

m(ω)为地体辐射比，定义为

其中，m_G(ω)和m_V(ω)分别表示地表层的表面散射和树冠层的体散射的幅度，γ_V为纯体散射，表达式如下，

其中，σ表示后向散射系数，θ表示入射角，k_z表示有效垂直波数，z表示垂直方向高度，h_V表示植被高度；

地形变化校正模块，用于在相干模型基础上，对地形引起的参数变化进行校正，实现如下，首先进行散射机制分解，

对于表面散射，采用Freeman-Durden方法进行散射机制分类；

对于体散射，选取满足下式的投影向量w₁和w₂，实现体散射机制分离，

其中，T_v为体散射单元，m_v表示每个散射体中体散射的幅度，η表示粒子的平均形状，C_v为系数，定义为

其中，θ₀表示局部入射角，e为自然常数，z'表示地面参考高度，j是虚数单位；

然后，针对地形引起的后向散射相位的变化，转换成拉格朗日优化的形式，通过求解矩阵特征值得到地形校正的结果，

F^-1Aw＝λw

其中，F、A为正半定Hermitian矩阵，λ为拉格朗日因子，w表示后向散射相应的归一化的权值向量；

仿真关系模块，用于在散射机制分解结果的基础上，结合地形仿真和植被仿真，建立植被与地形之间的定量描述关系；

校正输出模块，用于根据仿真关系模块所得结果，结合距离多普勒几何校正方式和归一化辐射校正方式进行校正，

所述归一化辐射校正方式实现如下，

cos(ψ)＝sin(θ₀)·cos(u)+cos(θ₀)·sin(u)·sin(v)＝A_area ^-1

其中，ψ为投影角，u和v分别表示传感器方向的地形坡度和方位角，A_area表示散射面积。

基于上述后向散射优化的森林复杂地形校正系统增加设置树高反演模块，即可实现树高反演系统。所述树高反演模块用于针对不同树种之间的结构和分布的差异，提出树种多样性因子，对不同的树种的植被高度分别进行估计然后获得整个森林场景的树高反演结果，

树种多样性参数如下式所示，

其中，S表述树的种类数，N表示每一类树的数目，α表示多样性参数，则树高估计的表达式为，

其中，为体散射估计值，k_z表示有效垂直波数，φ₀表示干涉相位角。

各模块具体实现可参见相应步骤，本发明不予赘述。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (4)

1.一种后向散射优化的森林复杂地形校正方法，其特征在于，包括以下步骤：

步骤1，针对森林场景的极化干涉SAR图像，根据森林场景的两层散射特点，构造两层结构的相干模型如下，

其中，ω表示某一给定的极化状态，γ(ω)为相应的干涉相干系数，表示极化方位角，i是虚数单位；

m(ω)为地体辐射比，定义为

其中，m_G(ω)和m_V(ω)分别表示地表层的表面散射和树冠层的体散射的幅度，

γ_V为纯的体散射，表达式如下，

${\mathrm{\γ}}_V=\frac{\underset{0}{\overset{h_v}{\∫}}\exp (2\mathrm{\σ}z/cos\mathrm{\θ})\exp \left({\mathrm{ik}}_{z}z\right)dz}{\underset{0}{\overset{h_v}{\∫}}\exp (2\mathrm{\σ}z/cos\mathrm{\θ})dz}$

其中，σ表示后向散射系数，θ表示入射角，k_z表示有效垂直波数，z表示垂直方向高度，h_v表示植被高度；

步骤2，在步骤1所得相干模型基础上，对地形引起的参数变化进行校正，实现如下，

首先进行散射机制分解，

对于表面散射，采用Freeman-Durden方法进行散射机制分类；

对于体散射，选取满足下式的投影向量w₁和w₂，实现体散射机制分离，

$T_v=m_v{C}_v{w}_1^H\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{\η}& 0\\ 0& 0& \mathrm{\η}\end{array}\right]w_2=0$

其中，T_v为体散射单元，m_v表示每个散射体中体散射的幅度，η表示粒子的平均形状，C_v为系数，定义为

其中，θ₀表示局部入射角，e为自然常数，z'表示地面参考高度，j是虚数单位；

然后，针对地形引起的后向散射相位的变化，转换成拉格朗日优化的形式，通过求解矩阵特征值得到地形校正的结果，

F^-1Aw＝λw

其中，F、A为正半定Hermitian矩阵，λ为拉格朗日因子，w表示后向散射相应的归一化的权值向量；

步骤3，在步骤2所得散射机制分解结果的基础上，结合地形仿真和植被仿真，建立植被与地形之间的定量描述关系；

步骤4，根据步骤3所得结果，结合距离多普勒几何校正方式和归一化辐射校正方式进行校正，

所述归一化辐射校正方式实现如下，

cos(ψ)＝sin(θ₀)·cos(u)+cos(θ₀)·sin(u)·sin(v)＝A_area ^-1

其中，ψ为投影角，u和v分别表示传感器方向的地形坡度和方位角，A_area表示散射面积。

2.根据权利要求1所述后向散射优化的森林复杂地形校正方法实现的树高反演方法，其特征在于：针对不同树种之间的结构和分布的差异，提出树种多样性因子，对不同的树种的植被高度分别进行估计然后获得整个森林场景的树高反演结果，

树种多样性参数如下式所示，

$S=\mathrm{\α}ln(1+\frac{N}{\mathrm{\α}})$

其中，S表述树的种类数，N表示每一类树的数目，α表示多样性参数，则树高估计的表达式为，

$h_v=\frac{1}{k_z}\{\arg \left({\widetilde{\mathrm{\γ}}}_V\exp \right(-{\mathrm{i\φ}}_0\left)\right)+\mathrm{\α}(\mathrm{\π}-2{\sin }^{-1}(|{\widetilde{\mathrm{\γ}}}_V{|}^{0.8}\left)\right)\}$

其中，为体散射估计值，k_z表示有效垂直波数，φ₀表示干涉相位角。

3.一种后向散射优化的森林复杂地形校正系统，其特征在于：包括以下模块，

相干模型构建模块，用于针对森林场景的极化干涉SAR图像，根据森林场景的两层散射特点，构造两层结构的相干模型如下，

其中，ω表示某一给定的极化状态，γ(ω)为相应的干涉相干系数，表示极化方位角，i是虚数单位；

m(ω)为地体辐射比，定义为

其中，m_G(ω)和m_V(ω)分别表示地表层的表面散射和树冠层的体散射的幅度，

γ_V为纯的体散射，表达式如下，

${\mathrm{\γ}}_V=\frac{\underset{0}{\overset{h_v}{\∫}}\exp (2\mathrm{\σ}z/cos\mathrm{\θ})\exp \left({\mathrm{ik}}_{z}z\right)dz}{\underset{0}{\overset{h_v}{\∫}}\exp (2\mathrm{\σ}z/cos\mathrm{\θ})dz}$

其中，σ表示后向散射系数，θ表示入射角，k_z表示有效垂直波数，z表示垂直方向高度，h_v表示植被高度；

地形变化校正模块，用于在相干模型基础上，对地形引起的参数变化进行校正，实现如下，

首先进行散射机制分解，

对于表面散射，采用Freeman-Durden方法进行散射机制分类；

对于体散射，选取满足下式的投影向量w₁和w₂，实现体散射机制分离，

$T_v=m_v{C}_v{w}_1^H\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{\η}& 0\\ 0& 0& \mathrm{\η}\end{array}\right]w_2=0$

其中，T_v为体散射单元，m_v表示每个散射体中体散射的幅度，η表示粒子的平均形状，C_v为系数，定义为

其中，θ₀表示局部入射角，e为自然常数，z'表示地面参考高度，j是虚数单位；

然后，针对地形引起的后向散射相位的变化，转换成拉格朗日优化的形式，通过求解矩阵特征值得到地形校正的结果，

F^-1Aw＝λw

其中，F、A为正半定Hermitian矩阵，λ为拉格朗日因子，w表示后向散射相应的归一化的权值向量；

仿真关系模块，用于在散射机制分解结果的基础上，结合地形仿真和植被仿真，建立植被与地形之间的定量描述关系；

校正输出模块，用于根据仿真关系模块所得结果，结合距离多普勒几何校正方式和归一化辐射校正方式进行校正，

所述归一化辐射校正方式实现如下，

cos(ψ)＝sin(θ₀)·cos(u)+cos(θ₀)·sin(u)·sin(v)＝A_area ^-1

其中，ψ为投影角，u和v分别表示传感器方向的地形坡度和方位角，A_area表示散射面积。

4.根据权利要求3所述后向散射优化的森林复杂地形校正系统实现的树高反演系统，其特征在于：设置树高反演模块，用于针对不同树种之间的结构和分布的差异，提出树种多样性因子，对不同的树种的植被高度分别进行估计然后获得整个森林场景的树高反演结果，

树种多样性参数如下式所示，

$S=\mathrm{\α}ln(1+\frac{N}{\mathrm{\α}})$

其中，S表述树的种类数，N表示每一类树的数目，α表示多样性参数，则树高估计的表达式为，

$h_v=\frac{1}{k_z}\{\arg \left({\widetilde{\mathrm{\γ}}}_V\exp \right(-{\mathrm{i\φ}}_0\left)\right)+\mathrm{\α}(\mathrm{\π}-2{\sin }^{-1}(|{\widetilde{\mathrm{\γ}}}_V{|}^{0.8}\left)\right)\}$

其中，为体散射估计值，k_z表示有效垂直波数，φ₀表示干涉相位角。