一种基于两景SAR影像的植被高度反演方法及装置与流程

本发明涉及极化合成孔径雷达干涉测量(polinsar)在植被参数反演领域，具体涉及一种基于两景sar影像和傅里叶-勒让德多项式植被高度反演模型的植被高度反演方法及装置。

背景技术：

森林高度以及植被垂直结构是研究全球碳循环、生态环境变化、森林蓄积量等的重要指标。近年来，极化合成孔径雷达干涉测量(polinsar)技术快速发展，其可以判断同一分辨单元内不同散射机制的相位中心位置。理论上，polinsar可以通过体散射和表面散射相位中心在植被垂直方向上的高度做差来提取森林高度。

目前应用polinsar技术进行植被高度提取应用最为广泛的散射模型为1996年提出的随机地体二层散射(randomvolumeoverground,rvog)模型。但是，rvog模型在建模时，认为植被体为一个均匀的散射体，这显然与现实情况中复杂的植被垂直结构不符。

2006年，极化相干层析(polarizationcoherencetomography，pct)技术的提出，为描述复杂的植被垂直结构带来了可能，pct技术在单基线情况下，只能展开到2阶多项式，如需展开到更高阶次，只能采用多景覆盖同一区域到sar影像构成多基线数据；同时，应用pct技术时，需要准确的植被高度和地表相位先验信息。

专利申请文献cn109738895a，其中公开的基于傅里叶-勒让德多项式的植被高度反演模型，在具有多景sar影像情况下，可以描述复杂的植被垂直结构。但是在仅有2景sar影像情况下，此模型只能展开到0阶多项式，相当于将植被垂直结构描述为固定的常数模型，这显然与现实不符。而就目前技术来说，覆盖同一区域的多景sar影像的准确获取还比较困难。

因此，针对仅有两景sar影像的情况下，有必要设计一种基于两景sar影像和傅里叶-勒让德多项式植被高度反演模型的植被高度反演方法。

技术实现要素：

本发明要解决的技术问题在于，提供一种基于两景sar影像和傅里叶-勒让德多项式的植被高度反演方法及装置，解决传统单基线傅里叶-勒让德多项式模型只能展开到0阶多项式，将植被垂直结构描述为固定的线性模型，与现实不符的问题。

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

一种基于两景sar影像的植被高度反演方法，包括以下步骤：

步骤s10，获取关于植被覆盖层垂直结构的两幅sar影像，进行配准、去平地、多视以及极化干涉处理，生成p个极化复相干系数ω表示极化矢量，p表示极化方式；；

步骤s20，设置植被参数反演的空间域为n行m列个像素点，且空间域内所有像素点对应的植被高度hv、地表高程hg相同；使用傅里叶-勒让德多项式表示植被垂直结构，除傅里叶-多项式最后一阶展开项的系数不同外，其他傅里叶-勒让德多项式展开项的系数均相同；

步骤s30，使用傅里叶-勒让德多项式表示植被垂直结构，并将空间域内像素点的极化复相干系数γ(ω)作为观测值，构建基于傅里叶-勒让德多项式的植被高度反演模型；空间域共得到mn个观测方程，且观测方程的表达式为：

式中，q为傅里叶-勒让德多项式展开项的阶数，i为复数虚部标识，hg表示地表高程，kv为中间参数，且有kv＝kzhv/2，hv表示植被高度；kz表示垂直向有效波数，且区域内每个像素点对应的垂直向有效波数不同；f0,f1,f2,f3,…分别为傅里叶-勒让德多项式的第0、1、2、3、……个展开项，均是有关中间参数kv的已知函数，a10,a20,a30,…分别为f1,f2,f3,…的系数；

步骤s40，利用rvog模型和三阶段植被高度反演算法，计算空间域内植被高度hv的初值和地表高程hg的初值；再利用步骤s30得到的观测方程，以及植被高度hv的初值和地表高程hg的初值，计算傅里叶-勒让德多项式各展开项系数的初值；

步骤s50，将步骤s40得到的植被高度、地表高程以及傅里叶-勒让德多项式各展开项系数的初值，均作为每个像素点观测方程中的对应参数初值，并采用非线性迭代算法求解观测方程中的未知参数：植被高度、地表高程以及傅里叶-勒让德多项式各展开项系数。

在更优的技术方案中，利用rvog模型和三阶段植被高度反演算法计算地表高程初值的方法为：

步骤1)，直线拟合：

对于mn个像素构成的空间域中每个像素点，均利用p个极化复相干系数观测值，对rvog模型在复平面内进行直线拟合，得到公式(3)所示直线：

表示地表相位，是有关地表高程hg的函数；μ(ω)表示极化矢量ω所对应的地体幅度比；γv表示纯体去相干系数，是有关植被高度hv和消光系数σ的函数；

步骤2)，计算地表相位和地表高程初值：

对于空间域内每个像素点，计算公式(3)所示直线与复平面单位圆的2个交点a和b，将其中与体散射占优极化方式点之间的距离大于与地面散射占优极化方式点之间的距离的交点，即为地表相位点然后利用地表相位按以下公式(4)计算地表高程初值：

步骤3)，将空间域mn个像素点的地表高程初值的平均值作为空间域的地表高程初值。

在更优的技术方案中，利用rvog模型和三阶段植被高度反演算法计算植被高度初值的方法为：

步骤1)，设像素点的体散射占优极化通道仅包含植被层散射贡献，其地体幅度比为0，则纯体相干性的表达式为表示体散射占优通道复相干系数，是有关植被高度hv以及消光系数σ的函数，表示为

采用查找表的方法，通过给定一系列植被高度hv以及消光系数σ，根据公式(5)构建二维查找表，分别计算然后在表格中找到与相差最小的一组，作为像素点对应的植被高度初值和消光系数初值：

步骤2)，将空间域mn个像素点的植被高度初值的平均值作为空间域的植被高度初值。

在更优的技术方案中，n和m均大于1。

在更优的技术方案中，傅里叶-勒让德多项式的最高展开阶数为(2pmn-2)/p。

在更优的技术方案中，傅里叶-勒让德多项式的第0、1、2、3、4、5、6个项展开项分别为：

本发明还提供一种基于两景sar影像的植被高度反演装置，包括观测值生成模块、空间域设置模块、观测方程构建模块、参数初值获取模块和未知参数计算模块；

所述观测值生成模块用于：获取关于植被覆盖层垂直结构的两幅sar影像，进行配准、去平地、多视以及极化干涉处理，生成极化复相干系数ω表示极化矢量，p表示极化方式；设共采用p种极化方式进行干涉处理，得到p个极化复相干系数；

所述空间域设置模块用于：设置植被参数反演的空间域为n行m列个像素点，且空间域内所有像素点对应的植被高度hv、地表高程hg以及植被体垂直结构相同；

所述观测方程构建模块用于：使用傅里叶-勒让德多项式表示植被垂直结构，并将空间域内像素点的极化复相干系数γ(ω)作为观测值，构建基于傅里叶-勒让德多项式的植被高度反演模型；空间域共得到mn个观测方程，且观测方程的表达式为：

式中，q为傅里叶-勒让德多项式展开项的阶数，i为复数虚部标识，hg表示地表高程，kv为中间参数，且有kv＝kzhv/2，hv表示植被高度；kz表示垂直向有效波数，且区域内每个像素点对应的垂直向有效波数不同；f0,f1,f2,f3分别为傅里叶-勒让德多项式的第0、1、2、3项，均是有关中间参数kv的已知函数，a10,a20,a30分别为f1,f2,f3的系数；

所述参数初值获取模块用于：利用rvog模型和三阶段植被高度反演算法，计算空间域内植被高度hv的初值和地表高程hg的初值；再利用观测方程构建模块得到观测方程，以及植被高度hv的初值和地表高程hg的初值，计算傅里叶-勒让德多项式各展开项系数的初值；

所述未知参数计算模块用于：将参数初值获取模块得到的植被高度、地表高程以及傅里叶-勒让德多项式各展开项系数的初值，均作为每个像素点观测方程中的对应参数初值，并采用非线性迭代算法求解观测方程中的未知参数：植被高度、地表高程以及傅里叶-勒让德多项式各展开项系数。

有益效果

本发明的有益效果包括：

1、本发明将植被垂直结构使用多阶傅里叶-勒让德多项式描述，可以精确地描述植被层的垂直结构，更加符合现实场景，提高植被高度反演的准确性；

2、正因为本发明采用傅里叶-勒让德多项式描述的垂直结构更符合现实场景，因此本发明的植被高度反演方法也就可适用于各种植被垂直结构进行植被高度反演，因而本发明的适用范围广，应用更灵活；

3、通过设置区域大小，并根据植被垂直结构的特点设定区域内所有像素点对应的植被高度、地表高程以及植被体垂直结构相同，从而使本发明建立的傅里叶-勒让德多项式植被高度反演模型，在仅有两景sar影像的情况下仍能解算得到植被高度，因此本发明实现容易，可以更好地服务于林业估产、全球碳循环研究等。

附图说明

图1为本发明实施例所述方法的流程示意图；

图2为本发明实施例所述方法的的植被高度反演结果；

图3为本发明实施例所述方法的植被高度反演结果与lidar测量结果交叉验证图。

具体实施方式

下面对本发明的实施例作详细说明，本实施例以本发明的技术方案为依据开展，给出了详细的实施方式和具体的操作过程，对本发明的技术方案作进一步解释说明。

本实施例提供一种基于两景sar影像和傅里叶-勒让德多项式的植被高度反演方法，如图1所述，包括以下步骤：

步骤s10，对2幅sar影像进行预处理，获取全极化复相干系数观测值；

首先获取关于植被覆盖层垂直结构的2幅sar影像，然后依次进行配准、平地相位去除、多视处理、极化干涉等步骤，最终获取sar影像每个像素点的p个复相干系数ω表示极化矢量，p表示极化方式。

步骤s20，设置植被参数反演的空间域为n行m列个像素点，且空间域内所有像素点对应的植被高度hv、地表高程hg相同；使用傅里叶-勒让德多项式表示植被垂直结构，除傅里叶-多项式最后一阶展开项的系数不同外，其他傅里叶-勒让德多项式展开项的系数均相同。

步骤s30，构建基于两景sar影像空间域的傅里叶-勒让德多项式植被高度反演模型；

传统rvog模型认为植被体为均匀层，将植被垂直结构描述为指数函数，这显然不符合现实复杂的植被场景。因此本实施例将植被垂直结构使用傅里叶-勒让德多项式描述，可以精确地描述植被层的垂直结构，更加符合现实场景。然后再将空间域内像素点的极化复相干系数作为观测值，构建基于傅里叶-勒让德多项式的植被高度反演模型；空间域共得到mn个复数观测方程，且观测方程的表达式为，以3阶傅里叶-勒让德多项式为例：

式中，i为复数虚部标识，hg表示地表高程；

kz表示垂直向有效波数，具体：kz＝4πb⊥/λrsinθ，b⊥表示垂直基线长度，λ表示微波波长，r表示雷达到目标的距离，θ表示入射角；由于区域内每个像素点，其对应的距离r和入射角θ不同，因此其垂直向有效波数不同；

kv为中间参数，没有物理意义，具体：kv＝kzhv/2，hv表示植被高度；

f0,f1,f2,f3，…分别为傅里叶-勒让德多项式的第0、1、2、3、……个展开项，均是有关中间参数kv的已知函数，a10,a20,a30，…分别为f1,f2,f3，…的系数，

傅里叶-勒让德多项式的第0至6项展开项分别为：

f0＝sinkv/kv

步骤s40，计算空间域范围内植被高度和地表高程初值：

首先，利用rvog模型和三阶段植被高度反演算法，计算空间域内植被高度hv的初值和地表高程hg的初值，具体方法如下：

步骤1)，直线拟合：

对于mn个像素构成的空间域中每个像素点，均利用p个极化复相干系数观测值，对rvog模型在复平面内进行直线拟合，得到公式(15)所示直线：

表示地表相位，是有关地表高程hg的函数；μ(ω)表示极化矢量ω所对应的地体幅度比；γv表示纯体去相干系数，是有关植被高度hv和消光系数σ的函数；

步骤2)，计算地表相位和地表高程初值：

对于空间域范围内每个像素点，通过步骤1)拟合得到的直线与复平面单位圆均存在2个交点，其中1个为像素点对应的地表相位点。根据地表相位点距离体散射占优极化方式点的距离要大于距离地面散射占优极化方式点的距离这一原则，通过判断公式(13)判断得到地表相位点，具体：

其中：分别代表体散射占优通道复相干系数、地面散射占优通道复相干系数。

因此，计算每个像素点的地表相位的方法为：计算公式(12)所示直线与复平面单位圆的2个交点a和b，将其中与体散射占优极化方式点之间的距离大于与地面散射占优极化方式点之间的距离的交点，即为该像素点的地表相位点

再利用每个像素点的地表相位按以下公式(14)计算该像素点的地表高程初值：

步骤3)，植被高度估计：

设像素点的体散射占优极化通道仅包含植被层散射贡献，其地体幅度比为0，则纯体相干性的表达式为表示体散射占优通道复相干系数，是有关植被高度hv以及消光系数σ的函数，表示为

采用数值计算方法，通过给定一系列植被高度hv以及消光系数σ，根据公式(13)建立二维查找表，找到差异最小的一组值，即为植被高度hv以及消光系数σ的估计值：

步骤4)，空间域范围内植被高度和地表高程初值估计：

将空间域mn个像素点的地表高程初值的平均值和植被高度初值的平均值，分别作为空间域的地表高程初值和植被高度初值。

然后，利用步骤s30得到的观测方程，以及前面得到的植被高度初值和地表高程初值，计算傅里叶-勒让德多项式各展开项系数的初值。具体计算可将公式(13)拆分实部和虚部并转换为如下公式(19)所示，以3阶傅里叶-勒让德多项式为例：

根据前述获得的植被高度初值和地表高程初值，公式(19)中的未知数仅为傅里叶-勒让德多项式各展开项系数通过现有已知的矩阵除法即可计算得到其中的各系数初值，在此不再赘述。

步骤s50，基于空间域的单基线极化干涉sar植被高度反演；

将公式(13)所示的观测方程拆分为实部和虚部，得到公式(20)所示的观测方程组，以3阶傅里叶-勒让德多项式为例：

将步骤s40得到的植被高度、地表高程以及傅里叶-勒让德多项式各展开项系数的初值，均作为每个像素点观测方程中的对应参数初值，并采用非线性迭代算法求解观测方程中的未知参数：植被高度、地表高程以及傅里叶-勒让德多项式各展开项系数。

由于采用p种极化方式进行极化干涉处理，因此：区域内每个像素点对应得到实部与虚部共2p个观测值，区域内nm个像素点即得到2pnm个观测值；再由于每个像素点的傅里叶-勒让德多项式对应展开项不同，因此2pnm个观测值可对应构建2pnm个观测方程。

本发明实施例区域内所有像素点对应的植被高度hv、地表高程hg以及植被体垂直结构相同，故有：(1)在同一极化方式下，所有像素点的傅里叶-勒让德多项式的对应项系数中，除最后一阶展开项的系数不同外，其他展开项的系统都相同，以3阶傅里叶-勒让德多项式为例，所有像素点的a10相同，所有像素点的a20相同，所有像素点的a30不同；(2)在不同极化方式下，傅里叶-勒让德多项式的对应项系数不同。因此2pnm个独立的观测方程中共有2+p(q-1)+pnm个未知参数(q为傅里叶-勒让德多项式的阶数)：植被高度hv、地表高程hg、p种极化方式的傅里叶-勒让德多项式各展开项系数。只要满足独立观测方程的数量大于未知参数的数量，即2pnm＞＝2+p(q-1)+pnm，即可求解得到植被高度、地表高程以及傅里叶-勒让德多项式各展开项系数。因此，在实际应用时，只要满足条件2pnm＞＝2+p(q-1)+pnm，均可以根据需要选择展开项的阶数q。

本发明实施例还提供一种植被高度反演装置，包括观测值生成模块、空间域设置模块、观测方程构建模块、参数初值获取模块和未知参数计算模块；其中，

所述观测值生成模块用于：获取关于植被覆盖层垂直结构的两幅sar影像，进行配准、去平地、多视以及极化干涉处理，生成极化复相干系数ω表示极化矢量，p表示极化方式；设共采用p种极化方式进行干涉处理，得到p个极化复相干系数；

所述空间域设置模块用于：设置植被参数反演的空间域为n行m列个像素点，且空间域内所有像素点对应的植被高度hv、地表高程hg以及植被体垂直结构相同；

所述观测方程构建模块用于：使用傅里叶-勒让德多项式表示植被垂直结构，并将空间域内像素点的极化复相干系数γ(ω)作为观测值，构建基于傅里叶-勒让德多项式的植被高度反演模型；空间域共得到mn个观测方程，且观测方程的表达式为，以3阶傅里叶-勒让德多项式为例：

式中，i为复数虚部标识，hg表示地表高程，kv为中间参数，且有kv＝kzhv/2，hv表示植被高度；kz表示垂直向有效波数，且区域内每个像素点对应的垂直向有效波数不同；f0,f1,f2,f3分别为傅里叶-勒让德多项式的第0、1、2、3项，均是有关中间参数kv的已知函数，a10,a20,a30分别为f1,f2,f3的系数；

所述参数初值获取模块用于：利用rvog模型和三阶段植被高度反演算法，计算空间域内植被高度hv的初值和地表高程hg的初值；再利用观测方程构建模块得到观测方程，以及植被高度hv的初值和地表高程hg的初值，计算傅里叶-勒让德多项式各展开项系数的初值；

所述未知参数计算模块用于：将参数初值获取模块得到的植被高度、地表高程以及傅里叶-勒让德多项式各展开项系数的初值，均作为每个像素点观测方程中的对应参数初值，并采用非线性迭代算法求解观测方程中的未知参数：植被高度、地表高程以及傅里叶-勒让德多项式各展开项系数。

为进一步清晰说明本发明一种基于基于两景sar影像和傅里叶-勒让德多项式植被高度反演模型的植被高度反演方法，采用afrisar项目提供的非洲加蓬mabounie地区的两幅全极化p波段sar影像对本专利模型和算法进行验证。注意，此处仅用于举例说明，本发明并不限定数据源。实验区域为非洲热带森林区域，区域植被覆盖主要为成熟的热带森林和部分退化的森林，森林高度主要在40米到50米之间。

实验数据采用2景机载p波段全极化数据组成干涉对，利用欧空局发布的polsarpro软件进行数据预处理。

具体步骤如下：

s1，极化干涉处理：

仅有2景sar影像构成一条单基线，应用polsarpro软件分别对主辅影像进行去除平地效应、多视处理(方位向4：距离向2)，然后获取极化复相干系数。

本实施例采用2个相位最大分离相干最优极化方式，即pdhigh与pdlow极化方式。

s2，生成植被高度和地表高程参数初值：

此步骤主要是为非线性迭代算法提供未知参数的初始值。

首先设置空间域由3行3列个像素点构成(其他空间域合理即可，不限于3行3列)，则根据本实施例，傅里叶-勒让德多项式的植被参数反演模型最多可以展开到9阶，但是本实施例算法只采用3阶傅里叶-勒让德多项式模型进行植被高度反演(阶数只要不大于9即可，不限于3阶)。

通过目前polinsar植被参数反演领域应用较为广泛的三阶段算法，来计算植被高度、地表高程的初值，具体方法如下：

1)直线拟合：在rvog模型框架下，选择2个相位最大分离相干最优极化方式(dhigh与pdlow)对空间域中9个像素点分别进行直线拟合。

2)地表相位估计：步骤1)中拟合得到的直线与复平面单位圆存在2个交点，其中1个为地表相位点。由于在复平面单位圆内，地表相位点距离体散射占优极化方式点的距离要大于距离地面散射占优极化方式点的距离，因此可采用如下判断公式确定每个像素的地表相位

其中，γpdhigh、γpdlow：分别代表体散射占优通道复相干系数、地面散射占优通道复相干系数。

每个像素的地表相位确定后，则可以直接通过以下公式计算每个像素的地表高程hg：

3)植被高度：

采用数值计算方法，实现给定一系列植被高度hv以及消光系数σ，根据公式(24)建立二维查找表，找到差异最小的一组值，即为植被高度hv估计值，公式(24)为：

4)植被高度ini_hv和地表高程初值ini_hg估计：

本实施例将步骤2)和3)空间域中9个像素点得到的植被高度和地表高程的平均值作为非线性迭代算法的植被高度和地表高程初值。

s3，生成傅里叶-勒让德多项式展开项系数的初值

步骤s2已经获取植被高度和地表高程的初始值，代入公式(25)，计算傅里叶-勒让德多项式展开项系数的初值，公式(25)如下：

其中：

傅利叶-勒让德多项式展开项通式为：

其中，上标11代表空间域中第一行第一列个像素，以此类推，上标33代表第三行第三列个像素。

s4，基于两景sar影像的傅里叶-勒让德多项式植被高度反演模型与空间域解算方法

根本实施例假设，基于傅里叶-勒让德多项式植被高度反演模型，构建非线性观测方程反演待求参数：

其中

傅里叶-勒让德多项式展开项通式为：

本实施例在参数反演的未知参数为共24个未知参数。利用步骤s2和s3中的植被高度初值、地表高程初值、傅里叶-勒让德多项式展开项系数的初值，并采用非线性迭代算法，最终得到待求取的植被高度等植被参数。

其中，本发明中采用的非线性迭代算法为现有技术，通过给定一个解空间以及参数初值，从参数初值开始，在给定的解空间内寻找一个最优解。

图2给出了采用本实施例得出的植被高度反演结果。为定量分析，本文采用lidar植被高度产品作为参考。在实验区域内选取51×51像素大小的1170块样地，并计算对应的平均植被高度用于精度验证。图3给出了本发明实施例植被高度反演结果与lidar植被高度产品的样地交叉验图，采用本发明算法得出的植被高度结果，其对应的rmse为4.32米，精度较高。

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