CN113341410B - 一种大范围林下地形估计方法、装置、设备及介质 - Google Patents

一种大范围林下地形估计方法、装置、设备及介质 Download PDF

Info

Publication number
CN113341410B
CN113341410B CN202010771942.5A CN202010771942A CN113341410B CN 113341410 B CN113341410 B CN 113341410B CN 202010771942 A CN202010771942 A CN 202010771942A CN 113341410 B CN113341410 B CN 113341410B
Authority
CN
China
Prior art keywords
phase
forest
terrain
height
interference
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010771942.5A
Other languages
English (en)
Other versions
CN113341410A (zh
Inventor
朱建军
王会强
付海强
李志伟
胡俊
刘志卫
张兵
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202010771942.5A priority Critical patent/CN113341410B/zh
Publication of CN113341410A publication Critical patent/CN113341410A/zh
Application granted granted Critical
Publication of CN113341410B publication Critical patent/CN113341410B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques

Abstract

本发明公开了一种大范围林下地形估计方法、装置、设备及介质，其方法为：步骤1，获取TanDEM‑X双站SLC数据进行处理得到干涉相位和相干幅度；步骤2，使用外部DEM数据计算地形残差，并对外部DEM数据进行补偿，得到TanDEM‑X双站干涉条件下的高精度的InSAR DEM；步骤3，根据干涉幅度计算森林研究区域全场景的穿透深度d，然后利用引入的稀疏分布的先验点森林高度信息，拟合有关于穿透深度与森林高度之间的线性关系，从而计算研究区域全场景的森林高度以及散射相位中心高度；步骤4，从步骤2得到的TanDEM‑X双站的InSAR DEM中扣除步骤3得到的散射相位中心高度，即实现大范围林下地形重建。

Description

一种大范围林下地形估计方法、装置、设备及介质
技术领域
本发明涉及合成孔径雷达干涉测量林下地形技术领域，具体涉及一种基于TanDEM-X双站单基线干涉数据的大范围林下地形估计方法、装置、设备及介质。
背景技术
近些年来，合成孔径雷达干涉测量(InSAR)技术已经具备了大范围高精度地形测绘的能力，比如目前已经发布了覆盖全球的高精度DEM数据，如STRM DEM、TanDEM-X DEM等。特别是TanDEM-X双站干涉模式，由于不受时间失相关和大气延迟误差的影响，TanDEM-XDEM可达到约为10米的全球绝对高程精度，被公认为是对地观测精度最高且最连续的全球DEM数据源之一。然而，为了实现高质量的气候变化监测、地质灾害监测和自然资源调查等目的，即便TanDEM-X DEM也是达不到要求，特别是在自然媒介覆盖区，如冰雪、沙漠和植被等。以植被覆盖为例，InSAR生产的数字高程模型(Digital Elevation Model，简称DEM)表征的是介于植被冠层顶部和层下地表之间的高程，既不是数字表面模型(DSM)，也不是数字地面模型(DTM)。要满足高质量的应用需求，开展林下地形(即DTM)测绘是很有必要的。考虑到单极化的TanDEM-X双站干涉数据主导了TanDEM-X DEM的全球测量，这为开展全球林下地形测绘提供了一个很好的契机。
当前，估计林下地形的方法可以分为四类。第一类是基于植被相关的产品校正InSAR DEM。最初的策略比较简单，即根据先验信息从DEM中减去植被高的平均值。显然，这种方式是粗糙的，因为它没有考虑植被密度和地形等因素。因此，一个改进方法被提出，主要思路是基于LiDAR、ICESAT和低分辨率植被产品(如冠层密度产品等)对InSAR DEM数据进行动态改正。然而，相对于大范围、高精度的林下地形需求，LiDAR和ICESAT数据的空间分布是相对稀疏的，且所采用的植被产品空间分辨率较低。而且，由于没有考虑SAR数据的未知穿透深度，这种方法是有偏的。第二类方法依赖于表征不同相位中心的干涉复相干系数。该复相干系数可以由全极化干涉(PolInSAR)数据提供，也可以由子孔径干涉数据提供。在复数单位圆中通过对离散的相位中心进行直线拟合探测地表相位。虽然直线拟合方法已经广泛被采用，但是它的性能与相位中心的离散度密切相关。在近似相等的相位中心或者严重的时间去相干条件下，直线拟合方法难以得到令人满意的结果。而且，该方法更多的用于机载全极化干涉数据，而在星载数据中应用很少，主要原因在于同时具备多极化功能和合适的InSAR观测几何的星载SAR系统很少。所以，大范围的林下地形反演受限。第三类方法是使用层析SAR的方法构建层析谱提取林下地形。为了达到此目的，植被的三维垂直结构需要先被重建，这依赖于多观测角度干涉数据。由于层析SAR技术对于观测数据要求比较严格，例如足够的高度向分辨率、时间去相干和大气影响较小等，因此该方法在星载情况下同样是严重受限的。而且，当地形存在明显起伏时，三维垂直结构的重建会受到影响，因此林下地形提取也同样受到影响。第四类方法被称为InSAR后向散射系数-高程映射法，这种方法类似于层析谱和LiDAR的全波形，通过捕捉在高度向的后向散射系数分布测量获取森林高度和林下地形。尽管该类方法不需要假设和散射模型，但是对使用数据同样要求比较严格，除了需要多基线数据，还需要满足低信噪比和低斑点噪声，足够的垂直向有效波数和可允许探测的坡度等。
综上所述，当前能够用于开展林下地形测绘的星载SAR系统是很有限的，要么是因为不合适的观测几何，要么是严重的时间去相干，要么就是不满足全极化获取。考虑到TanDEM-X双站数据的高质量干涉以及全球可用性，有必要设计一种专门适合于TanDEM-X双站干涉的大范围林下地形估计方法，从而为全球林下地形获取奠定基础。
发明内容
本发明所要解决的技术问题是突破InSAR技术难以进行林下地形测量的瓶颈，提供一种基于TanDEM-X双站单基线干涉数据的大范围林下地形估计方法、装置、设备及介质，可以实现大范围、高精度的林下地形估计，有助于实现全球林下地形的高精度重建。
为实现上述技术目的，本发明采用如下技术方案：
一种基于TanDEM-X双站单基线干涉数据的大范围林下地形估计方法，包括以下步骤：
步骤1，获取TanDEM-X双站SLC数据进行配准、重采样和干涉处理，生成干涉相位和相干幅度；
步骤2，使用外部DEM数据模拟平地相位和地形相位，并利用模拟平地相位和地形相位与步骤1得到的干涉相位之间的关系，计算地形残差Δh；然后使用得到的地形残差Δh对外部DEM数据进行补偿，得到TanDEM-X双站干涉条件下的高精度的InSAR DEM；
所述地形残差Δh是指外部DEM数据与InSAR真实测量的DEM之间的差异；
步骤3，首先，根据步骤1得到的相干幅度，计算X波段在森林研究区域全场景的穿透深度d；然后，从外部引入稀疏分布的先验点森林高度信息，并利用先验点的森林高度及穿透深度，拟合有关于穿透深度与森林高度之间的线性关系；再后，基于拟合的线性关系和研究区域全场景的穿透深度d，计算研究区域全场景的森林高度；最后，将研究区域全场景的森林高度与穿透深度做差计算散射相位中心高度；
步骤4，从步骤2得到的TanDEM-X双站干涉条件下的高精度的InSAR DEM中扣除步骤3得到的散射相位中心高度，即得到森林研究区域全场景的林下地形。
进一步的，步骤1得到的干涉相位与平地相位和地形相位之间的关系表达式为：
Figure GDA0003792870990000031
式中，
Figure GDA0003792870990000032
表示TanDEM-X双站干涉相位，
Figure GDA0003792870990000033
和
Figure GDA0003792870990000034
表示平地相位和地形相位，
Figure GDA0003792870990000035
表示轨道误差，
Figure GDA0003792870990000036
表示噪声相位；
步骤2中，计算地形残差Δh的具体过程为：
步骤2.1，使用外部DEM数据同时模拟平地相位和地形相位之和：
Figure GDA0003792870990000037
Figure GDA0003792870990000038
Figure GDA0003792870990000039
式中，
Figure GDA00037928709900000310
表示模拟的平地相位
Figure GDA00037928709900000311
和模拟的地形相位
Figure GDA00037928709900000312
之和，λ表示波长；R1和R20分别表示主、从影像的SAR传感器天线中心到地面目标的斜距；B代表基线长；α代表基线倾角；θ代表考虑地形的视角，θ0代表参考视角，Δθ代表地形相关的视角与参考视角之差；h代表地形高程，由外部DEM数据提供；
步骤2.2，从干涉相位中去除模拟相位
Figure GDA00037928709900000313
得到差分干涉相位，并进行相位解缠，得到解缠差分干涉相位
Figure GDA00037928709900000314
即
Figure GDA00037928709900000315
式中，
Figure GDA00037928709900000316
是地形残差相位，
Figure GDA00037928709900000317
由外部DEM与InSAR真实测量的DEM之间的地形残差Δh引起，可以表达为：
Figure GDA00037928709900000318
式中，B⊥代表与地形相关的垂直基线；KZ代表垂直向有效波数；
步骤2.3，根据InSAR几何关系在每一距离行构建轨道误差与基线参数(B，α)误差之间的函数关系如公式(7)，并逐行开展基线参数误差估计，进而得到轨道误差相位
Figure GDA00037928709900000319
式中，ΔB代表基线长误差；Δα代表基线倾角误差；ε代表模型误差，与相位解缠引起的相位偏差有关；
步骤2.4，从解缠差分相位
Figure GDA0003792870990000041
中改正轨道误差相位
Figure GDA0003792870990000042
得到包括地形残差Δh的相位，然后通过滤波处理去除噪声相位
Figure GDA0003792870990000043
最终通过相位-高度转换得到Δh，如公式(8)所示：
Figure GDA0003792870990000044
式中，filter(·)代表滤波处理。
进一步的，步骤3中，X波段在森林研究区域全场景的穿透深度d的计算方法为：
TanDEM-X双站干涉条件下，X波段与森林冠层的相互作用，其复相干系数γ符合以下干涉散射模型：
Figure GDA0003792870990000045
式中，
Figure GDA0003792870990000046
代表地表相位，γvol代表纯体相干性，表达为
Figure GDA0003792870990000047
式中，d代表穿透深度，hv代表森林高度，KZ表示垂直向有效波数；
基于X波段在森林中的穿透性，即X波段在森林的散射相位中心大概率集中在森林冠层内，而且相位中心位置几乎在植被的一半高度以上，可得森林高度hv与穿透深度d的比值RFP≥2，故可假设干涉散射模型中的森林高度为无限大；通过对公式(9)两侧取模可得相干幅度与穿透深度之间的一一对应的函数关系，即为穿透深度模型，
Figure GDA0003792870990000048
最终，使用垂直向有效波数KZ和步骤1得到的相干幅度|γ|，根据公式(11)所示的穿透深度模型计算得到穿透深度d。
进一步的，根据外部DEM数据对垂直向有效波数Kz进行校正后，再用于穿透深度模型中计算穿透深度。
进一步的，所述从外部引入稀疏分布的先验森林高度信息，由ICESAT-2或者激光雷达测量数据提供。
进一步的，所述外部DEM数据为全球发布的DEM数据，可以是ASTER DEM、SRTM DTM、TanDEM-X DEM中的任意一种。
本发明还提供一种基于TanDEM-X双站单基线干涉数据的大范围林下地形估计装置，包括：
干涉相位和相干幅度获取模块，用于：获取TanDEM-X双站SLC数据进行配准、重采样和干涉处理，生成干涉相位和相干幅度；
高精度InSAR DEM数据计算模块，用于：使用外部DEM数据模拟平地相位和地形相位，并利用模拟平地相位和地形相位与干涉相位之间的关系，计算地形残差Δh；然后使用得到的地形残差Δh对外部DEM数据进行补偿，得到TanDEM-X双站干涉条件下的高精度的InSAR DEM；
其中，所述地形残差Δh是指外部DEM数据与InSAR真实测量的DEM之间的差异；
散射相位中心高度计算模块，用于：首先，提取X波段对应的相干幅度，计算X波段在森林研究区域全场景的穿透深度d；然后，从外部引入稀疏分布的先验点森林高度信息，并利用先验点的森林高度及穿透深度，拟合有关于穿透深度与森林高度之间的线性关系；再后，基于拟合的线性关系和研究区域全场景的穿透深度d，计算研究区域全场景的森林高度；最后，将研究区域全场景的森林高度与穿透深度做差计算散射相位中心高度；
林下地形重建模块，用于：将TanDEM-X双站干涉条件下的高精度的InSAR DEM数据扣除散射相位中心高度，即得到森林研究区域全场景的林下地形。
本发明还提供一种设备，包括处理器和存储器；其中：所述存储器用于存储计算机指令；所述处理器用于执行所述存储器存储的计算机指令，具体执行上述任一技术方案所述的方法。
本发明还提供一种计算机存储介质，用于存储程序，所述程序被执行时，用于实现上述任一技术方案所述的方法。
有益效果
本发明的有益效果包括：
(1)通过先验点的森林高度与穿透深度在坐标平面内的线性函数关系，揭示了二者之间存在显著的线性函数关系，从而实现根据穿透深度计算研究区域全场景的森林高度；
(2)通过融合双站单基线InSAR干涉对和稀疏分布的先验点森林高度信息，实现对大范围的林下地形进行高精度估计；
(3)仅通过一个双站干涉对，然后融合稀疏分布的ICESAT-2(或LiDAR)先验数据，可以同时获取InSAR DEM产品、森林高度产品、DSM产品和DTM产品4个产品，这些都是当前国土、林业、测绘、矿业以及交通等多个部门比较迫切需求的。
附图说明
图1为本发明实施例所述方法的流程图；
图2为本发明实施例中ICESAT-2点的空间分布情况；
图3中，(a)为本发明实施例中步骤3中得到的穿透深度，(b)为先验点的森林高度与穿透深度的分布以及拟合得到的线性关系；
图4中，(a)为本发明实施案例中步骤2重建得到的高精度InSAR DEM图，(b)为本发明实施案例中步骤4得到的森林研究区域全场景的林下地形图；
图5中，(a)为本发明实施案例中步骤2重建得到的高精度InSAR DEM与LiDAR林下地形验证数据的交叉对比，(b)本发明实施案例中步骤4重建得到的高精度林下地形与LiDAR林下地形验证数据的交叉对比；
图6为从本发明实施例得到的林下地形结果图中选取的剖面展示图。
具体实施方式
为了使本发明的目的、技术方案及其优点更加清楚明白，采用德国宇航局通过TerrSAR-X/TanDEM-X系统获取的双站干涉数据，具体参数见附表1，对本发明进行深入详细地说明。TanDEM-X双站数据研究区(见附图2)的地形起伏为海拔100-400m，且该地区为典型的北方针叶森林覆盖区，森林平均树高为18m，最高的树高为30m。附图2中展示了ICESAT-2点的空间分布情况，同时为了验证林下地形重建结果，也展示了LiDAR验证数据的覆盖区域。此处所描述的具体实施案例仅用以解释本发明，并不用于限定本发明。
基于上述所采用的TanDEM-X双站单基线干涉数据，本发明提供的基于TanDEM-X双站单基线干涉数据的大范围林下地形估计方法实施，参考图1所示，其具体实施步骤如下：
步骤1：先对获取的TanDEM-X双站SLC数据进行配准、多视处理(方位向x距离向为6x6多视)使其空间分辨率约为12米x12米，重采样和干涉处理，生成干涉相位和相干幅度。
其中，得到的干涉相位与平地相位和地形相位之间的关系表达式为：
Figure GDA0003792870990000061
式中，
Figure GDA0003792870990000062
表示TanDEM-X双站干涉相位，
Figure GDA0003792870990000063
和
Figure GDA0003792870990000064
表示平地相位和地形相位，
Figure GDA0003792870990000065
表示轨道误差，
Figure GDA0003792870990000066
表示噪声相位。
步骤2：TanDEM-X双站干涉条件下的高精度InSAR DEM的生成。
在该步骤2中，使用外部DEM数据模拟平地相位和地形相位，并利用模拟平地相位和地形相位与步骤1得到的干涉相位之间的关系，计算地形残差Δh；然后使用得到的地形残差Δh对外部DEM数据进行补偿，得到TanDEM-X双站干涉条件下的高精度的InSAR DEM；具体地，可通过改正TanDEM-X双站干涉图中的轨道误差，从而获取准确的DEM残差(即下文所述的地形残差Δh)，进而生成高精度的InSAR DEM。
计算地形残差Δh的具体过程为：
步骤2.1，使用外部DEM数据同时模拟平地相位和地形相位之和：
Figure GDA0003792870990000071
Figure GDA0003792870990000072
Figure GDA0003792870990000073
式中，
Figure GDA0003792870990000074
表示模拟的平地相位
Figure GDA0003792870990000075
和模拟的地形相位
Figure GDA0003792870990000076
之和，λ表示波长；R1和R20分别表示主、从影像的SAR传感器天线中心到地面目标的斜距；B代表基线长；α代表基线倾角；θ代表考虑地形的视角，θ0代表参考视角，Δθ代表地形相关的视角与参考视角之差；h代表地形高程，由外部DEM数据提供；
其中的外部DEM数据为全球发布的DEM数据，可以是ASTER DEM、SRTM DTM、TanDEM-X DEM等任意一种。
步骤2.2，从干涉相位中去除模拟相位
Figure GDA0003792870990000077
得到差分干涉相位，并进行相位解缠，得到解缠差分干涉相位
Figure GDA0003792870990000078
即
Figure GDA0003792870990000079
式中，
Figure GDA00037928709900000710
是地形残差相位，
Figure GDA00037928709900000711
由外部DEM与InSAR真实测量的DEM之间的地形残差Δh引起，可以表达为：
Figure GDA00037928709900000712
式中，B⊥代表与地形相关的垂直基线；KZ代表垂直向有效波数；
步骤2.3，根据InSAR几何关系在每一距离行构建轨道误差与基线参数(B，α)误差之间的函数关系如公式(7)，并逐行开展基线参数误差估计，进而得到轨道误差相位：
Figure GDA00037928709900000713
式中，ΔB代表基线长误差；Δα代表基线倾角误差；ε代表模型误差，与相位解缠引起的相位偏差有关；
基线参数(B，α)误差估计属于现有技术，比如论文《Modeling and RobustEstimation for the Residual Motion Error in Airborne SAR Interferometry》中记载的方法，本发明不做具体阐述。
步骤2.4，从解缠差分相位
Figure GDA00037928709900000714
中改正轨道误差相位
Figure GDA00037928709900000715
得到包括地形残差Δh的相位，然后通过滤波处理去除噪声相位
Figure GDA00037928709900000716
最终通过相位-高度转换得到Δh，如公式(8)所示：
Figure GDA00037928709900000717
式中，filter(·)代表滤波处理。
最终将得到地形残差Δh加到外部DEM上完成InSAR DEM重建，即得到TanDEM-X双站干涉条件下的高精度InSAR DEM，见附图4(a)。
步骤3：穿透深度、森林高度估计以及散射相位中心高度的计算。具体如下，
首先，穿透深度的计算：
TanDEM-X双站干涉条件下，X波段与森林冠层的相互作用，且几乎很少穿过冠层到达地表，因此复相干系数γ符合干涉散射模型：
Figure GDA0003792870990000081
式中，
Figure GDA0003792870990000082
代表地表相位，γvol代表纯体相干性，表达为
Figure GDA0003792870990000083
式中，d代表穿透深度，hv代表森林高度，KZ表示垂直向有效波数；
基于X波段在森林中的穿透性，即X波段在森林的散射相位中心大概率集中在森林冠层内，而且相位中心位置几乎在植被的一半高度以上，可得森林高度hv与穿透深度d的比值RFP≥2，故可假设干涉散射模型中的森林高度为无限大；通过对公式(9)两侧取模可得相干幅度与穿透深度之间的一一对应的函数关系，即为穿透深度模型，
Figure GDA0003792870990000084
因此，使用垂直向有效波数KZ和步骤1得到的相干幅度|γ|，根据公式(11)所示的穿透深度模型计算得到穿透深度d，见附图3(a)。
在本实施例中，为了削弱受地形坡度影响的相干性对穿透深度估计产生的偏差，垂直向有效波数Kz可以根据所用的TanDEM-XDEM数据进行校正，校正方法可采用现有技术，比如论文《Canopypenetration depth estimation with TanDEM-X and itscompensation in temperateforests》中记载的方法。
然后，线性关系的拟合：
考虑到穿透深度相对于森林高度呈现一定的比例，本实施例引入稀疏分布的ICESAT-2(或者激光雷达测量数据，如LiDAR)提供部分先验森林高度信息，通过穿透深度与森林高度坐标平面内的先验分布拟合线性关系，见附图3(b)。
再后，森林高度的估计：
基于前述拟合的有关于穿透深度与森林高度坐标平面内的先验分布拟合线性关系，根据研究区域全场景的穿透深度d，完成对研究区域全场景的森林高度估计。
最后，将研究区域全场景的森林高度与穿透深度做差计算散射相位中心高度。
步骤4：林下地形的重建；
TanDEM-X双站干涉条件下的高精度的InSARDEM中扣除步骤3得到的散射相位中心高度，即得到森林研究区域全场景的林下地形，见附图4(b)。
相应于上述所述的大范围林下地形估计方法实施例，本发明还提供一种基于TanDEM-X双站单基线干涉数据的大范围林下地形估计装置，包括：
干涉相位和相干幅度获取模块，用于：获取TanDEM-X双站SLC数据进行配准、重采样和干涉处理，生成干涉相位和相干幅度；
高精度InSAR DEM数据计算模块，用于：使用外部DEM数据模拟平地相位和地形相位，并利用模拟平地相位和地形相位与干涉相位之间的关系，计算地形残差Δh；然后使用得到的地形残差Δh对外部DEM数据进行补偿，得到TanDEM-X双站干涉条件下的高精度的InSAR DEM；
其中，所述地形残差Δh是指外部DEM数据与InSAR真实测量的DEM之间的差异；
散射相位中心高度计算模块，用于：首先，提取X波段对应的相干幅度，计算X波段在森林研究区域全场景的穿透深度d；然后，从外部引入稀疏分布的先验点森林高度信息，并利用先验点的森林高度及穿透深度，拟合有关于穿透深度与森林高度之间的线性关系；再后，基于拟合的线性关系和研究区域全场景的穿透深度d，计算研究区域全场景的森林高度；最后，将研究区域全场景的森林高度与穿透深度做差计算散射相位中心高度；
林下地形重建模块，用于：从TanDEM-X双站干涉条件下得到的高精度的InSAR DEM中扣除散射相位中心高度，即得到森林研究区域全场景的林下地形。
本发明还提供一种设备实施例，包括处理器和存储器；其中：所述存储器用于存储计算机指令；所述处理器用于执行所述存储器存储的计算机指令，具体执行上述的大范围林下地形估计方法。
本发明还提供一种计算机存储介质，用于存储程序，所述程序被执行时，用于实现上述的大范围林下地形估计方法。
为了验证重建的林下地形的可靠性，本发明引入LiDAR林下地形作为验证数据，交叉验证结果见附图5。结果显示，在没有进行散射相位中心高度改正时，即TanDEM-X InSARDEM相对于LiDAR林下地形的验证精度为7.76m，如图5(a)所示；经过散射相位中心改正后，重建的林下地形精度为3.31m，如图5(b)所示，总体精度提高了57.3％。为了进一步验证林下地形结果的可靠性，在验证区内选取了一个剖面，结果展示见附图6。结果显示估计的林下地形结果与LiDAR林下地形最为接近，证实了估计的林下地形的可靠性。
表1 TanDEM-X双站单基线干涉数据的干涉几何参数
Figure GDA0003792870990000101
以上实施例为本申请的优选实施例，本领域的普通技术人员还可以在此基础上进行各种变换或改进，在不脱离本申请总的构思的前提下，这些变换或改进都应当属于本申请要求保护的范围之内。

Claims (9)

1.一种基于TanDEM-X双站单基线干涉数据的大范围林下地形估计方法，其特征在于，包括以下步骤：
步骤1，获取TanDEM-X双站SLC数据进行配准、重采样和干涉处理，生成干涉相位和相干幅度；
步骤2，使用外部DEM数据模拟平地相位和地形相位，并利用模拟平地相位和地形相位与步骤1得到的干涉相位之间的关系，计算地形残差Δh；然后使用得到的地形残差Δh对外部DEM数据进行补偿，得到TanDEM-X双站干涉条件下的高精度的InSAR DEM；
所述地形残差Δh是指外部DEM数据与InSAR真实测量的DEM之间的差异；
步骤3，首先，根据步骤1得到的相干幅度，计算X波段在森林研究区域全场景的穿透深度d；然后，从外部引入稀疏分布的先验点森林高度信息，并利用先验点的森林高度及穿透深度，拟合有关于穿透深度与森林高度之间的线性关系；再后，基于拟合的线性关系和研究区域全场景的穿透深度d，计算研究区域全场景的森林高度；最后，将研究区域全场景的森林高度与穿透深度做差计算散射相位中心高度；
步骤4，从步骤2得到的TanDEM-X双站干涉条件下的高精度的InSAR DEM中扣除步骤3得到的散射相位中心高度，即得到森林研究区域全场景的林下地形。
2.根据权利要求1所述的方法，其特征在于，步骤1得到的干涉相位与平地相位和地形相位之间的关系表达式为：
Figure FDA0003792870980000011
式中，
Figure FDA0003792870980000012
表示TanDEM-X双站干涉相位，
Figure FDA0003792870980000013
和
Figure FDA0003792870980000014
表示平地相位和地形相位，
Figure FDA0003792870980000015
表示轨道误差，
Figure FDA0003792870980000016
表示噪声相位；
步骤2中，计算地形残差Δh的具体过程为：
步骤2.1，使用外部DEM数据同时模拟平地相位和地形相位之和：
Figure FDA0003792870980000017
Figure FDA0003792870980000018
Figure FDA0003792870980000019
式中，
Figure FDA00037928709800000110
表示模拟的平地相位
Figure FDA00037928709800000111
和模拟的地形相位
Figure FDA00037928709800000112
之和，λ表示波长；R1和R20分别表示主、从影像的SAR传感器天线中心到地面目标的斜距；B代表基线长；α代表基线倾角；θ代表考虑地形的视角，θ0代表参考视角，Δθ代表地形相关的视角与参考视角之差；h代表地形高程，由外部DEM数据提供；
步骤2.2，从干涉相位中去除模拟相位
Figure FDA0003792870980000021
得到差分干涉相位，并进行相位解缠，得到解缠差分干涉相位
Figure FDA0003792870980000022
即
Figure FDA0003792870980000023
式中，
Figure FDA0003792870980000024
是地形残差相位，
Figure FDA0003792870980000025
由外部DEM与InSAR真实测量的DEM之间的地形残差Δh引起，表达为：
Figure FDA0003792870980000026
式中，B⊥代表与地形相关的垂直基线；Kz代表垂直向有效波数；
步骤2.3，根据InSAR几何关系在每一距离行构建轨道误差与基线参数(B，α)误差之间的函数关系如公式(7)，并逐行开展基线参数误差估计，进而得到轨道误差相位
Figure FDA0003792870980000027
式中，ΔB代表基线长误差；Δα代表基线倾角误差；ε代表模型误差，与相位解缠引起的相位偏差有关；
步骤2.4，从解缠差分相位
Figure FDA0003792870980000028
中改正轨道误差相位
Figure FDA0003792870980000029
得到包括地形残差Δh的相位，然后通过滤波处理去除噪声相位
Figure FDA00037928709800000210
最终通过相位-高度转换得到Δh，如公式(8)所示：
Figure FDA00037928709800000211
式中，filter(·)代表滤波处理。
3.根据权利要求1所述的方法，其特征在于，步骤3中，X波段在森林研究区域全场景的穿透深度d的计算方法为：
TanDEM-X双站干涉条件下，X波段与森林冠层的相互作用，其复相干系数γ符合以下干涉散射模型：
Figure FDA00037928709800000212
式中，
Figure FDA00037928709800000213
代表地表相位，γvol代表纯体相干性，表达为
Figure FDA00037928709800000214
式中，d代表穿透深度，hv代表森林高度，Kz表示垂直向有效波数；
基于X波段在森林中的穿透性，即X波段在森林的散射相位中心大概率集中在森林冠层内，而且相位中心位置在植被的一半高度位置的上方，可得森林高度hv与穿透深度d的比值RFP≥2，故可假设干涉散射模型中的森林高度为无限大；通过对公式(9)两侧取模可得相干幅度与穿透深度之间的一一对应的函数关系，即为穿透深度模型，
Figure FDA0003792870980000031
最终，使用垂直向有效波数Kz和步骤1得到的相干幅度|γ|，根据公式(11)所示的穿透深度模型计算得到穿透深度d。
4.根据权利要求3所述的方法，其特征在于，根据外部DEM数据对垂直向有效波数Kz进行校正后，再用于穿透深度模型中计算穿透深度。
5.根据权利要求1所述的方法，其特征在于，所述从外部引入稀疏分布的先验森林高度信息，由ICESAT-2或者激光雷达测量数据提供。
6.根据权利要求1所述的方法，其特征在于，所述外部DEM数据为全球发布的DEM数据，可以是ASTER DEM、SRTM DTM、TanDEM-X DEM中的任意一种。
7.一种基于TanDEM-X双站单基线干涉数据的大范围林下地形估计装置，其特征在于，包括：
干涉相位和相干幅度获取模块，用于：获取TanDEM-X双站SLC数据进行配准、重采样和干涉处理，生成干涉相位和相干幅度；
高精度InSAR DEM数据计算模块，用于：使用外部DEM数据模拟平地相位和地形相位，并利用模拟平地相位和地形相位与干涉相位之间的关系，计算地形残差Δh；然后使用得到的地形残差Δh对外部DEM数据进行补偿，得到TanDEM-X双站干涉条件下的高精度的InSARDEM；
其中，所述地形残差Δh是指外部DEM数据与InSAR真实测量的DEM之间的差异；
散射相位中心高度计算模块，用于：首先，提取X波段对应的相干幅度，计算X波段在森林研究区域全场景的穿透深度d；然后，从外部引入稀疏分布的先验点森林高度信息，并利用先验点的森林高度及穿透深度，拟合有关于穿透深度与森林高度之间的线性关系；再后，基于拟合的线性关系和研究区域全场景的穿透深度d，计算研究区域全场景的森林高度；最后，将研究区域全场景的森林高度与穿透深度做差计算散射相位中心高度；
林下地形重建模块，用于：将TanDEM-X双站干涉条件下的高精度的InSAR DEM数据扣除散射相位中心高度，即得到森林研究区域全场景的林下地形。
8.一种设备，其特征在于，包括处理器和存储器；其中：所述存储器用于存储计算机指令；所述处理器用于执行所述存储器存储的计算机指令，具体执行如权利要求1-6任一所述的方法。
9.一种计算机存储介质，其特征在于，用于存储程序，所述程序被执行时，用于实现如权利要求1-6任一所述的方法。
CN202010771942.5A 2020-08-04 2020-08-04 一种大范围林下地形估计方法、装置、设备及介质 Active CN113341410B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010771942.5A CN113341410B (zh) 2020-08-04 2020-08-04 一种大范围林下地形估计方法、装置、设备及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010771942.5A CN113341410B (zh) 2020-08-04 2020-08-04 一种大范围林下地形估计方法、装置、设备及介质

Publications (2)

Publication Number Publication Date
CN113341410A CN113341410A (zh) 2021-09-03
CN113341410B true CN113341410B (zh) 2022-11-04

Family

ID=77467475

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010771942.5A Active CN113341410B (zh) 2020-08-04 2020-08-04 一种大范围林下地形估计方法、装置、设备及介质

Country Status (1)

Country Link
CN (1) CN113341410B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117452432B (zh) * 2023-12-21 2024-03-15 西南林业大学 一种基于森林穿透补偿的森林冠层高度估测方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5552787A (en) * 1995-10-10 1996-09-03 The United States Of America As Represented By The Secretary Of The Navy Measurement of topography using polarimetric synthetic aperture radar (SAR)
CN101078769B (zh) * 2006-05-25 2010-06-16 中国科学院中国遥感卫星地面站 单次全极化合成孔径雷达图像反演数字高程模型的方法
JP4982769B2 (ja) * 2009-03-25 2012-07-25 国立大学法人長岡技術科学大学 作物の生育診断方法及び生育診断システム
ES2384922B1 (es) * 2010-06-07 2013-06-11 Universitat Politècnica De Catalunya Procedimiento para la estimación de la topografía de la superficie de la tierra en áreas con cobertura vegetal.
CN105005047B (zh) * 2015-07-17 2017-06-09 武汉大学 后向散射优化的森林复杂地形校正及树高反演方法、系统
CN105204079B (zh) * 2015-08-31 2018-05-22 中国科学院测量与地球物理研究所 一种利用TanDEM-X双站InSAR提取地震滑坡体积的方法
CN111273293B (zh) * 2020-03-03 2021-11-23 中南大学 一种顾及地形起伏的InSAR残余运动误差估计方法及装置
CN111398957B (zh) * 2020-04-01 2022-08-02 中国林业科学研究院资源信息研究所 改进相干性计算方法的短波长双天线InSAR森林高度反演方法

Also Published As

Publication number Publication date
CN113341410A (zh) 2021-09-03

Similar Documents

Publication Publication Date Title
CN106772342B (zh) 一种适用于大梯度地表沉降监测的时序差分雷达干涉方法
Liu et al. Estimating Spatiotemporal Ground Deformation With Improved Persistent-Scatterer Radar Interferometry $^\ast$
Samsonov et al. A simultaneous inversion for deformation rates and topographic errors of DInSAR data utilizing linear least square inversion technique
CN111273293B (zh) 一种顾及地形起伏的InSAR残余运动误差估计方法及装置
CN111239736B (zh) 基于单基线的地表高程校正方法、装置、设备及存储介质
CN108663678B (zh) 基于混合整数优化模型的多基线InSAR相位解缠算法
Tang et al. Atmospheric correction in time-series SAR interferometry for land surface deformation mapping–A case study of Taiyuan, China
Wang et al. Demonstration of time-series InSAR processing in Beijing using a small stack of Gaofen-3 differential interferograms
CN112711021A (zh) 一种多分辨率InSAR交互干涉时序分析方法
CN113341410B (zh) 一种大范围林下地形估计方法、装置、设备及介质
Mao et al. Estimation and compensation of ionospheric phase delay for multi-aperture InSAR: An azimuth split-spectrum interferometry approach
Yu et al. Digital Elevation Model generation using ascending and Descending multi-baseline ALOS/PALSAR radar images
CN112363165B (zh) 一种林下地形反演方法、装置、设备及介质
Haque et al. Time series analysis of subsidence in Dhaka City, Bangladesh using Insar
Kim et al. Ground deformation tracking over Mt. Baekdu: A pre-evaluation of possible magma recharge by D-InSAR analysis
Yamashita et al. Mitigation of ionospheric noise in azimuth offset based on the split-spectrum method
Thapa et al. Estimation of atmospheric effects of RADARSAT-2 D-InSAR product using groundbasedand spaceborne meterological data
CN111580101A (zh) 一种基于外部DEM的InSAR基线误差无控改正方法及装置
Wieczorek EVALUATION OF DEFORMATIONS IN THE URBAN AREA OF OLSZTYN USING SENTINEL-1 SAR INTERFEROMETRY.
Lombardi et al. Accuracy of high resolution CSK interferometric Digital Elevation Models
Biswas et al. Investigation of ground deformation with PSInSAR approach in an unstable urban area Naples, Italy using X-band SAR images
Wu et al. Regression analysis of errors of sar-based dems and controlling factors
Slacikova et al. Evaluation of interpolation methods in InSAR DEM derivation from ERS tandem data
Lanari et al. A differential SAR interferometry approach for monitoring urban deformation phenomena
Nitti et al. DEM generation by using COSMO/SkyMed tandem pairs and weather models

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant