测绘学报 | 张涛：单基线TanDEM-X InSAR相干性反演森林高度

本文内容来源于《测绘学报》2022年第9期（审图号GS京（2022）0824号）

单基线TanDEM-X InSAR相干性反演森林高度

张涛, 朱建军, 付海强, 汪长城     

中南大学地球科学与信息物理学院, 湖南 长沙 410083

基金项目：国家自然科学基金(4181001033;41820104005;41904004;42030112);中南大学研究生自主探索创新项目(2021zzts0829)

摘要：由RVoG模型演变的SINC函数模型运用到TanDEM-X InSAR数据反演森林高度时面临观测信息不足的问题, 且无法适应不同的森林场景。本文针对这一问题, 提出小范围LiDAR数据辅助下基于TanDEM-X InSAR相干性的森林高度反演方法, 此方法无须全局高精度DEM产品提供地表相位信息而仅利用相干性即可获取森林高度信息。为验证此方法的可行性, 采用西班牙境内两个典型试验区的TanDEM-X InSAR数据进行验证, 并利用LiDAR树高数据进行精度评定。结果表明, 单基线TanDEM-X InSAR相干性可用于反演森林高度, 两个试验区的精度分别为2.34 m和1.74 m。

关键词：森林高度    合成孔径雷达干涉测量    TanDEM-X    SINC函数模型    C-SINC模型   

引文格式：张涛, 朱建军, 付海强, 等. 单基线TanDEM-X InSAR相干性反演森林高度[J]. 测绘学报，2022，51(9)：1931-1941. DOI: 10.11947/j.AGCS.2022.20210298

ZHANG Tao, ZHU Jianjun, FU Haiqiang, et al. Forest height inversion with single-baseline TanDEM-X InSAR coherence[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(9): 1931-1941. DOI: 10.11947/j.AGCS.2022.20210298 

阅读全文：http://xb.chinasmp.com/article/2022/1001-1595/20220910.htm

引 言

森林高度是生物量反演、积蓄量估算和碳循环分析的重要输入数据^[1]。合成孔径雷达(synthetic aperture radar, SAR)技术因其具备穿透观测的特点成为获取大范围、高分辨率及高精度森林高度的有力手段^[2-3]。单一SAR观测模式结合干涉测量技术可得到InSAR技术^[4]，该能够极大提高对目标体垂直结构的敏感性。为实现InSAR信号与森林参数之间的关联，文献[5—6]提出的随机地体散射模型(random volume over ground, RVoG)，并被广泛运用。由于RVoG模型的复杂性，InSAR技术无法提供足够的观测量用于模型反演，需借助先验信息。为解决该问题，文献[5, 7—8]引入极化测量技术，即PolInSAR技术。然而，PolInSAR观测数据目前主要通过机载SAR系统获取，无法实现大范围的森林高度信息提取。现有的星载SAR系统大多采用重轨飞行干涉模式，时间去相干影响严重，导致无法提供高质量的PolInSAR数据^[9-10]。

德国宇航局发射的TanDEM-X卫星采用双站干涉模式，不受时间去相干的影响，保证了较高的干涉质量^[11-12]，对大范围尺度的森林高度制图具有巨大潜力^[13]。然而，TanDEM-X系统主要获取单极化InSAR数据(极化方式为HH)^[12-13]，无法提供足够的观测量用于解算RVoG模型^[13-16]，已有研究主要通过引入外部林下地形数据辅助获取森林高度^[13, 16]。但是，很多地区高精度的林下地形数据不可用。为了解决该问题，考虑到X波段的穿透能力较弱，文献[14—15, 17]忽略地表散射贡献，并假设消光系数为0，RVoG模型被简化成SINC函数模型，可基于InSAR相干性实现森林高度反演。但是，SINC函数模型存在以下问题：①对稀疏或低矮植被覆盖区不能忽略地表散射贡献；②X波段SAR信号穿透植被覆盖层过程中衰减现象明显，假设消光系数为0存在不合理性^[18]。针对上述问题，已有研究提出多种改进模型用于补偿SINC函数模型假设带来的估计偏差。在众多模型中，文献[19]提出的C-SINC模型最具代表性，其在SINC函数模型的基础上，引入两个未知参数建立了半经验模型，成为已有改进SINC函数模型的通式表达^[19-21]。基于C-SINC模型，文献[22—23]验证了该模型反演森林高度的可行性。然而，已有工作中C-SINC模型参数的确定主要依靠经验假设，通用性不强，如何借助少量外部数据实现模型参数的可靠估计值得深入研究。

鉴于此，本文基于C-SINC模型建立森林高度与InSAR观测量之间的关联；在此基础之上，引入少量LiDAR树高数据，确定C-SINC模型中未知参数；之后，构建基于单基线TanDEM-X InSAR相干性的森林高度制图的技术框架；最后，利用西班牙地区森林数据对算法进行验证。

1 InSAR技术反演森林高度的模型与方法

1.1 TanDEM-X InSAR散射模型

RVoG模型已被广泛运用于森林高度反演，该模型建立了InSAR观测量与森林生物物理参数的有效关联，将植被覆盖区建模为各向同性的粒子构成的均匀植被层和电磁波信号不可穿透的地表层^[5-6](图 1)。

图 1 InSAR技术测量森林高度时的RVoG模型

Fig. 1 RVoG model in InSAR measurement of forested area

图选项 

RVoG模型中，观测相干性表示为^[7-8]

 (1)

 (2)

 (3)

式中，f(z)为表征体散射贡献的垂直结构函数；μ为依赖于极化方式ω的地体幅度比；φ₀为地表相位；σ为消光系数；h_v为散射体高度，即树高；θ为电磁波信号入射角；Δθ为干涉影像的入射角差异；λ为雷达波工作波长；R为斜距；B_⊥为垂直基线。k_z为有效垂直波数^[24]，单站模式下，d=2；双站模式下，d=1。

利用TanDEM-X单极化干涉数据来反演森林高度时观测信息不足无法支撑RVoG模型解算。当忽略X波段与森林覆盖区的地表散射贡献(μ≈0)^[21]时，纯体去相干γ_ν即为观测相干性。同时，假设森林区域电磁波的消光系数σ趋近于零，此时纯体去相干仅为森林高度h_v的函数，称为SINC函数模型^[17]，表达式为

 (4)

至此，当使用单基线观测数据时，纯体去相干γ_ν仅为树高h_v的函数关系，故可舍弃纯体去相干γ_ν的相位信息，并代入HoA=2π/k_z，其中HoA为测高模糊度(height of ambiguity)，表征InSAR信号的测高敏感度。可得SINC函数模型表达式如下^[17, 21]

 (5)

式(5)可视为y=SINC(x)，即y=sin(x)/x。

1.2 少量LiDAR数据辅助的TanDEM-X InSAR相干性反演森林高度方法

由上述分析可知，简化为SINC函数模型后，单基线TanDEM-X InSAR相干性仅为树高的函数关系，为了补偿模型忽略地体幅度比及消光系数引起的树高估计偏差，已有研究提出C-SINC模型，其表达式为^{[19-20, 22-23]}

 (6)

式中，C₁和C₂是引入的半经验参数，起到调节相干性与树高之间的关系，使得二者更满足函数关系，将相干性与森林高度在像素级别上进行画图，其示意图如图 2所示。图 2中不同颜色代表了点云的密度，颜色由蓝变黄再到红色，说明该平面坐标系内该位置的像素变多。曲线包围部分为参与估计参数的部分，SINC函数模型用红色虚线表示，C-SINC模型用蓝色虚线表示。另外，为避免不同基线引起的差异，在估计参数时采用h_v/HoA作为自变量进行归一化。当配置了合适的半经验参数时，C₁和C₂具有明确的物理意义。首先，由式(5)可知，当裸露的地面上不存在树高时，InSAR观测相干性应为1。然而在包括系统去相干因子γ_system等未补偿的残余去相干影响下|γ_v|小于1，因此增加半经验参数C₁补偿残余失相干因子。C₁值的大小在图 2纵轴上表现为由相干性为1处向下移动，一般位于0.9~1之间。C₂被认为进一步调节相干性与树高的关系以贴合点云密度较大的趋势，使得反演结果误差更小。如使用较长或较短的基线获取观测数据时，因测高敏感度存在差异而导致SINC函数模型无法准确描述相干性与树高之间的关系，借助C₂可纠正不同的基线配置、森林类型或随机噪声等非主导因素导致的模型匹配误差。C₂值的大小在图 2横轴上表现为在h_v/HoA为1处左右浮动，一般位于0.8~1.5之间。此时，借助C₁和C₂的共同调节，相较于SINC函数模型，C-SINC模型能较好地描述两者的关系。

图 2 C-SINC模型与SINC函数模型

Fig. 2 C-SINC model and SINC function model

图选项 

为确定合适的半经验参数C₁和C₂，选取小区域LiDAR树高数据作为树高真值辅助进行半经验参数的估计。需注意的是，虽然对有效垂直波数k_z进行局部校正，然而SAR系统的电磁波信号在坡度起伏区域的传输过程较为复杂，为使估计参数具有稳健性，估计时需排除较大坡度的影响，本文采取的准则为：视0坡度为中心，取正负坡度按等长度进行探测，设置50%为剩余像素百分比阈值，并将此区间外坡度较大的像素掩膜。同时，鉴于低相干性区域存在的不确定性和未知残余误差，掩膜相干性低于0.3的像素，以提高估计参数的稳健性^[14, 16]。其次，为避免误差较大的像素对参数估计准确性的影响，需排除图 2中距点云分布密集区域较远的像素，剩余参与估计参数点云部分如图 2曲线包围部分。最后，以非线性最小二乘估计准则为约束，使得观测相干性γ_vi^obs与模型期望相干性γ_vi^mod之差的平方和最小，其准则为

 (7)

式中，C₁=0.8~1；C₂=0.8~1.5。

建立C-SINC模型后，认为在SAR观测数据的影像覆盖区域森林属性具有相似性，利用小范围LiDAR数据估计参数得到单基线TanDEM-X InSAR相干性森林高度反演准则为

 (8)

另外，实际情况下，森林大多位于地形起伏的区域，在此场景下电磁波在森林散射体中传播路径改变，从而需要改进模型以适应地形起伏区域的森林高度反演。借助于外部DEM产品，在计算整个区域的森林高度时，基于式(9)对入射角进行校正计算局部k_z^[25-26]

 (9)

式中，θ′=θ-α。

至此，总结本文的反演森林高度流程如图 3所示。

图 3 少量LiDAR数据辅助的TanDEM-X InSAR相干性森林高度反演流程

Fig. 3 Flowchart of forest height inversion with TanDEM-X InSAR coherence aided a subset of LiDAR data

图选项 

2 试验结果及分析

2.1 试验区及数据介绍

本文选取位于西班牙境内的两个试验区，试验区为典型的地中海气候，位置如图 4所示。两个试验区内森林高度相对较低，主要集中在3~20 m之间，其LiDAR树高分布如图 5所示。另外，森林覆盖区域地形起伏较大，主要山脉沿东西走向，其中试验区Ⅰ的地形条件较为复杂，海拔高度变化集中在500~2000 m之间；试验区Ⅱ相对平坦，海拔高度主要集中在1200~2000 m之间，试验区的坡度统计如图 6所示。

图 4 试验区位置

Fig. 4 Location of the study sites

图选项 

图 5 试验区内LiDAR森林高度分布直方图

Fig. 5 Histogram of LiDAR forest height in the study sites

图选项 

图 6 两个试验区的地形坡度统计

Fig. 6 Statistics of terrain slope of the study sites

图选项 

获取覆盖两个试验区的单基线TanDEM-X InSAR影像，每个试验区干涉对包含了配准的主从复数影像，用于估计得到相干性，观测数据基本信息见表 1。

表 1 两个试验区选取的TanDEM-XInSAR数据基本信息Tab. 1 Information of TanDEM-X InSAR data in two study sites

区域	获取时间	模糊高/m	入射角/(°)	极化方式	kz/(rad/m)	分辨率(Rg×Az)/m
试验区Ⅰ	2012-12-28	36.6	35	HH	0.18	1.36×2.01
试验区Ⅱ	2012-10-23	32.3	38	HH	0.19	1.36×2.21

表选项 

另外，西班牙林业部门在2009—2015年间获取了高精度机载LiDAR数据，处理得到与干涉影像对获取时间同期的高分辨率(10×10 m)的LiDAR树高产品和DEM产品，以验证本文TanDEM-X InSAR相干性反演森林高度的性能与精度。

2.2 试验结果与讨论

2.2.1 反演结果

选取上述两个试验区内小范围LiDAR树高区域，如图 7中红色矩形框。其中，两个试验区的单个TanDEM-X InSAR干涉对覆盖面积分别约为2009 km²和1970 km²，而LiDAR树高覆盖面积约为129 km²和115 km²，分别占InSAR观测数据覆盖面积的6.4%和5.8%。绘制相干性与森林高度的点云密度图，分析不同高度的森林与相干性之间的关系，如图 7所示。由图 7可看出，SINC函数模型未穿过点云密度最大的区域，无法准确描述相干性与树高间的关系。其主要原因为基于SINC函数模型中，假设电磁波信号穿透至森林冠层高度一半^[13, 17]，而当森林冠层高度较矮时，电磁波信号可能穿透至冠层高度一半以下位置，在实际森林场景中相干性与树高间的关系未能严格满足SINC函数模型。为修正以上偏差建立C-SINC模型，需先通过小范围的LiDAR森林高度信息估计半经验参数C₁和C₂，利用像素之间的对应关系进行点云密度图绘制，根据式(7)约束进行参数估计的结果如图 7所示。

图 7 两个试验区的半经验参数估计结果

Fig. 7 Results of fitting the semi-empirical parameters in two study sites

图选项 

由图 7可看出，具有半经验参数的C-SINC模型能更准确描述相干性与森林高度之间的关系。两个试验区中森林高度超过HoA/2的点较少，其点云主要集中在0.1~0.5倍HoA处，与试验区内的森林高度主要集中在3~20 m的实际场景对应。另外，每个试验区的散点分布存在显著差异，试验区Ⅱ的散点分布相较于试验区Ⅰ较离散，而两个试验区超过HoA/2部分的点云都较为发散，其主要原因是复杂地形条件下，随着森林高度的增加，电磁波与森林的相互作用变得更为复杂，相干性的不确定性增加。由两个试验区估计得到的C₁和C₂可知，C₁为0.90和0.93，C₂为1.01和1.14，都为1附近的值，符合模型的预期假设。利用估计得到的参数进行整个SAR影像覆盖区域的森林高度反演，对得到的结果中相干性低于0.3的像素和无森林覆盖的区域进行掩膜，得到的树高制图结果如图 8所示，其中红色矩形框区域使用于拟合参数C₁和C₂。

图 8 两个试验区的森林高度制图结果

Fig. 8 The results of forest height mapping in two study sites

图选项 

由图 8可看出，在视觉效果上两个试验区的反演结果和LiDAR树高的形态特征总体上吻合较好，但都出现了高植被区被低估的现象，主要原因为图 8(a)中的C-SINC模型曲线在高植被区偏离了点云集中区域，其点云分布在高于HoA/2部分较为发散，经过最小二乘非线性约束得到的参数不能较好地表征所有像素的映射关系。为对反演结果进行定量分析，利用已有LiDAR数据进行验证与统计，计算得到的森林高度与LiDAR森林高度在样地水平下的均方根误差(root mean square error, RMSE)和相关系数R²如图 9所示。由图 9可以看出，两个试验区的RMSE为2.34 m和1.74 m，R²为0.80和0.82，从统计散点图的分布和两个统计指标来看，试验区Ⅰ的反演结果较差，其结果在低植被区存在明显的低估；试验区Ⅱ的反演结果相对较好，散点统计图较为收敛，相关性也较高。这是由于两方面因素影响：①由图 5可以看出，试验区Ⅱ内高植被覆盖较试验区Ⅰ较多，在图 7(b)点云中高植被覆盖较多，而随着森林高度的增加，电磁波与森林的相互作用过程变得更为复杂，相干性的不确定性增加，因此在图 7中试验区Ⅱ的点云图相对离散；②试验区Ⅰ的地形起伏较复杂，而在选取小范围LiDAR数据估计半经验参数时，为确保半经验参数的稳健型，在参数估计时仅使用了坡度较小的像素，估计得到半经验参数建立C-SINC模型后，用于反演所有坡度下的森林高度。因此，在反演所有坡度下的森林区域时，试验区Ⅰ的结果统计更分散，试验区Ⅱ因坡度相对较小反演结果更好，相关性更高。

图 9 两个试验区的森林高度反演结果与LiDAR森林高度的验证

Fig. 9 Validation of the forest height inversion results vs. LiDAR forest height of two study sites

图选项 

采用相同的试验区且控制其余流程一致，与SINC函数模型的反演结果进行对比，结果见表 2。由表 2可知，相比于SINC函数模型，引入半经验参数的C-SINC模型更适合描述相干性与森林高度之间的关系，其反演结果R²和RMSE取得了显著提升。

表 2 本文方法结果与SINC函数模型反演结果对比Tab. 2 Comparison of the results of the proposed method and the SINC function model

试验区	反演模型	(半经验)参数		精度验证指标
C1	C2		R2	RMSE/m
Ⅰ	SINC	1.00	1.00		0.63	4.6
本文C-SINC	0.90	1.01		0.80	2.34
Ⅱ	SINC	1.00	1.00		0.72	3.97
本文C-SINC	0.93	1.14		0.82	1.74

表选项 

2.2.2 坡度对于森林高度反演的影响

基线配置和局部地形条件决定InSAR系统获取观测数据的干涉质量和测高能力^[24]：当基线较长时，测高敏感度较强而相干性损失较多；当基线较短时，相干性能保持较高水平而测高敏感度不足；另外，不同的入射角对地形起伏区域也会使垂直波数k_z产生较大偏差^[25]。为分析坡度引起的估计偏差，不同坡度引起的局部垂直波数k_z变化情况如图 10蓝色虚线所示(红色虚线表示坡度为0，黑色实线表示名义k_z的大小)。

图 10 局部垂直波数k_z随坡度的变化关系

Fig. 10 The relationship between the local vertical wave number and slope

图选项 

由外部DEM产品统计试验区坡度分布信息可知，两个试验区无坡度低于-30°的区域；同时，由图 10可看出，当试验区Ⅰ的坡度大于20°、试验区Ⅱ的坡度大于30°时，其局部k_z的绝对值已超过名义k_z的两倍，测高敏感度过度饱和无法提供准确的测量结果，故将这部分极端起伏区域掩膜^[27]。另外，由图 8可知，树高制图结果中部分区域存在明显低估，由以上分析初步推测可能是由于局部极端坡度造成的测高失准导致。为此，利用图 8中的红色矩形部分将坡度、LiDAR树高和估计森林高度进行对比，如图 11所示，其中椭圆标注部分为存在明显低估区域。

图 11 图 8红色矩形区域的坡度、LiDAR树高与估计树高

Fig. 11 The slope information, LiDAR forest height and estimated height of the area covered by red rectangle in fig. 8

图选项 

从视觉效果上可看出，在正坡度区域存在明显低估，其主要原因与局部k_z的测高敏感度饱和有关。已有研究展开关于甄选最优基线进行TanDEM-X InSAR数据的树高估计探索，建议当提供多余观测的基线数据时，可选择HoA为高植被高度的两倍基线作为合适观测数据用以反演森林高度，其在坡度较大区域仍能提供合适的测高能力以避免起伏区域造成的估计偏差^[24, 28]。

另外，充分利用已有外部辅助LiDAR树高数据(图 8中的红色矩形区域)，分析坡度与反演结果偏差的关系，并建立树高与估计偏差之间的回归函数关系(图 12)。由图 12可以看出，两个试验区在正坡区域存在明显的低估趋势；试验区Ⅰ在负坡区域存在离散的高估和低估现象，而试验区Ⅱ在负坡区域存在显著的高估趋势。

图 12 坡度与树高估计偏差的关系

Fig. 12 The relationship between slope and forest height estimation bias

图选项 

借助已建立的坡度与估计偏差之间的关系，对上述反演结果进行补偿，再次统计反演精度，结果如图 13所示。

图 13 两个试验区的森林高度反演结果与LiDAR森林高度的验证

Fig. 13 Validation of the forest height inversion results vs. LiDAR forest height of two study sites

图选项 

由图 13可以看出，当把坡度较大的区域掩膜和进行坡度估计偏差补偿后，两个试验区的精度都得到了明显的提高，RMSE分别提高0.23 m和0.45 m，R²分别提高0.04和0.08，散点分布相比所有坡度的统计结果较为收敛，尤其是试验区Ⅱ，正负坡度像素均能得到有效补偿，其精度提高较多；而试验区Ⅰ负坡度主导区域内的地形趋势，其估计偏差较为发散，无法准确削弱负坡度的影响。

总结当所有坡度参与精度统计时结果较差的原因有如下两点：其一，当估计半经验参数C₁和C₂时，为避免坡度对于估计参数准确性造成粗差影响，仅提取坡度较小的像素用于估计，对于坡度较大的像素适用性较差；其二，坡度对于SAR系统的侧视成像造成了相干性的损失，由此造成的失相干在模型中被参数C₁仅吸收一部分，但不能准确补偿。解决坡度使得反演森林高度造成误差较大的问题，可从两方面解决：一方面是提供多角度的观测数据进行联合解算，以解决极端坡度区域的估计偏差问题；另一方面可通过建模补偿坡度区域造成的非体去相干，这是在后续工作中亟须解决的问题。

3 结论

本文提出了少量外部数据辅助的单基线TanDEM-X InSAR相干性森林高度反演方法，该方法不依赖于全局高精度的DEM提供地表相位信息，通过两个不同场景下的试验区进行验证表明：

(1) 加入适应场景和基线配置的半经验参数后，C-SINC模型能够在复杂森林场景中具备较好的反演性能，能修正SINC函数模型中的模型假设和场景差异引入的误差，较平坦试验区能够达到1.29 m的精度，且具有较好相关性。

(2) 坡度是制约高精度反演森林高度的主要因素之一，通过分析坡度对反演结果的影响表明：通过小范围的LiDAR数据辅助估计得到参数能够抵抗较小坡度的影响，但在坡度起伏较大的区域反演性能下降，通过外部数据的回归分析可进一步减小由坡度引起的估计偏差，下一步的工作将围绕深入探索坡度引起的去相干问题展开研究。

(3) 本文方法的生产实际意义在于：当使用机载LiDAR获取局部小范围的森林高度时，可辅助单基线TanDEM-X InSAR数据进行大范围森林高度反演。具体地，可反演LiDAR数据约20倍面积的单个TanDEM-X InSAR干涉对面积的森林高度。同时，当获取到同时期的相同基线配置相同区域的TanDEM-X InSAR观测数据时，所估计的半经验参数可继续辅助相似干涉对反演大范围的森林高度。

致谢

感谢德国宇航局(DLR)为本文提供的TanDEM-X InSAR观测数据(OTHER7349)；感谢西班牙林业部门和阿利坎特大学(University of Alicante)提供的地面实测林业验证数据。

作者简介

第一作者简介：张涛(1996—), 男, 博士生, 研究方向为PolInSAR数据处理与森林参数反演。E-mail: [e-mail]

通信作者：朱建军, E-mail: [e-mail]

初审：张艳玲

复审：宋启凡
终审：金   君

往期推荐

资讯