CN110988879B

CN110988879B - 一种植被参数反演方法、终端设备及存储介质

Info

Publication number: CN110988879B
Application number: CN201911345531.3A
Authority: CN
Inventors: 朱建军; 张兵; 付海强; 李志伟; 汪长城
Original assignee: Central South University
Current assignee: Central South University
Priority date: 2019-12-24
Filing date: 2019-12-24
Publication date: 2022-08-12
Anticipated expiration: 2039-12-24
Also published as: CN110988879A

Abstract

本发明公开了一种植被参数反演方法、终端设备及存储介质，其中方法包括以下步骤：步骤1，对2幅SAR影像进行预处理和极化干涉处理，以获取单基线全极化复相干系数观测值；步骤2，确定参与联合解算的区块大小，即联合解算的像素点的个数，并设置未知参数反演的初始值；步骤3，基于RVoG模型与各像素的全极化复相干系数，构建观测方程；步骤4，根据前述获得的参数初值，采用非线性迭代算法进行植被高度反演。本发明解决了现有技术中应用RVoG模型反演植被高度时的秩亏问题，提高参数反演求解时的解算稳定性。

Description

一种植被参数反演方法、终端设备及存储介质

技术领域

本发明涉及极化干涉合成孔径雷达(PolInSAR)在植被参数反演领域，具体涉及一种联合多像素的单基线极化干涉SAR植被高度反演方法、终端设备及存储介质。

背景技术

地球陆表约32.3％的面积被植被覆盖。准确提取森林高度对于研究生态环境变化、全球碳循环和森林资源管理具有重要意义。PolInSAR技术能够在一个分辨率单元中分辨出不同散射机制的相位中心，理论上可以通过确定表面散射和体积散射的相位中心来反演森林高度。近年来，PolInSAR的发展为在全球范围内提取森林高度提供了契机。

在采用PolInSAR技术进行植被高度反演领域，目前应用最为广泛的模型为随机地体二层相干散射(Random Volume over Ground,RVoG)模型，其建立了InSAR观测值与植被参数之间的函数关系。RVoG模型由于建模过程中存在过度参数化的情况导致未知数个数大于独立观测值个数，未知参数解算秩亏。

当仅存在单基线PolInSAR数据时，基于RVoG模型反演植被参数的观测信息更少，导致未知参数解算秩亏。因此，有必要设计一种仅存在单基线PolInSAR数据时的植被参数反演方法。

发明内容

本发明所要解决的技术问题在于，当仅存在单基线PolInSAR数据时，提供一种联合多像素的PolInSAR植被参数反演方法、终端设备及存储介质，可以解决基于RVoG模型的解算过程中观测信息不足导致的未知参数解算的秩亏问题。

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

一种植被参数反演方法，包括以下步骤：

步骤1，对2幅SAR影像进行配准、去平地、多视以及极化干涉处理，生成单基线全极化复相干系数γ_(ω)；

步骤2，设置植被参数反演的细分区块大小为N行M列，获取区块每个像素点的复相干系数，与RVoG模型进行匹配，得到观测方程为：

其中，ω表示与基线和极化方式相关的极化状态；

表示地表相位，是有关地表高程h_g的函数；μ_(ω)表示极化状态为ω所对应的地体幅度比；γ_v表示纯体去相干系数，是有关植被高度h_v和消光系数σ的函数；k_z为垂直向有效波数；

针对每个像素的多种极化方式，均可得到至少2个复数观测方程，则区块NM个像素点可得到至少2NM个复数观测方程，即4NM个实数观测方程；

步骤3，计算区块内每个像素点的参数初值：将地体幅度比μ_(ω)的归一化参数x的初值设为0，然后利用RVoG模型和三阶段植被高度反演算法计算植被高度初值、消光系数初值和地表高程初值；

其中，归一化参数x是指：在复平面内，体散射占优通道的极化复相干系数到纯体去相干系数的距离；

步骤4，根据步骤2得到的4NM个实数观测方程以及步骤3得到的各初值，采用非线性迭代算法求解观测方程中的2NM+2个未知参数：植被高度、消光系数、地表高程、地体幅度比归一化参数；

其中，区块内NM个像素点所对应的植被高度和消光系数均恒定，地表高程和地体幅度比的归一化参数x均随像素点的变化而变化。

进一步地，利用RVoG模型和三阶段植被高度反演算法计算植被高度初值的方法为：

步骤1)，直线拟合：

对于区块内每个像素点，均利用在该基线2种极化状态下的复相干系数观测值，对RVoG模型在复平面内进行直线拟合，得到公式(3)所示直线：

步骤2)，计算地表相位：

对于区块内每个像素点，计算公式(3)所示直线与复平面单位圆的2个交点A和B，将其中与体散射占优极化方式点之间的距离大于与地面散射占优极化方式点之间的距离的交点，即为地表相位点

然后利用地表相位

按以下公式(4)计算地表高程初值：

进一步地，利用RVoG模型和三阶段植被高度反演算法计算消光系数初值、地表高程初值的方法为：

选择区块的中心像素点计算整个区块的植被高度和消光系数初始值：

设所选择像素点的体散射占优极化通道仅包含植被层散射贡献，其地体幅度比为0，则纯体相干性的表达式为

表示体散射占优通道复相干系数，是有关植被高度h_v以及消光系数σ的函数，表示为

采用查找表的方法，通过给定一系列植被高度h_v以及消光系数σ，根据公式(5)构建二维查找表，分别计算

然后在表格中找到

与

相差最小的一组，作为植被高度初值和消光系数初值：

进一步地，区块的行数N和列数M均为大于1的奇数。

本发明还提供一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述任一项所述的方法。

本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现上述任一项所述的方法。

有益效果

本发明的有益效果包括：

1、对于单基线数据，联合多像素进行参数反演，解决了现有技术中应用RVoG模型反演植被高度时的秩亏问题，提高参数反演求解时的解算稳定性；

2、植被参数反演方法简单明了、易于实现，未来可以广泛应用于长时间基线下的大范围植被高度、全球碳储量估计等多方面。

附图说明

图1为本发明植被高度反演方法的流程示意图；

图2为RVoG模型框架下的地体幅度比几何表示图；

图3为本发明实施例的植被高度反演结果；

图4为本发明实施例植被高度反演结果与LiDAR测量结果交叉验证图。

具体实施方式

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

本发明提供一种联合多像素的单基线极化干涉SAR植被高度反演方法，如图1所示，包括以下步骤：

步骤1，极化干涉处理；

此步骤主要获取单基线全极化复相干系数观测值，包括对覆盖植被区域的主、辅2幅SAR影像进行配准、去平地相位、多视去除斑点噪声和极化干涉处理几个步骤，获取覆盖植被区域的单基线全极化复相干系数观测值。

其中：

ω表示与基线和极化方式相关的极化状态；

μ_(ω)表示极化状态为ω所对应的地体幅度比，不同像素的同一种极化方式、同一像素的不同种极化方式均存在差异；

表示地表相位，与地表高程和基线几何相关：

h_g：地表高程，本发明中任意像素的地表高程可能存在差异；

k_z：垂直向有效波数，与干涉基线几何构型有关，其具体表达如下：

其中：

B_⊥：垂直基线长度；

R：雷达卫星到观测目标之间的距离；

λ：微波波长；

γ_v表示纯体去相干系数，具体表达式如下：

其中：

h_v：植被高度，本发明区块中NM个像素的植被高度保持一致；

σ：植被体平均消光系数；在RVoG模型中，设植被体为均匀介质，微波穿透植被过程中服从指数分布衰减，故植被体消光系数为固定的实数；本专利认为，区块中NM个像素的消光系数保持一致，不同区块消光系数可能存在变化；

θ：微波入射角。

由于每个像素的2种极化方式，均可得到2个如通式(6)所示的复数观测方程，因此区块NM个像素点可得到2NM个复数观测方程；再拆分复数的实部和虚部，区块NM个像素点可得到4NM个实数观测方程。即4NM个实数观测方程。

步骤2.5，简化RVoG模型，减少观测方程中的参数个数；

对于单基线干涉图内任意像素点，将极化相关的地体幅度比μ_(ω)进行归一化处理，减少至唯一的归一化参数x；其中，归一化参数x是指，在复平面内，该像素中体散射占优极化方式的复相干系数到纯去相干系数的距离，具体方式如下：

如图2，每个像素在RVoG模型框架下的相干迹，在复平面单位圆内均是一条直线，地体幅度比

的作用为确定不同极化方式的复相干性在相干直线上的位置，显然其存在一定的线性相关。因此，地体幅度比表示为公式(11)：

其中：

任意极化方式的复相干性系数；

体散射占优方式的复相干系数；

γ_v：纯体去相干系数；

x：体散射占优方式的复相干系数

到纯体去相干系数γ_v的距离。

其中，公式(11)中的

D均可以通过PolInSAR估计的复相干系数进行求解，因此通过公式(11)可将区块中第i行第j列像素点的2种极化方式的地体幅度比归一化为未知参数x^ij的表达式。

步骤3，计算区块内每个像素点的参数初值：

由于理想状态下的地体幅度比的归一化参数为0，且本发明认为不同像素点的同一种极化方式、同一像素点的不同种极化方式均存在差异，故在设置初值时，将每个像素点的地体幅度比归一化参数初值均设置为0；

并利用RVoG模型和三阶段植被高度反演算法计算植被高度初值、消光系数初值和地表高程初值，具体方法如下：

步骤1)，直线拟合：

对于区块内每个像素点，均利用在该基线2种极化状态下的复相干系数观测值，对RVoG模型在复平面内进行直线拟合，得到公式(12)所示直线：

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

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

其中：

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

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

然后利用每个像素点的地表相位

按以下公式(14)计算该像素点的地表高程初值：

步骤3)，植被高度与消光系数估计：

本发明认为，区块中NM个像素点的植被高度和消光系数保持一致，不同区块间存在变化。本发明选择区块的中心像素点的植被高度以及消光系数初值，作为整个区块的植被高度和消光系数初值：

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

其中，区块中心像素点的选择是指选择区块处于最中心位置的像素点：如果最中心为1个像素点，则选择该像素点的复相干系数作为观测值以建立观测方程；如果最中心位置为2个像素点，则选择该2个像素点的复相干系数的平均值作为观测值以建立观测方程；如果最中心位置为4个像素点，则选择该4个像素点的复相干系数的平均值作为观测值以建立观测方程。

在本实施例中，优选区块的行数N和列数M均为大于1的奇数，满足区块最中心为1个像素点，可以避免复相干系数的平均值计算。

步骤4，联合多像素的单基线极化干涉SAR植被高度反演；

根据步骤2得到的4NM个实数观测方程以及步骤3得到的各初值，采用非线性迭代算法求解观测方程中的2NM+2个未知参数：植被高度、消光系数、地表高程、地体幅度比归一化参数。

本发明的理论上，RVoG模型相干迹为复平面单位圆内的一条直线。对于区块的任意像素，均可通过纯体散射点和地表相位点可唯一确定这条直线在复平面单位圆内的位置和方向：

(1)对于纯体散射点，

可以提供实部和虚部2个独立的实数观测值；

(2)对于地表相位点，

由于地表相位点的实部与虚部的平方和恒为1，因此只能提供1个独立的实数观测值。

因此对于任意像素，无论存在多少个极化复相干系数，即采用多少种极化方式，RVoG模型本身只能提供3个独立观测值，因此上述4NM个实数观测方程中，独立观测方程的个数为3NM。

本发明认为，区块内NM个像素点所对应的植被高度和消光系数均恒定，因此3NM个独立观测方程中包括2个固定的未知参数；地表高程和地体幅度比的归一化参数x均随像素点的变化而变化，因此3NM个独立观测方程中还包括2NM个未知参数。因此，3NM个独立观测方程中共有2+2NM个未知参数，当独立观测方程的个数大于未知参数的个数时，即3NM≥2+2NM时，即可根据步骤2得到的各参数初值、采用非线性迭代算法解算3NM个独立观测方程，最终得到2+2NM个未知参数：区块恒定的植被高度和消光系数、区块每个像素的地表高程和地体幅度比归一化参数。

为进一步清晰说明本发明联合多像素点的单基线极化干涉SAR植被高度反演方法，采用BioSAR 2008项目提供的瑞典北部Krycklan地区(64°14′N，19°46′E)的两幅全极化L波段SAR影像对本文算法进行验证。注意，此处仅用于举例说明，本发明并不限定数据源。

实验区域海拔在150m-380m米；植被主要以人工林为主，树种有云杉、松树、桦树等，植被高度平均18m。

采用E-SAR L波段SAR数据进行实验。实验数据采用2景机载P波段全极化数据，以编号0201影像为主影像、编号0205影像为辅影像组成干涉对。利用欧空局发布的PolSARpro软件进行数据预处理。

具体步骤如下：

S1，极化干涉处理：

以编号0201影像为主影像、编号0205影像为辅影像组成干涉对，仅有一条单基线，应用PolSARpro软件分别对主辅影像进行去除平地效应、多视处理(方位向2：距离向1)，然后获取极化复相干系数。

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

S2，生成参数初值：

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

首先设置植被参数反演的细分区块大小为N行M列(此处设置N＝M＝3，其他区块大小合理即可，不限于N＝M＝3)，将每个像素的地体幅度比归一化参数初值设置为0。然后通过目前PolInSAR植被参数反演领域应用较为广泛的三阶段算法，来计算植被高度、消光系数、地表高程的初值，具体方法如下：

1)直线拟合：在RVoG模型框架下，选择PDHigh和PDLow两种极化方式对区块中9个像素点分别进行直线拟合。

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

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

每个像素的地表相位确定后，则可以直接通过以下公式计算每个像素的地表高程h_g作为非线性迭代算法的初值：

3)植被高度与消光系数估计：

本专利假设，区块9个像素点内，植被高度与消光系数恒定相同。本发明计算区块中心像素点(即区块第2行第2列的像素点)的植被高度与消光系数值作为非线性迭代算法的初值。假设中心像素点的体散射占优极化通道γ_PDHigh仅包含植被层散射贡献，其地体幅度比为0，则其纯体相干性的表达式为

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

S3，联合多像素的单基线极化干涉SAR植被高度反演方法

根据RVoG模型和本发明假设，并对地体幅度比μ_(ω)进行归一化为地体幅度比归一化参数x后，构建非线性观测方程反演待求参数：

其中，上标11代表区块中第一行第一列个像素，以此类推，上标33代表区块中第三行第三列个像素。本实例中，采用3×3＝9个像素进行多像素联合解算植被高度。

从公式(19)可以得到，本实施例在参数反演的未知参数为固定的σ,h_v以及9个像素的地表高程

和地体幅度比归一化因子x¹¹,x¹²,x¹³,…,x³³，共20个未知参数。利用步骤2中的植被高度初值、消光系数初值、地表高程初值、地体幅度比归一化参数初值，并采用非线性迭代算法，最终得到待求取的植被高度、消光系数、地表高程、地体幅度比归一化参数等植被参数。

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

本发明实施例还提供一种终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述实施例所述方法。

本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现上述实施例所述方法。

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

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