CN112363161B - 基于散射机制分解的植被垂直结构及林下地形反演方法及装置 - Google Patents

Abstract

本发明公开了一种基于散射机制分解的植被垂直结构及林下地形反演方法及装置，其方法为：获取森林反演区域多基线多极化的SAR数据集，并进行配准、去平地处理，得到极化层析SAR数据；对极化层析SAR数据计算其多基线多极化协方差矩阵，并分解计算冠层散射和地表散射的干涉协方差矩阵；针对冠层散射和地表散射的干涉协方差矩阵，均进行矢量化和稀疏表达，得到对应的稀疏贝叶斯学习模型；采用EM算法迭代更新稀疏贝叶斯学习模型中的超参数，进而获得冠层散射和地表散射的层析谱；根据地表散射和冠层散射的层析谱，分别提取森林反演区域的林下地形和植被高度。本发明可以利用少量、非均匀分布的SAR数据高精度地反演植被高和林下地形。

Description

基于散射机制分解的植被垂直结构及林下地形反演方法及 装置

技术领域

本发明涉及雷达遥感对地观测技术，具体涉及一种在少量非均匀采样数据条件下基于散射机制分解的超分辨率极化SAR层析植被垂直结构及林下地形反演方法及装置。

背景技术

生态文明建设是我国发展的重要布局之一，森林资源的调查与保护在其中扮演着重要的角色。森林中包含着种类繁多的植被，是地球上最宝贵的自然资源之一。森林在维持全球的碳氧平衡、维护生态***的多样性、调节大气湿度等方面都有着巨大的作用。森林中植被的垂直结构是量化分析森林的重要参数，对全球生物量的估计、林下地形测量、森林动态监测等有着重要的意义。

长波段合成孔径雷达技术(Synthetic Aperture Radar,SAR)具有较强的穿透力，在森林的透视测绘中有着重要的应用。合成孔径雷达层析技术(Synthetic ApertureRadar Tomography,TomoSAR)通过传感器多次近似平行的飞行，在高度向和方位向上分别形成合成孔径，具备三维成像能力(图1)。在雷达信号中，森林植被与林下地形有着不同的散射机制，而SAR影像不同极化通道对不同散射机制的敏感性不同。因此极化层析SAR技术(Polarimetric Synthetic Aperture Radar Tomography,Pol-TomoSAR)正成为森林植被参数反演的重要技术手段。

在植被的SAR层析领域，周期图法、自适应波束形成法、多重信号分类法等已广泛用于植被的三维成像。但这类常规谱分析方法的分辨率和重构精度往往受到采样数量和采样分布的限制。观测数量少、基线分布不均匀往往会使得层析结果出现严重的旁瓣效应、分辨率低、受噪声影响大。而增加观测数据则意味着大量的观测成本，同时带来时间失相干的影响。压缩感知方法可以很好地克服数据采样的限制，但传统基于l1范数凸优化的压缩感知方法(Compressive Sensing,CS)需要设置用户参数，很难选择一个普遍适用于所有模型和场景的参数。此外，凸优化运算量巨大，限制了该方法的应用。

发明内容

基于上述现有技术存在的技术问题，本发明提供一种基于散射机制分解的植被垂直结构及林下地形反演方法及装置，将森林反演区域冠层和地表两种不同散射机制的信号进行稀疏表达，进而将其纳入协方差矢量的稀疏贝叶斯学习模型中以求解对应的层析谱，从而实现少量的观测数据重构高精度的植被垂直结构和林下地形。本发明具有更接近于l0范数的全局最优解，无需设置用户参数，具有更广泛的实用价值。

为实现上述技术目的，本发明采用如下技术方案：

一种基于散射机制分解的植被垂直结构及林下地形反演方法，包括以下步骤：

步骤1，获取森林反演区域多基线多极化的SAR数据集，并进行配准、去平地处理，得到极化层析SAR数据；

步骤2，对步骤1得到的极化层析SAR数据计算其多基线多极化协方差矩阵，并分解计算冠层散射和地表散射的干涉协方差矩阵；

步骤3，针对冠层散射和地表散射的干涉协方差矩阵，均进行矢量化和稀疏表达，分别得到冠层散射和地表散射的稀疏贝叶斯学习模型；

步骤4，采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数，进而获得冠层散射和地表散射的层析谱；

步骤5，根据地表散射的层析谱提取森林反演区域的林下地形，根据冠层散射的层析谱提取森林反演区域的植被高度。

在更优的技术方案中，步骤2中多基线多极化协方差矩阵的计算方法为：

式中，y为步骤1得到的极化层析SAR数据，

g_HV＝g_VH，g_HH、g_HV、g_VH、g_VV分别为HH、HV、VH、VV极化方式对应的SAR数据构成的N维SAR数据观测矩阵，N为基线的数量；(*)^H为共轭转置操作，E[·]为求期望运算；W为计算得到的多基线多极化协方差矩阵；C_g和R_g分别为地表散射的极化协方差矩阵和干涉协方差矩阵，C_v和R_v分别为冠层散射的极化协方差矩阵和干涉协方差矩阵，

为克罗内克积；

对多基线多极化协方差矩阵W进行SVD分解，使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g。

在更优的技术方案中，使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g的方法为：

首先，取SVD分解得到的前两个最大特征值对应的特征向量重塑矩阵得到R₁和R₂，使用R₁和R₂表示冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g为：

式中，a、b为一对实数；

然后，基于冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g的半正定约束，即R_v和R_g的特征值非负，求解实数a、b的取值区间；

再后，根据以下散射机制最大化分离的准则，求解实数a、b：

其中，trace()为求矩阵的迹，||·||_F为矩阵的Frobenius范数。最后，根据求解得到的实数a、b和矩阵R₁和R₂，即可得到冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g。

在更优的技术方案中，由步骤3得到地表散射和冠层散射的稀疏贝叶斯学习模型分别为：

式中，R_g、R_v分别为地表散射和冠层散射的干涉协方差矩阵，vec(·)为矢量化操作，y_g、y_v分别为地表散射和冠层散射的干涉协方差矩阵的矢量表示；

B为中间参数，且B＝A⊙conj(A)，⊙为Khatri-Rao积，A为中间参数，A＝exp(j2πξ_ns_d)，s_d为散射体的空间高度，ξ_n为空间频率；Ψ_g表示地表散射的小波稀疏基，采用单位正交基，Ψ_v表示冠层散射的小波稀疏基；

p_g、p_v分别为地表散射和冠层散射的后向散射功率矢量；

σ_g、σ_v分别为地表散射和冠层散射对应的噪声功率，t_g＝vec(σ_gI)、t_v＝vec(σ_vI)分别为对应噪声功率矩阵的矢量化操作，I为单位矩阵；δ_g和δ_v分别为地表散射和冠层散射的协方差矩阵与对应真实协方差矩阵之间的估计误差。

在更优的技术方案中，步骤4中采用EM算法迭代更新稀疏贝叶斯学习模型中的超参数并获得冠层散射和地表散射的层析谱，具体为：

步骤a，定义先验分布：

为了便于阐述，将步骤3得到的地表散射和冠层散射的稀疏贝叶斯学习模型在形式上统一表示为：Y＝Gx+δ；其中，G＝[BΨ^H,I]表示在稀疏基下新的测量矩阵，x＝[Ψp,t]为待求稀疏域下的后向散射功率和噪声矢量，Y为地表散射或冠层散射的干涉协方差矩阵的矢量表示；

根据稀疏操作后的散射信号是稀疏这一先验信息，定义其服从高斯分布，即x～N(0,Λ)；其中Λ＝diag(Γ)，

为稀疏贝叶斯学习模型的超参数矢量，表征了x的空间分布；前D个超参数τ₁,...τ_D控制高度向上D个散射体后向散射功率的稀疏性，后N²个超参数

为***噪声；则Y关于Γ的概率密度函数为：

p(Y|Γ)＝∫p(Y|x)p(x|Γ)dx；

步骤b，参数初始值定义：

设置稀疏贝叶斯学习模型的超参数向量Γ初始值和迭代终止条件；

步骤c，计算后验概率分布：

基于贝叶斯准则和步骤a定义的先验分布，计算x的均值矢量和协方差矩阵：

Σ_Y＝V+GΛ⁽ⁱ⁾G^H

式中，上标i表示当前迭代次数，μ为x的均值矢量，Σ_x和Σ_Y分别为x和Y的协方差矩阵，V为样本协方差矩阵估计误差δ的二阶统计。

步骤d，基于步骤c得到的均值矢量和协方差矩阵，采用EM算法学习更新超参数向量：

步骤e，重复步骤c、d，直至满足迭代终止条件，使用超参数向量Γ中的前D个超参数Γ_(1…D)计算冠层散射和地表散射的层析谱：

在更优的技术方案中，迭代终止条件为||Γⁱ⁺¹-Γⁱ||₂/||Γⁱ||₂≤tol，和/或达到最大迭代次数。

在更优的技术方案中，根据地表散射的层析谱提取森林反演区域的林下地形的方法为：提取地表散射的层析谱中功率最大值所在的位置作为地面高度的估计值。

在更优的技术方案中，根据冠层散射的层析谱提取森林反演区域的植被高度的方法为：通过功率衰减法确定冠层顶部的位置，即：以冠层相位中心处为起始点，以预设的步长向上衰减，当功率衰减处提取的高度值和Lidar数据提取的冠层高度(DSM)的均方根误差达到最小时，将此时的高度值作为冠层高度的估计值；通过冠层高度减去地面高度即为植被高度。

本发明还提供一种基于散射机制分解的植被垂直结构及林下地形反演装置，包括：

数据预处理模块，用于：获取森林反演区域多基线多极化的SAR数据集，并进行配准、去平地处理，得到极化层析SAR数据；

协方差计算分解模块，用于：对数据预处理模块得到的极化层析SAR数据计算其多基线多极化协方差矩阵，并分解计算冠层散射和地表散射的干涉协方差矩阵；

模型构建模块，用于：针对冠层散射和地表散射的干涉协方差矩阵，均进行矢量化和稀疏表达，分别得到冠层散射和地表散射的稀疏贝叶斯学习模型；

层析谱计算模块，用于：采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数，进而获得冠层散射和地表散射的层析谱；

植被垂直结构及林下地形反演模块，用于：根据地表散射的层析谱提取森林反演区域的林下地形，根据冠层散射的层析谱提取森林反演区域的植被高度。

有益效果

本发明使用多基线全极化的SAR数据，对冠层和地表两种不同散射机制进行了分离与层析。一方面，可以获得更加准确的相位中心估计值，提高林下地形与植被高度的反演精度；另一方面，使用不同稀疏基对不同散射机制进行稀疏表达，使得本发明方法能在森林植被区域合理地运用。

而且，本发明方法是一种近似于l0范数的稀疏重构方法，相较于传统基于l1范数凸优化的压缩感知方法，可以获得更稀疏的全局最优解，能够以较少的测量数据满足RIP条件(即感知矩阵G列不相关)；无需艰难地设置噪声上限参数，可以自适应噪声的影响，具备更高的植被高和林下地形重构精度；同时，本发明的方法具备更高的运算效率，适用于大范围、大尺度的林下地形测绘以及植被高度反演,具有更高的生产实用价值。

另外，本发明顾及了层析SAR数据样本协方差矩阵的估计误差，该方法更加稳健。

附图说明

图1是层析SAR森林三维成像示意图；

图2是本发明实施例所述方法的流程图；

图3是本发明实施例所述方法各步骤子流程图。

具体实施方式

为了更加清楚地说明本发明的目的、实施流程和技术优势，本实施例选取Tropi-SAR 2009项目ONERA SETHI飞行***获得的P波段全极化的SAR数据，对本发明所提方法进行进一步的详细说明。应当理解，此处所描述的具体实施仅用于解释本发明，并不用于限定本发明。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

本发明基于散射机制分解的植被垂直结构及林下地形反演方法，为超分辨率的稀疏重构方法。为了说明方法的优越性能，具体实施中分析了不同层析SAR方法在不同数目的全极化SAR数据下的性能，然后仅采用3条基线的全极化SAR数据反演林下地形和植被高度。

本实施例的具体实施步骤参照图2、图3所示，包括：

步骤1，获取森林反演区域多基线多极化的SAR数据集，并进行配准、去平地处理，得到极化层析SAR数据。具体为：

(1)选取TropiSAR2009项目全极化SAR影像作为实验数据。

(2)配准：对于指定的极化通道(例如：HH极化通道)，在获取的N景SAR影像中，选择其中一景影像为参考影像，其余N-1景影像应用最大化相干系数配准到参考影像的坐标框架中。然后将剩余极化通道(HV、VH、VV)的SAR影像采用与HH极化通道相同的多项式进行配准。

(3)去平地：即去除由平地效应引起的干涉相位。

(4)数据离散化：经上述处理后，对于选定极化方式的极化SAR层析数据，距离-方位向平面像素(r,x)聚焦后的复数值为：

式中b_n为垂直基线，Δs为高度向的采样范围，γ(r,x,s)为高度向上的反射率函数，ξ_n＝2πb_n/λr为空间频率。

HH、HV、VH、VV极化方式对应的SAR数据构成的N维SAR数据观测值分别为g_HH、g_HV、g_VH、g_VV；N为基线的数量；在单站后向散射体制下g_HV＝g_VH，因此可将步骤1最终得到的极化层析SAR数据表示为

对于多基线的层析SAR数据集，对高度向上连续的信号进行D次离散的采样后，可将(1)表达为以下矩阵形式：

g＝Aγ+e (2)

其中g是N维观测值。A＝exp(j2πξ_ns_d)，s_d为散射体的空间高度，e是N维的噪声向量。

步骤2，对极化层析SAR数据计算其多基线多极化协方差矩阵，并分解计算冠层散射和地表散射的干涉协方差矩阵。

对于森林植被的SAR层析而言，往往利用不同极化通道对不同散射机制的敏感性不同来提取林下地形和植被高度，例如在HH极化通道的层析谱中提取林下地形，在HV极化通道的层析谱中提取冠层高度。但由于森林场景较为复杂，各极化后向散射能量并不完全具备如此规律的分布，前述提取不同散射机制的策略可能会造成一定的偏差。此外，本发明基于散射机制分解的极化SAR层析技术要求所测量的信号是稀疏的(信号矢量中只有少量元素不为零)。在森林区域，植被信号的散射是一个连续的过程，而地面信号本身具备稀疏性，因此本发明对冠层和地表两种不同散射机制的信号采用不同的稀疏基进行稀疏表达。

为了更加准确地提取出不同散射机制相位中心的位置，以及对不同的散射机制进行稀疏表达，本发明对步骤1得到的极化层析SAR数据计算其协方差矩阵，即得到多基线多极化协方差矩阵，进而使用代数合成的方法对两种不同的散射机制(地表散射、冠层散射)进行分离：

(*)^H为共轭转置操作，E[·]为求期望运算；W为计算得到的多基线多极化协方差矩阵；C_g和R_g分别为地表散射的极化协方差矩阵和干涉协方差矩阵，C_v和R_v分别为冠层散射的极化协方差矩阵和干涉协方差矩阵，

为克罗内克积。

得到多基线多极化协方差矩阵W后，对其进行SVD分解，使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g：

首先，取SVD分解得到的前两个最大特征值对应的特征向量重塑矩阵得到R₁和R₂，使用矩阵R₁和R₂表示冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g为：

式中，a、b为一对实数；

然后，基于冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g的半正定约束，即R_v和R_g的特征值非负，求解实数a、b的取值区间；

再后，根据以下散射机制最大化分离的准则，求解实数a、b：

其中，trace()为求矩阵的迹，||·||_F为矩阵的Frobenius范数。最后，根据求解得到的实数a、b和矩阵R₁和R₂，即可得到冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g。

步骤3，针对冠层散射和地表散射的干涉协方差矩阵，均进行矢量化和稀疏表达，分别得到冠层散射和地表散射的稀疏贝叶斯学习模型。

在步骤2分解得到冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g后，分别进行矢量化和稀疏表达，得到地表散射和冠层散射的稀疏贝叶斯学习模型分别为：

式中，R_g、R_v分别为地表散射和冠层散射的干涉协方差矩阵，vec(·)为矢量化操作，y_g、y_v分别为地表散射和冠层散射的干涉协方差矩阵的矢量表示；

B为中间参数，且B＝A⊙conj(A)，⊙为Khatri-Rao积，A为中间参数，A＝exp(j2πξ_ns_d)，s_d为散射体的空间高度，ξ_n为空间频率；Ψ_g表示地表散射的小波稀疏基，采用单位正交基，Ψ_v表示冠层散射的小波稀疏基；

p_g、p_v分别为地表散射和冠层散射的后向散射功率矢量；

σ_g、σ_v分别为地表散射和冠层散射对应的噪声功率，t_g＝vec(σ_gI)、t_v＝vec(σ_vI)分别为对应噪声功率矩阵的矢量化操作，I为单位矩阵；δ_g和δ_v分别为地表散射和冠层散射的协方差矩阵与对应真实协方差矩阵之间的估计误差，其二阶统计量分别为V_g,V_v。

步骤4，采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数，进而获得冠层散射和地表散射的层析谱。具体为：

步骤a，定义先验分布：

为了便于阐述，将步骤3得到的地表散射和冠层散射的稀疏贝叶斯学习模型在形式上统一表示为：Y＝Gx+δ；其中，G＝[BΨ^H,I]表示在稀疏基下新的测量矩阵，x＝[Ψp,t]为待求稀疏域下的后向散射功率和噪声矢量，Y为地表散射或冠层散射的干涉协方差矩阵的矢量表示；

根据稀疏操作后的散射信号(即x)是稀疏这一先验信息，定义其服从高斯分布，即x～N(0,Λ)；其中Λ＝diag(Γ)，

为稀疏贝叶斯学习模型的超参数矢量，表征了x的空间分布；前D个超参数τ₁,...τ_D控制高度向上D个散射体后向散射功率的稀疏性，后N²个超参数

为***噪声；则Y关于Γ的概率密度函数为：

p(Y|Γ)＝∫p(Y|x)p(x|Γ)dx；

步骤b，参数初始值定义：

设置稀疏贝叶斯学习模型的超参数向量Γ初始值和迭代终止条件，如：(Γ)⁰是一个D+N²维的随机非负矢量，Λ⁽⁰⁾＝diag(Γ⁽⁰⁾)；迭代终止误差||Γⁱ⁺¹-Γⁱ||₂/||Γⁱ||₂≤tol＝10^-4，最大迭代次数为10次；

步骤c，计算后验概率分布：

基于贝叶斯准则和步骤a定义的先验分布，计算x的均值矢量和相关协方差矩阵：

Σ_Y＝V+GΛ⁽ⁱ⁾G^H

式中，上标i表示当前迭代次数，μ为x的均值矢量，Σ_x和Σ_Y分别为x和Y的协方差矩阵，V为样本协方差矩阵估计误差δ的二阶统计；

步骤d，基于步骤c得到的均值矢量和协方差矩阵，采用EM算法(即期望最大化算法)学习更新超参数向量：

超参数向量的更新同对以下目标函数进行最小化进行：

步骤e，重复步骤c、d，直至满足任一个迭代终止条件，使用超参数向量Γ中的前D个超参数Γ_(1…D)计算冠层散射和地表散射的层析谱(即空间功率谱)：

由于已对不同散射机制的信号进行了稀疏化操作，故用于表达后向散射功率的小波系数在稀疏域中大部分位置的值都是0值，故本发明中假设稀疏化后的森林观测信号x服从高斯分布符合先验信息，具有合理性。同时，高斯先验分布十分有利于稀疏贝叶斯学习的开展。

步骤5，根据地表散射的层析谱计算森林反演区域的林下地形，根据冠层散射的层析谱计算森林反演区域的植被高度。具体为：

林下地形的计算方法为：将地表散射的层析谱中功率最大值所在的位置作为地面高度的估计值。

植被高度的计算方法为：通过功率衰减法确定冠层顶部的位置，即：以冠层相位中心处为起始点，以预设的步长-0.1db向上衰减，当功率衰减处提取的高度值和Lidar数据DSM的均方根误差达到最小时，将此时的高度值作为冠层高度的估计值；通过冠层高度减去地面高度即为植被高度。

多基线多极化的SAR数据中的每个像元，均按上述方法遍历，即可获得整个森林反演区域的林下地形和植被高度。

通过对实例数据反演的不同散射层析谱、林下地形和植被高度进行分析，本发明方法能够在少量观测数据下获得清晰的层析谱，几乎不受噪声影像，分辨率高。而传统谱分析方法Beamforming和Capon即使使用了全极化的数据，仍然不能完全地区分不同散射机制的散射体。随着基线数量的下降，两种方法的性能受到了很大的影响，以致不能区分不同高度的散射体。而本发明方法受影像数量的影响较小，即使在仅使用3条基线SAR数据的情况下，仍然保持着较高的重构性能。此外，本发明方法仅在仅使用3条基线全极化SAR数据的情况下，也获得了可靠的林下地形和植被高反演结果。

本发明还提供一种基于散射机制分解的植被垂直结构及林下地形反演装置，与上述反演方法对应，包括：

数据预处理模块，用于：获取森林反演区域多基线多极化的SAR数据集，并进行配准、去平地处理，得到极化层析SAR数据；

协方差计算分解模块，用于：对数据预处理模块得到的极化层析SAR数据计算其多基线多极化协方差矩阵，并分解计算冠层散射和地表散射的干涉协方差矩阵；

模型构建模块，用于：针对冠层散射和地表散射的干涉协方差矩阵，均进行矢量化和稀疏表达，分别得到冠层散射和地表散射的稀疏贝叶斯学习模型；

层析谱计算模块，用于：采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数，进而获得冠层散射和地表散射的层析谱；

植被垂直结构及林下地形反演模块，用于：根据地表散射的层析谱提取森林反演区域的林下地形，根据冠层散射的层析谱提取森林反演区域的植被高度。

以上实施例为本申请的优选实施例，本领域的普通技术人员还可以在此基础上进行各种变换或改进，在不脱离本申请总的构思的前提下，这些变换或改进都应当属于本申请要求保护的范围之内。

Claims (7)

1.一种基于散射机制分解的植被垂直结构及林下地形反演方法，其特征在于，包括以下步骤：

步骤1，获取森林反演区域多基线多极化的SAR数据集，并进行配准、去平地，得到极化层析SAR数据；

步骤2，对步骤1得到的极化层析SAR数据计算其多基线多极化协方差矩阵，并分解计算冠层散射和地表散射的干涉协方差矩阵；

其中，多基线多极化协方差矩阵的计算方法与分解方法为：

式中，y为步骤1得到的极化层析SAR数据，

g_HV＝g_VH，g_HH、g_HV、g_VH、g_VV分别为HH、HV、VH、VV极化方式对应的SAR数据构成的N维SAR数据观测值；(*)^H为共轭转置操作，E[·]为求期望运算；W为计算得到的多基线多极化协方差矩阵；C_g和R_g分别为地表散射的极化协方差矩阵和干涉协方差矩阵，C_v和R_v分别为冠层散射的极化协方差矩阵和干涉协方差矩阵，

为克罗内克积；

对多基线多极化协方差矩阵W进行SVD分解，使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g：

首先，取SVD分解得到的前两个最大特征值对应的特征向量重塑矩阵得到R₁和R₂，使用R₁和R₂表示冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g为：

式中，a、b为一对实数；

然后，基于冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g的半正定约束，即R_v和R_g的特征值非负，求解实数a、b的取值区间；

再后，根据以下散射机制最大化分离的准则，求解实数a、b：

其中，trace()为求矩阵的迹，||·||_F为矩阵的Frobenius范数；

最后，根据求解得到的实数a、b和矩阵R₁和R₂，即可得到冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g；

步骤3，针对冠层散射和地表散射的干涉协方差矩阵，均进行矢量化和稀疏表达，分别得到冠层散射和地表散射的稀疏贝叶斯学习模型；

步骤4，采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数，进而获得冠层散射和地表散射的层析谱；

步骤5，根据地表散射的层析谱提取森林反演区域的林下地形，根据冠层散射的层析谱提取森林反演区域的植被高度。

2.根据权利要求1所述的方法，其特征在于，由步骤3得到地表散射和冠层散射的稀疏贝叶斯学习模型分别为：

式中，R_g、R_v分别为地表散射和冠层散射的干涉协方差矩阵，vec(·)为矢量化操作，y_g、y_v分别为地表散射和冠层散射的干涉协方差矩阵的矢量表示；

B为中间参数，且B＝A⊙conj(A)，⊙为Khatri-Rao积，A为中间参数，A＝exp(j2πξ_ns_d)，s_d为散射体的空间高度，ξ_n为空间频率；Ψ_g表示地表散射的小波稀疏基，采用单位正交基，Ψ_v表示冠层散射的小波稀疏基；

p_g、p_v分别为地表散射和冠层散射的后向散射功率矢量；

σ_g、σ_v分别为地表散射和冠层散射对应的噪声功率，t_g＝vec(σ_gI)、t_v＝vec(σ_vI)分别为对应噪声功率矩阵的矢量化操作，I为单位矩阵；δ_g和δ_v分别为地表散射和冠层散射的协方差矩阵与对应真实协方差矩阵之间的估计误差。

3.根据权利要求2所述的方法，其特征在于，步骤4中采用EM算法迭代更新稀疏贝叶斯学习模型中的超参数并获得冠层散射和地表散射的层析谱，具体为：

步骤a，定义先验分布：

为了便于阐述，将步骤3得到的地表散射和冠层散射的稀疏贝叶斯学习模型在形式上统一表示为：Y＝Gx+δ；其中，G＝[BΨ^H,I]表示在稀疏基下新的测量矩阵，x＝[Ψp,t]为待求稀疏域下的后向散射功率和噪声矢量，Y为地表散射或冠层散射的干涉协方差矩阵的矢量表示；

根据待求稀疏域下的后向散射功率和噪声矢量x是稀疏这一先验信息，定义其服从高斯分布，即x～N(0,Λ)；其中Λ＝diag(Γ)，

为稀疏贝叶斯学习模型的超参数矢量，表征了x的空间分布；前D个超参数τ₁,...τ_D控制高度向上D个散射体后向散射功率的稀疏性，后N²个超参数

为***噪声；则Y关于Γ的概率密度函数为：

p(Y|Γ)＝∫p(Y|x)p(x|Γ)dx；

步骤b，参数初始值定义：

设置稀疏贝叶斯学习模型的超参数向量Γ初始值和迭代终止条件；

步骤c，计算后验概率分布：

基于贝叶斯准则和步骤a定义的先验分布，计算x的均值矢量和协方差矩阵：

Σ_Y＝V+GΛ⁽ⁱ⁾G^H

式中，上标i表示当前迭代次数，μ为x的均值矢量，Σ_x和Σ_Y分别为x和Y的协方差矩阵，V为样本协方差矩阵估计误差δ的二阶统计；

步骤d，基于步骤c得到的均值矢量和协方差矩阵，采用EM算法学习更新超参数向量：

步骤e，重复步骤c、d，直至满足迭代终止条件，使用超参数向量Γ中的前D个超参数Γ₍₁…D)计算冠层散射和地表散射的层析谱：

4.根据权利要求3所述的方法，其特征在于，迭代终止条件为||Γⁱ⁺¹-Γⁱ||₂/||Γⁱ||₂≤tol，和/或达到最大迭代次数。

5.根据权利要求1所述的方法，其特征在于，根据地表散射的层析谱提取森林反演区域的林下地形的方法为：提取地表散射的层析谱中功率最大值所在的位置作为地面高度的估计值。

6.根据权利要求3所述的方法，其特征在于，根据冠层散射的层析谱提取森林反演区域的植被高度的方法为：通过功率衰减法确定冠层顶部的位置，即：以冠层相位中心处为起始点，以预设的步长向上衰减，当功率衰减处提取的高度值和Lidar数据提取的冠层高度DSM的均方根误差达到最小时，将此时的高度值作为冠层高度的估计值；通过冠层高度减去地面高度即为植被高度。

7.一种基于散射机制分解的植被垂直结构及林下地形反演装置，其特征在于，包括：

数据预处理模块，用于：获取森林反演区域多基线多极化的SAR数据集，并进行配准、去平地，得到极化层析SAR数据；

协方差计算分解模块，用于：对数据预处理模块得到的极化层析SAR数据计算其多基线多极化协方差矩阵，并分解计算冠层散射和地表散射的干涉协方差矩阵；

其中，多基线多极化协方差矩阵的计算方法与分解方法为：

式中，y为步骤1得到的极化层析SAR数据，

g_HV＝g_VH，g_HH、g_HV、g_VH、g_VV分别为HH、HV、VH、VV极化方式对应的SAR数据构成的N维SAR数据观测值；(*)^H为共轭转置操作，E[·]为求期望运算；W为计算得到的多基线多极化协方差矩阵；C_g和R_g分别为地表散射的极化协方差矩阵和干涉协方差矩阵，C_v和R_v分别为冠层散射的极化协方差矩阵和干涉协方差矩阵，

为克罗内克积；

对多基线多极化协方差矩阵W进行SVD分解，使用分解得到的特征矢量计算冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g：

首先，取SVD分解得到的前两个最大特征值对应的特征向量重塑矩阵得到R₁和R₂，使用R₁和R₂表示冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g为：

式中，a、b为一对实数；

然后，基于冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g的半正定约束，即R_v和R_g的特征值非负，求解实数a、b的取值区间；

再后，根据以下散射机制最大化分离的准则，求解实数a、b：

其中，trace()为求矩阵的迹，||·||_F为矩阵的Frobenius范数；

最后，根据求解得到的实数a、b和矩阵R₁和R₂，即可得到冠层散射的干涉协方差矩阵R_v和地表散射的干涉协方差矩阵R_g；

模型构建模块，用于：针对冠层散射和地表散射的干涉协方差矩阵，均进行矢量化和稀疏表达，分别得到冠层散射和地表散射的稀疏贝叶斯学习模型；

层析谱计算模块，用于：采用EM算法迭代更新两个稀疏贝叶斯学习模型中的超参数，进而获得冠层散射和地表散射的层析谱；

植被垂直结构及林下地形反演模块，用于：根据地表散射的层析谱提取森林反演区域的林下地形，根据冠层散射的层析谱提取森林反演区域的植被高度。

Images