CN112363165B - 一种林下地形反演方法、装置、设备及介质 - Google Patents

Abstract

本发明公开了一种林下地形反演方法、装置、设备及介质，方法为：获取林下地形反演区域的双站TanDEM‑X数据并进行干涉处理，获取干涉相位和体失相干系数；对干涉相位去除平地、模拟地形、噪声、轨道误差等相位，得到残余地形相位；使用残余地形相位和模拟地形相位重建高精度InSAR地形相位；构建林下反演模型，联合ICESat‑2地面测高点进行平差计算，求解模型中的未知参数；使用体失相干系数和高精度InSAR地形相位信息，根据林下反演模型计算整个反演区域的林下地形数据。本发明是一种稳健可行的林下地形估算方法，可以广泛应用于TanDEM‑X数据测量大范围的林下地形高度。

Description

一种林下地形反演方法、装置、设备及介质

技术领域

本发明涉及基于合成孔径雷达干涉测量(InSAR)的大地测量领域，具体是指一种基于 双站TanDEM-X数据的林下地形反演方法、装置、设备及介质。

背景技术

数字高程模型(digital elevation model,DEM)是国家经济建设、社会发展和国家安全必 不可少的基础图件。尤其，高精度的覆盖层下地形高度信息(如：林下地形)对于自然资源 调查、基础设施建设、国家经济发展以及国防安全等领域均具有十分重要的意义。目前，为 了获取高精度的地形高度信息，合成孔径雷达(InSAR)、LiDAR技术以及光学遥感技术是目 前常用的技术手段。其中，InSAR技术具备全天候、全天时、观测周期短、大范围的优势， 成为了一种获取全球DEM的强有力工具，如：SRTM DEM和TanDEM-X DEM等全球或近全球覆盖的DEM均是采用InSAR技术所获取。此外，TanDEM-X采用了TerraSAR-X和 TanDEM-X组成的单轨、双星串行星座，这一特殊的数据获取模式最大程度上减少了时间去 相关、大气延迟以及地表型形变等因素多获取高精度地形的影响。基于上述优势，TanDEM-X DEM有望成为覆盖范围最广、精度最高的全球DEM数据，代表了全球DEM产品的发展方 向。然而，尽管微波信号具备一定的穿透能力，但是在森林区，短波长(X波段)的微波信 号难以直接“穿透”植被覆盖层直接获取林下地形。为了获取大范围、高精度的林下地形高 度，必须对干涉相位中包含的植被高度信号进行有效的去除。

目前，常用的林下地形获取方法主要包括以下两种：1)基于外部辅助数据校正现有DSM 产品中的树高信号；2)从极化数据中分离地表散射的贡献。其中，第一类方法通常需要依赖 大量的外部辅助数据(如：温度数据、降水数据、低分辨率的植被高度数据)进行建模，且 为了获取高精度的林下地形结果，通常需要具备研究区域一定的先验信息，并且辅助数据的 精度和模型的选择也进一步限制了该方法的实用性；相较于第一类方法，第二类方法要求多 极化数据进行地表散射信息的分离，由于目前尚缺少可获取多极化数据的星载SAR系统，目 前的方法只能通过机载SAR数据进行相关的实验，这就限制了大范围林下地形的获取。因此， 有必要设计一种基于星载SAR数据的大范围、高精度林下地形获取方法。

发明内容

本发明所要解决的技术问题在于，提供一种基于双站TanDEM-X数据的林下地形反演方 法、装置、设备及介质，可以实现快速、大范围、高精度的林下地形反演。

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

一种基于双站TanDEM-X数据的林下地形反演方法，包括以下步骤：

步骤1，获取林下地形反演区域的双站TanDEM-X数据并进行干涉处理，获取干涉相位 和相干系数，并对相干系数进行校正获取相应的体失相干系数；

步骤2，对步骤1得到的干涉相位，首先进行去平地相位处理和使用外部DEM数据进行 去地形处理得到差分干涉相位，然后对差分干涉相位依次进行滤波和解缠处理得到解缠差分 干涉相位，再从解缠差分相位去除轨道误差相位，最终得到残余地形相位；

步骤3，使用外部DEM数据模拟地形相位，使用模拟地形相位和残余地形相位重建高精 度InSAR地形相位；

步骤4，根据林下地形高度、植被体散射相位中心与InSAR地形相位之间的关系，构建 林下地形反演模型：

式中，φ_topo为真实地表相位，γ_vol为体失相干系数，λ为雷达信号波长，R表示雷达到目 标的斜距，θ表示主传感器对应的侧视角，B_⊥代表垂直基线长度，k_z表示垂直向有效波数； s和m均为模型的待估参数，分别表示植被高度与穿透深度之间的线性系数和误差系数；

步骤5，从星载LiDAR数据中提取反演区域的若干地面控制点，获取每个地面控制点的 林下地形高度、InSAR地形相位、对应的体失相干系数，并采用最小二乘法计算林下地形反 演模型中的待估参数s和m；

步骤6，根据步骤1得到的体失相干系数和步骤3得到的高精度InSAR地形相位，使用 估算参数s和m的林下地形反演模型，计算整个反演区域的林下地形数据。

在更优的技术方案中，采用稳健最小二乘法计算林下地形反演模型中的待估参数s和m， 具体为：

步骤5.1，使用n个若干地面控制点的林下地形高度、InSAR地形相位、对应的体失相干 系数代入到林下地形反演模型，得到以下矩阵表达式：

式中，h_i、φ_{topo_i}、

分别表示星载LiDAR数据中第i个地面控制点的林下地形高度、 InSAR地形相位、体失相干系数，i＝1,2,…,n；

步骤5.2，采用稳健最小二乘估计对未知参数X进行估计：

X＝(B^TPB)^-1(B^TPL)，

其中，

式中，X为待估参数向量，B为系数矩阵，L为观测向量，P为权阵。通过三段法的IGGIII方案等价权函数确定每个地面控制点在每次迭代的权，计算表达式为：

其中，

代表第i个地面控制点经迭代后得到的数据标准化残差；

表示经过迭代更新后 第i个地面控制点的权；k₀,k₁代表两个常量因子，标准化残差表达式如下：

式中，σ₀为单位权中误差；q_vi为权倒数；

步骤5.3，由所有n个地面控制点在每次迭代的权，得到其权阵为：

式中，c为迭代次数，从而可得每次迭代的X和相应的标准化残差V分别为：

X^c＝(B^TP^cB)^-1(B^TP^cL)，

V^c＝BX^c-L；

步骤5.4，重复步骤5.2、步骤5.3，直到满足迭代运算的终止条件为|X^c+1-X^c|<ε，即得 到X中的待估参数s和m。

在更优的技术方案中，步骤1得到的干涉相位φ包括真实地形相位φ_topo、平地相位φ_flat、 轨道误差相位φ_orbit和噪声相位φ_noi，表示为：

φ＝φ_topo+φ_flat+φ_oribit+φ_noi；

使用外部DEM数据进行去地形处理得到差分干涉相位具体为：使用外部DEM数据模拟 地形相位φ_{sim_topo}，然后将得到的模拟地形相位φ_{sim_topo}从干涉相位φ中扣除，得到差分干涉相 位φ_diff：

φ_diff＝φ_{res_topo}+φ_orbit+φ_noi；

式中，h是从外部DEM中获取的地面点高度，φ_{res_topo}为外部DEM数据模拟地形相位后 剩余的残余地形相位。

在更优的技术方案中，去除平地相位采用的平地相位计算公式为：

式中，φ_flat为平地相位，R₁和R₂分别表示主、从雷达天线中心到地面目标点的斜距。

在更优的技术方案中，轨道误差相位的去除方法为：使用融入地表高程h的多项式模型 拟合轨道误差相位，通过最小二乘法解算多项式模型中的未知参数，然后将轨道误差相位从 解缠差分干涉相位去除。

在更优的技术方案中，采用SNR校正方法对相干系数进行失相关补偿，得到体失相干系 数。

一种基于双站TanDEM-X数据的林下地形反演装置，包括：

数据预处理模块，用于：获取林下地形反演区域的双站TanDEM-X数据并进行干涉处理， 获取干涉相位和相干系数，并对相干系数进行校正获取相应的体失相干系数；

残余地形相位提取模块，用于：对数据预处理模块得到的干涉相位，首先进行去平地相 位处理和使用外部DEM数据进行去地形处理得到差分干涉相位，然后对差分干涉相位依次 进行滤波和解缠处理得到解缠差分干涉相位，再从解缠差分相位去除轨道误差相位，最终得 到残余地形相位；

高精度地形相位重构模块，用于：使用外部DEM数据模拟地形相位，使用模拟地形相 位和残余地形相位重建高精度InSAR地形相位；

林下地形反演模型构建模块，用于：根据林下地形高度、植被体散射相位中心与InSAR 地形相位的关系，构建林下地形反演模型：

式中，φ_topo为真实地表相位，γ_vol为体失相干系数，λ为雷达信号波长，R表示雷达到目 标的斜距，θ表示主传感器对应的侧视角，B_⊥代表垂直基线长度，k_z表示垂直向有效波数； s和m均为模型的待估参数，分别表示植被高度与穿透深度之间的线性系数和误差系数；

待估参数计算模块，用于：从星载LiDAR数据中提取反演区域的若干地面控制点，获取 每个地面控制点的林下地形高度、InSAR地形相位、对应的体失相干系数，并采用最小二乘 法计算林下地形反演模型中的待估参数s和m；

林下地形计算模块，用于：根据数据预处理模块得到的体失相干系数和高精度地形相位 重构模块得到的高精度InSAR地形相位，使用已知参数s和m的林下地形反演模型，计算整 个反演区域的林下地形数据。

一种设备，包括处理器和存储器；其中：所述存储器用于存储计算机指令；所述处理器 用于执行所述存储器存储的计算机指令，具体执行上述任一技术方案所述的方法。

一种计算机存储介质，用于存储程序，所述程序被执行时，用于实现上述任一技术方案 所述的方法。

有益效果

本发明的基于双站TanDEM-X数据的林下地形提取方法，该方法将InSAR数据与LiDAR 数据进行有机融合：首先利用常用常规的InSAR数据处理方法，通过基于外部DEM数据的 InSAR处理技术重建高精度的InSAR地形相位；随后利用简化的RVoG模型(即SINC模型) 并结合得到的体失相干系数对植被引起的散射相位中心高度进行建模，构建本发明提出的林 下地形反演模型；利用部分ICESat-2提供的地表高度数据作为林下地形参考数据，采用稳健 最小二乘法估算上述模型中的未知参数；最后，根据得到的体失相干系数和重建的InSAR地 形相位，计算整个区域的林下地形。该方法根据体失相干系数和SINC模型对InSAR散射相 位中心高度进行建模，构建了基于双站TanDEM-X数据的林下地形反演模型，突破了传统 InSAR技术无法去除干涉相位中的植被高度的限制；并通过易于获取的星载LiDAR数据作为 参考数据进行模型的计算，实现了InSAR与LiDAR数据的深度融合。其原理更为直观，且 易于编程实现和广泛应用，是一种可行的林下地形获取方法，可以实现快速、大范围、高精 度的林下地形反演。

附图说明

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

图2为本发明所选择的验证区域位置及SAR影像的覆盖范围。

图3为本发明残余地形相位提取模块涉及的差分干涉相位示意图。其中，图3(a)代表初 始的差分干涉相位，图3(b)代表经轨道误差去除后的差分干涉相位结果。

图4为生成的InSAR DEM(见图4(c))及采用本发明的方法与步骤得到的林下地形结果 (见图4(b))，以及上述结果与参考LiDAR DTM(图4(a))之间的误差图，其中图4(d)为图 4(b)与图4(a)之间的误差图，图4(e)为图4(c)与图4(a)之间的误差图。

具体实施方式

为了更好的说明本发明的方法与步骤，利用位于瑞典北部BioSAR2008试验区的双站 TanDEM-X雷达数据，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施仅 用以解释本发明，并不用于限定本发明。

本实施例提供一种基于双站TanDEM-X数据的林下地形反演方法，参照图1所示，包括 以下步骤：

步骤1，根据所选的林下地形反演研究区域的坐标范围，获取德国宇航中心对应的双站 TanDEM-X数据，对获取的双站TanDEM-X数据(DLR已经通过ITP数据处理软件对数据进行了一定的辐射校正和精配准)进行干涉处理，获取干涉相位和相干系数。

本实施例选择了位于瑞典北部BioSAR2008试验区的TanDEM-X数据，其数据获取时间 为2011年12月11日。试验区位置以及SAR影像的覆盖范围如图2所示。首先，对经过配 准后的数据进行共轭相乘，获得相应的干涉条纹图。

相干性估计是InSAR测量中的一项重要工作，相干系数的大小代表了InSAR相位质量的 好与坏。为了避免滤波对相干性估计的干扰，我们需要在对干涉图滤波之前采用窗口估计方 法获取相干系数，其公式为：

式中，M和N为计算相干性的数据块尺寸大小；m和n为数据块内行列号；μ₁(m,n)，μ₂(m,n) 为主、副影像数据块内影像坐标(m,n)处的复数值；|·|²为数据的二阶范数。

在得到相干系数后，本实施例采用SNR校正方法对相干系数进行失相关补偿，获取相应 的体失相干系数。SNR校正方法为现有技术，可参考文献Kugler F.,Schulze D.,Hajnsek I.,et al. TanDEM-X Pol-InSAR Performance for Forest Height Estimation[J].IEEE Transactions on Geoence&Remote Sensing,2014,52(10):6404-6422.，本实施例不再详细阐述。

步骤2，对步骤1得到的干涉相位，首先进行去平地相位处理和使用外部DEM数据进行 去地形处理得到差分干涉相位，，然后对差分干涉相位依次进行滤波和解缠处理得到解缠差分 干涉相位，再从解缠差分相位去除轨道误差相位，最终得到残余地形相位。

由于双站TanDEM-X数据的主、从两幅影像几乎是同时获得，因此可以排除大气和形变 等因素的影响，因此，通过常规干涉处理得到的上述干涉相位φ可表示为如下几个部分：

φ＝φ_topo+φ_flat+φ_oribit+φ_noi

式中，φ_topo为代表与地形测量有关的真实地形相位，包含“裸地球”高度和地面目标高 度(如：植被高度信号)；φ_flat为平地相位，φ_orbit为轨道误差相位，φ_noi为噪声相位。

因此，为了获取高精度的林下地形，只需要保留与“裸地球”高度相关的相位成分，而 其余的相位成分需要进行有效的去除或削弱。

首先，需要根据雷达成像的轨道状态矢量参数和成像几何关系计算平地相位，并将其从 干涉相位φ中去除。其中，平地相位的计算公式如下：

式中，B_⊥代表垂直基线长度，R₁和R₂分别表示主、从雷达天线中心到地面目标点的斜 距，λ是雷达信号的波长，θ表示主传感器对应的侧视角。

然后，使用外部DEM数据进行去地形处理得到差分干涉相位，具体为：使用外部DEM数据模拟地形相位φ_{sim_topo}，再将得到的模拟地形相位φ_{sim_topo}从干涉相位φ中扣除，得到差分 干涉相位φ_diff：

φ_diff＝φ_{res_topo}+φ_orbit+φ_noi；

式中，h是从外部DEM中获取的地面点高度，φ_{res_topo}为外部DEM数据模拟地形相位后 剩余的残余地形相位。

为了抑制噪声对干涉相位的影响，在对差分干涉相位φ_diff解缠之前需要对相位进行相应 的滤波处理，去除噪声相位。通常采用的滤波方式包括Glodstein滤波以及Non-Local(NL) 滤波。随后，就可以对滤波后的差分干涉相位进行解缠处理，得到解缠后的差分干涉相位。

由于星载SAR系统的轨道误差在雷达坐标系下通常表现为规则的条纹，如图3(a)所示， 因此可以采用融入高程信息的多项式曲面模型对上述得到的解缠差分干涉图进行拟合处理去 除所包含的轨道相位。其中，采用的多项式曲面模型如下式所示：

φ_orbit＝a₀+a₁x+a₂y+a₃xy+a₄h

式中，x,y分别表示雷达坐标系下像素点方位向和距离向下的坐标，a_i(i＝1,2,3,4)表示 模型中的待估参数，h代表像素点对应的地面高度。为了准确估计模型中的参数，需要采用 最小二乘的方法。

经过上述步骤2中所述的相位处理，最终得到残余地形相位φ_{diff_unw}。轨道误差相位去除 后的差分干涉相位如图3(b)所示。

步骤3，使用外部DEM数据模拟地形相位φ_{sim_topo}，使用模拟地形相位φ_{sim_topo}和残余地形 相位φ_{diff_unw}重建高精度InSAR地形相位：φ_topo＝φ_{sim_topo}+φ_{diff_unw}。

步骤4，根据林下地形高度、植被体散射相位中心与InSAR地形相位的关系，构建林下 地形反演模型：

式中，φ_topo为真实地表相位，γ_vol为体失相干系数，λ为雷达信号波长，R表示雷达到目 标的斜距，θ表示主传感器对应的侧视角，B_⊥代表垂直基线长度，k_z表示垂直向有效波数； s和m均为模型的待估参数，分别表示植被高度与穿透深度之间的线性系数和误差系数。

在森林区，由于地表植被的覆盖，上述得到的地形相位包含两部分：“裸地球”高度(即 林下地形高度)和部分植被的高度。因此，根据雷达成像几何可以将真实InSAR地形相位φ_topo表示为：

式中，h_DTM表示林下地形高度，h_forest代表体散射导致的高度，即散射相位中心的高度 (scatter phase center height:SPCH，代表散射点与地面间的高度差)。根据上式，林下地形高 度可以通过下式计算：

为了通过上式获取林下地形高度，我们需要已知体散射贡献所导致的相位中心高度 h_forest。而在通常情况下，仅通过干涉相位无法直接区分地形高度和散射相位中心的高度。但 是，森林的体散射所导致的体失相干γ_vol不仅提升了干涉相位中心的高度，也会导致相干性的 降低。在已有研究中，通过引入简化的RVoG(随机地体两层散射)模型，如：SINC模型， 体失相干已经被广泛应用于森林高度的估算，其公式如下：

式中，h_v代表植被高度；γ_vol代表体失相干系数，可通过对相干系数进行失相关补偿得 到(对TanDEM-X而言，主要是SNR校正)，k_z是垂直波数，并可以通过下式表示：

对于TanDEM-X而言，受限于X-波段的弱穿透性，所得到的通常是信号在森林区的穿透 深度(penetration depth:PD，表示树顶端值至顶端与地面间某一位置的高度)。且常用的SINC 模型指出，植被高度在数值上等于2倍的穿透深度，即

因此， 通过移除去除散射相位中心的高度，可以考虑采用如下公式获取林下地形，即：

该模型在不依赖外部数据的情况下，可以被用来估算林下地形。但是，采用SINC模型 无法被用来准确的描述模式复杂的森林场景，即植被高度等于2倍的穿透深度不能适用于所 有的森林场景。为了补偿上述模型在林下地形估算中存在的模型误差，本实施例引入部分ICESat-2地面高度点来解决。其核心思想是：已有研究表明，植被高度与穿透深度是线性相 关的，因此散射相位中心的高度可以通过与穿透深度相关的一个线性模型来表示。最终可得 本实施例建立的林下地形反演模型如下：

其中，s和m代表模型的待估参数，分别表示植被高度与穿透深度之间的线性系数和误 差系数，可以在一定程度上修正模型的匹配误差。

步骤5，从星载LiDAR数据中提取反演区域的若干地面控制点，获取每个地面控制点的 林下地形高度、InSAR地形相位、对应的体失相干系数，并采用最小二乘法计算林下地形反 演模型中的待估参数s和m。

根据上述建立的林下地形反演模型，本实施例引入一些地面控制点来求解模型中的待估 参数s和m，而近年来公开发布的星载LiDAR数据(如：ICESat-2、GEDI)为这一问题的解 决提供了契机。本实施例基于ICESat-2所包含的林下地形高度信息，采用最小二乘平差来估 算待估参数s和m，具体包括：

步骤5.1，使用n个若干地面控制点的林下地形高度、InSAR地形相位、对应的体失相干 系数代入到林下地形反演模型，得到以下矩阵表达式：

式中，h_i、φ_{topo_i}、

分别表示星载LiDAR数据中第i个地面控制点的林下地形高度、 InSAR地形相位、体失相干系数，i＝1,2,...,n；

步骤5.2，为了降低残余误差相位对最小二乘估计的干扰，采用稳健最小二乘估计对未知 参数X进行估计：

X＝(B^TPB)^-1(B^TPL)，

式中，X为待估参数向量，B为系数矩阵，L为观测向量，P为权阵。通过三段法的IGGIII方案等价权函数确定每个地面控制点在每次迭代的权，计算表达式为：

其中，

代表第i个地面控制点经迭代后得到的数据标准化残差；

表示经过迭代更新后 第i个地面控制点的权；k₀,k₁代表两个常量因子，且通常的经验取值为k₀＝1.5，k₁＝2.5。其 中，标准化残差表达式如下：

式中，σ₀为单位权中误差；q_vi为权倒数；

步骤5.3，由所有n个地面控制点在每次迭代的权，得到其权阵为：

式中，c为迭代次数，从而可得每次迭代的X和相应的标准化残差V分别为：

X^c＝(B^TP^cB)^-1(B^TP^cL)，

V^c＝BX^c-L；

步骤5.4，重复步骤5.2、步骤5.3，直到满足迭代运算的终止条件为|X^c+1-X^c|<ε，即得 到X中的待估参数s和m；其中，ε为一极其小的预设量值，一般设置为10^-6。

步骤6，根据经SNR校正得到的体失相干系数和步骤3得到的高精度InSAR地形相位， 使用已知参数s和m的林下地形反演模型，计算整个反演区域的林下地形数据。

基于上述分析，本实施例从NASA官网下载了覆盖该研究区的ICESat-2数据作为外部参 考数据，并采用python编写了读取和处理ICESat-2原始“h5”文件的相关程序。最后，将提 取出的ICESat-2数据作为地面控制点，并联合得到的体失相干系数和重建的高精度InSAR地 形相位进行建模，并采用稳健最小二乘估计获取相应的参数。随后便可获取整个研究区域的 林下地形数据，如图4(b)所示。为了验证所提方案的正确性和有效性，同时也采用得到的 地形相位直接计算得到了传统InSAR对应的地面高度，如图4(c)所示。最后，选择覆盖该 研究区的机载LiDAR地面数据(图4(a))对上述结果的精度进行了相关的验证，其误差分 布如图4(d)和(e)所示。可以看到，经过所提方法得到的DTM更接近与地面验证数据， 说明该方法可以有效的去除传统InSAR结果中的植被高信息，因此，可以应用于大范围林下 地形的获取。此外，随着ICESat-2地面点数据量的不断增加，为大范围、高精度林下地形提 取提供了可能。

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

Claims (9)

1.一种基于双站TanDEM-X数据的林下地形反演方法，其特征在于，包括以下步骤：

步骤1，获取林下地形反演区域的双站TanDEM-X数据并进行干涉处理，获取干涉相位和相干系数，并对相干系数进行校正获取相应的体失相干系数；

步骤2，对步骤1得到的干涉相位，首先进行去平地相位处理和使用外部DEM数据进行去地形处理得到差分干涉相位，然后对差分干涉相位依次进行滤波和解缠处理得到解缠差分干涉相位，再从解缠差分相位去除轨道误差相位，最终得到残余地形相位；

步骤3，使用外部DEM数据模拟地形相位，使用模拟地形相位和残余地形相位重建高精度InSAR地形相位；

步骤4，根据林下地形高度、植被体散射相位中心与InSAR地形相位之间的关系，构建林下地形反演模型：

式中，φ_topo为真实地表相位，γ_vol为体失相干系数，λ为雷达信号波长，R表示雷达到目标的斜距，θ表示主传感器对应的侧视角，B_⊥代表垂直基线长度，k_z表示垂直向有效波数；s和m均为模型的待估参数，分别表示植被高度与穿透深度之间的线性系数和误差系数；

步骤5，从星载LiDAR数据中提取反演区域的若干地面控制点，获取每个地面控制点的林下地形高度、InSAR地形相位、对应的体失相干系数，并采用最小二乘法计算林下地形反演模型中的待估参数s和m；

步骤6，根据步骤1得到的体失相干系数和步骤3得到的高精度InSAR地形相位，使用估算参数s和m的林下地形反演模型，计算整个反演区域的林下地形数据。

2.根据权利要求1所述的方法，其特征在于，采用稳健最小二乘法计算林下地形反演模型中的待估参数s和m，具体为：

步骤5.1，使用n个若干地面控制点的林下地形高度、InSAR地形相位、对应的体失相干系数代入到林下地形反演模型，得到以下矩阵表达式：

式中，h_i、φ_{topo_i}、

分别表示星载LiDAR数据中第i个地面控制点的林下地形高度、InSAR地形相位、体失相干系数，i＝1,2,…,n；

步骤5.2，采用稳健最小二乘估计对未知参数X进行估计：

X＝(B^TPB)^-1(B^TPL)，

其中，

X＝[s m]^T；

式中，X为待估参数向量，B为系数矩阵，L为观测向量，P为权阵，通过三段法的IGG III方案等价权函数确定每个地面控制点在每次迭代的权，计算表达式为：

其中，

代表第i个地面控制点经迭代后得到的数据标准化残差；

表示经过迭代更新后第i个地面控制点的权；k₀,k₁代表两个常量因子，标准化残差表达式如下：

式中，σ₀为单位权中误差；q_vi为权倒数；

步骤5.3，由所有n个地面控制点在每次迭代的权，得到其权阵为：

式中，c为迭代次数，从而可得每次迭代的X和相应的标准化残差V分别为：

X^c＝(B^TP^cB)^-1(B^TP^cL)，

V^c＝BX^c-L；

步骤5.4，重复步骤5.2、步骤5.3，直到满足迭代运算的终止条件为|X^c+1-X^c|<ε，即得到X中的待估参数s和m。

3.根据权利要求1所述的方法，其特征在于，步骤1得到的干涉相位φ包括真实地形相位φ_topo、平地相位φ_flat、轨道误差相位φ_orbit和噪声相位φ_noi，表示为：

φ＝φ_topo+φ_flat+φ_oribit+φ_noi；

使用外部DEM数据进行去地形处理得到差分干涉相位具体为：使用外部DEM数据模拟地形相位φ_{sim_topo}，然后将得到的模拟地形相位φ_{sim_topo}从干涉相位φ中扣除，得到差分干涉相位φ_diff：

φ_diff＝φ_{res_topo}+φ_orbit+φ_noi；

式中，h是从外部DEM中获取的地面点高度，φ_{res_topo}为外部DEM数据模拟地形相位后剩余的残余地形相位。

4.根据权利要求1所述的方法，其特征在于，去除平地相位采用的平地相位计算公式为：

式中，φ_flat为平地相位，R₁和R₂分别表示主、从雷达天线中心到地面目标点的斜距。

5.根据权利要求1所述的方法，其特征在于，轨道误差相位的去除方法为：使用融入地表高程h的多项式模型拟合轨道误差相位，通过最小二乘法解算多项式模型中的未知参数，然后将轨道误差相位从解缠差分干涉相位去除。

6.根据权利要求1所述的方法，其特征在于，采用SNR校正方法对相干系数进行失相关补偿，得到体失相干系数。

7.一种基于双站TanDEM-X数据的林下地形反演装置，其特征在于，包括：

数据预处理模块，用于：获取林下地形反演区域的双站TanDEM-X数据并进行干涉处理，获取干涉相位和相干系数，并对相干系数进行校正获取相应的体失相干系数；

残余地形相位提取模块，用于：对数据预处理模块得到的干涉相位，首先进行去平地相位处理和使用外部DEM数据进行去地形处理得到差分干涉相位，然后对差分干涉相位依次进行滤波和解缠处理得到解缠差分干涉相位，再从解缠差分相位去除轨道误差相位，最终得到残余地形相位；

高精度地形相位重构模块，用于：使用外部DEM数据模拟地形相位，使用模拟地形相位和残余地形相位重建高精度InSAR地形相位；

林下地形反演模型构建模块，用于：根据林下地形高度、植被体散射相位中心与InSAR地形相位的关系，构建林下地形反演模型：

式中，φ_topo为真实地表相位，γ_vol为体失相干系数，λ为雷达信号波长，R表示雷达到目标的斜距，θ表示主传感器对应的侧视角，B_⊥代表垂直基线长度，k_z表示垂直向有效波数；s和m均为模型的待估参数，分别表示植被高度与穿透深度之间的线性系数和误差系数；

待估参数计算模块，用于：从星载LiDAR数据中提取反演区域的若干地面控制点，获取每个地面控制点的林下地形高度、InSAR地形相位、对应的体失相干系数，并采用最小二乘法计算林下地形反演模型中的待估参数s和m；

林下地形计算模块，用于：根据数据预处理模块得到的体失相干系数和高精度地形相位重构模块得到的高精度InSAR地形相位，使用已知参数s和m的林下地形反演模型，计算整个反演区域的林下地形数据。

8.一种设备，其特征在于，包括处理器和存储器；其中：所述存储器用于存储计算机指令；所述处理器用于执行所述存储器存储的计算机指令，具体执行如权利要求1-6任一所述的方法。

9.一种计算机存储介质，其特征在于，用于存储程序，所述程序被执行时，用于实现如权利要求1-6任一所述的方法。

Images