基于TanDEM-X相干系数的森林高度估测方法*

2020-07-28范亚雄陈尔学李增元张王菲金玉栋蔡丽杰

林业科学 2020年6期

范亚雄陈尔学李增元赵磊张王菲金玉栋蔡丽杰

(1. 中国林业科学研究院资源信息研究所北京 100091； 2. 西南林业大学林学院昆明 650224； 3. 喀喇沁旗旺业甸实验林场赤峰 024423)

作为重要的森林垂直结构参数，林分优势木平均高(以下简称森林高度)是估测森林材积生产潜力的主要指标，且与森林地上生物量(above ground biomass，AGB)、碳储量等紧密联系，准确、大面积地获取森林高度信息对于森林的精细化经营管理、碳循环和气候变化科学研究等均具有重要意义(陈尔学等， 2007)。传统基于抽样标准样本调查森林高度的方法费时费力，也难以得到空间连续的森林高度测量结果，已无法满足现代森林资源经营管理和生态环境科学研究的需求。目前，具有区域森林高度制图潜力的遥感手段主要有激光雷达(light detection and ranging，LiDAR)、立体摄影测量和干涉合成孔径雷达(interferometric synthetic aperture radar，InSAR)3种，其中，LiDAR和立体摄影测量由于易受云雾等不利天气状况的影响，限制了其大范围、连续性森林高度制图能力，而InSAR具有SAR全天候、全天时成像的优势，同时对植被垂直结构敏感，是大面积估测森林高度的有效手段之一。

当前在轨星载SAR系统的InSAR数据通常需要重轨获取，植被区域受时间失相干影响严重，不利于森林高度估测。而TanDEM-X(TerraSAR-X add-on for digital elevation measurement)系统由2颗相似的卫星在轨道中以几百米间距并行，采用一发双收模式对地物进行重复观测，时间基线可以忽略，失相干主要由植被体散射效应引起(Kriegeretal.， 2007)，因此在森林垂直结构估测中具有很大潜力。国内外基于TanDEM-X InSAR数据进行森林高度估测的研究主要有2种思路： 1) 基于相位信息的DSM-DEM差分法，将InSAR获得的DSM与已知高精度的DEM相减得到有效散射中心高度，该高度与真实森林高度具有很好的相关性，但有效散射中心的具体位置与森林结构和微波频率有关，往往需要实测数据对其进行标定，才能得到森林高度估测结果(Sojaetal.， 2013； Sadeghietal.， 2014)； 2) 基于相干散射模型，通过物理模型将相干性与森林垂直结构相联系，由相干性求解森林高度，其中最常用的是随机体散射-地表散射模型(random volume over ground, RVoG)，其假设地表散射贡献可以忽略，首先利用已知DEM估算地表相位，然后构建森林高度-消光系数的二维查找表反演森林高度和消光系数(Praksetal.， 2012； Kugkeretal.， 2014)。除了RVoG模型外，Soja等(2014； 2015)还基于干涉水云模型(interferometric water cloud model，IWCM)建立了一种2层散射模型，在已知林下地表相位的情况下，该模型可同时反演森林高度和冠层密度。但需要说明的是，由于X-波段的微波对森林穿透能力较差，难以从干涉数据本身提取准确的林下地形，因此在估测森林高度时通常需要其他来源的高精度DEM，在应用中受到很大限制。

基于X-波段InSAR估测森林高度时，往往假设地表散射贡献为零，此时若进一步假设森林内部的平均消光系数为零，则可直接由相干系数得到森林高度(Cloude， 2010)，该方法称为SINC模型。SINC模型本质上是对RVoG模型的简化，但该模型不依赖高精度DEM，与以往研究方法相比具有较大优势。Cloude等(2014)基于TanDEM-X双极化数据信噪比优化后得到的相干系数，分别在北方针叶林和温带雨林验证SINC模型的可行性，结果发现SINC模型在北方针叶林可取得良好估测结果，但在温带雨林由于较高消光系数的影响，森林高度估测精度较低。冯琦等(2016)利用机载X-波段双天线InSAR数据估测森林高度，结果表明SINC模型可取得较好的估测精度，且相比差分法具有不依赖高精度DEM的优势。Chen等(2016)研究指出，TanDEM-X单极化数据相比双极化数据具有更高的信噪比，只利用单极化数据的相干系数就可取得较好的森林高度估测结果。Olesk等(2015)重点分析季节变化对估测结果的影响，结果表明，对于落叶树种，冬季数据的估测精度要优于夏季数据，而常绿树种对季节变化不敏感。在Olesk等(2016)的研究中，通过将SINC模型参数化，使其适应不同的季节条件和林分条件，在估测森林高度时更加稳健。

以往研究论述了SINC模型在森林高度估测中的可行性，但在计算局部垂直波数(kz)对地形进行补偿和地理编码时往往采用LiDAR提取的高分辨率DEM作为参考。目前，国内高分辨率DEM通常由机载LiDAR获取，覆盖范围有限且成本高昂，而中等分辨率SRTM DEM则覆盖全球大部分陆地，因此分析高分辨率与中等分辨率的DEM对SINC模型森林高度估测精度的影响很有必要。鉴于此，本研究采用TanDEM-X单极化InSAR数据，以内蒙古根河市北部天然林为研究区，研究基于相干系数的SINC模型森林高度估测方法，并分析高分辨率(5 m)的LiDAR DEM和中等分辨率(30 m)的SRTM DEM对模型估测精度和稳定性的影响。

1 研究区概况与数据

1.1 研究区概况

研究区位于内蒙古根河市北部，121.42°—121.57°E，50.91°—50.98°N，属高纬度、高寒地区。呈大兴安岭山地地貌特征，海拔800～1 100 m，地势起伏相对平缓，大部分坡度在10°以下。该区属大兴安岭北部针叶林生态区，森林覆盖率75%以上，主要树种有兴安落叶松(Larixgmelinii)、樟子松(Pinussylvestrisvar.mongolica)、白桦(Betulaplatyphylla)和山杨(Populusdavidiana)等。在图1所示的Google Earth多光谱遥感影像上，叠加着InSAR数据的强度影像(主图像)，其中黄色矩形框代表机载LiDAR数据范围，红色矩形框代表研究区范围。

图1 TanDEM-X InSAR数据和机载LiDAR数据的覆盖范围Fig.1 The coverage of TanDEM-X InSAR data and airborne LiDAR data

1.2 InSAR数据及预处理

TanDEM-X是德国宇航局发射的高精度、高分辨率X-波段雷达干涉测量系统。本研究获取TanDEM-X条带模式下的升轨InSAR数据，成像时间为2012年8月14日，空间有效基线202.28 m，中心入射角41.40°，距离向和方位向的采样间隔分别为1.4和2.0 m。数据覆盖范围见图1中的灰色图，东西向约33 km，南北向约56 km。

获取的TanDEM-X单视复数数据(single look complex, SLC)已经过亚像元级精确配准，还需要进行干涉图生成、基线估计、相干性估计、相位解缠和地理编码等处理，详细流程如图2所示。

图2 Tandem-X InSAR数据处理流程Fig.2 Processing flowchart of TanDEM-X InSAR data

1.3 机载LiDAR数据及处理

机载LiDAR数据获取时间为2012年8月，以“运5”飞机为平台，采用Leica LiDAR系统开展飞行试验。平均飞行高度2 700 m，航带扫描宽度约1 000 m，波长1 550 nm，获取的激光点云密度平均每平方米5.6个。基于TerraSolid软件对点云数据进行噪声点去除及地面点和非地面点分类，分类后的地面点云数据以不规则三角网邻近像元内插算法生成DEM，非地面点云数据以最大高程内插算法生成DSM，将DSM与DEM相减，得到研究区的CHM，三者空间分辨率均为5 m。CHM代表冠层顶部高度，根据当地树高分布特点，将CHM中小于2 m和大于40 m的异常点剔除。穆喜云(2015)利用该数据提取CHM的研究表明，样地的胸高断面积加权高(Lorey’s高)与基于LiDAR CHM数据的算术平均高具有很高相关性(R2=0.834)，因此，本研究以一定大小的采样窗口对CHM中的有效值进行算术平均作为森林高度，验证估测结果。

1.4 SRTM DEM数据

SRTM是NASA在2000年运行航天飞机雷达地形测绘任务获得的北纬60°到南纬54°之间大部分陆地地区的DEM数据。本研究获取覆盖研究区30 m空间分辨率的SRTM DEM，采用双线性内插法将其过采样到5 m像元大小，以与LiDAR数据的分辨率保持一致，并比较SRTM DEM、LiDAR DEM和LiDAR DSM之间的高程差异。以图3a中的黑色实线为剖面，绘制三者的剖面图如图4所示。可以看到，LiDAR DSM整体要比LiDAR DEM和SRTM DEM高，且高程波动较大，体现了森林高度信息。LiDAR DEM和SRTM DEM的高程变化趋势则更相符，但由于SRTM DEM受植被偏差影响，所以在大部分区域要高于LiDAR DEM。

图4 不同DEM的剖面线对比Fig.4 Comparison of selected profiles from SRTM DEM, LiDAR derived DEM and DSM

图3 SRTM DEM与LiDAR获取的DEM、DSMFig.3 SRTM DEM and LiDAR derived DEM, DSM

SINC模型估测森林高度中DEM的作用主要体现在2方面：首先，局部垂直波数(kz)的估算需要坡度信息来补偿地形的影响；其次，DEM用于相干性、kz和森林高度等估测结果的地理编码，其质量直接影响地理编码的精度。本研究将分析LiDAR DEM和SRTM DEM的差异对SINC模型森林高度估测结果的影响。

2 研究方法

2.1 非体散射失相干校正

在干涉测量中，相干性γObs是对2次回波信号相似程度的度量，定义为复信号s1、s2的归一化互相关，如式(1)所示：

(1)

式中： *表示复共轭；〈…〉表示空间平均；s1、s2分别为在空间基线两端接收到的复信号。

通常认为，观测到的γObs由式(2)中几项失相干组成：

γObs=γSNRγBaseγProcγTmpγVol。

(2)

式中：γSNR表示信噪比失相干；γBase是由2次观测几何差异引起的基线失相干；γProc是由数据处理造成的失相干，如主辅影像配准误差；γTmp是由2次数据采集期间散射体(植被-地表)发生变化引起的时间失相干；γVol是由植被-地表引起的失相干。

由于TanDEM-X系统无时间基线，且获得的SLC数据已经过亚像元级的配准和距离向、方位向的频谱滤波处理，因此可将γBase、γProc和γTmp忽略(假设均等于1.0)，数据中的非体散射失相干只剩下γSNR，经式(3)、(4)校正后可得到γVol：

(3)

(4)

式中： SNRTDX和SNRTSX分别表示2颗卫星的信噪比水平，可从头文件中获取。

2.2 SINC估测模型

模型是定量遥感的基础，在极化干涉SAR的发展过程中，陆续提出了多种将雷达观测参数与植被垂直结构联系起来的反演模型，其中最常用的是RVoG模型。该模型包含2层结构，上层为代表森林冠层结构的体散射层，下层为地表层，有效地将γVol与森林垂直结构相联系，且简化了散射过程的复杂性(Treuhaftetal.， 1996)。RVoG模型假设相干函数如式(5)所示：

(5)

式(6)为体散射失相干γV的表达式，当假设森林冠层结构由各项同性的随机散射体组成时，可采用指数函数对f(z)建模，如式(7)所示：

(6)

f(z)=exp(2σz/cosθ)。

(7)

式中：kz为局部垂直波数；hv为森林高度；f(z)为垂直结构函数；σ为体散射层中的消光系数；θ为入射角。

(8)

x∈[0,π],y∈[0,1]；

(9)

x≈π-2sin-1(y0.8)；

(10)

(11)

相干性γObs与局部垂直波数kz估算准确与否影响SINC模型森林高度估测结果的精度。本研究在估算γObs时，先对SLC数据采用3×3窗口多视生成干涉条纹，再以3×3滑动窗口估算相干性。相干性影像在距离-多普勒坐标下的像元大小约为6 m×6 m，经地理编码后，在地理坐标空间下的像元大小与参考DEM相同，即5 m×5 m。由于相干性是在约18 m×18 m面积区域平均得到的，而SLC数据距离向和方位向的分辨率分别为2.6和3.3 m，因此多视视数约为38视。根据Seymour等(1994)提出的Cramer-Rao bound方法，图5绘制了不同视数条件下由式(1)估算γObs时的方差随|γObs|大小变化的情况，可见当多视视数约为40视时，方差总是小于0.02，因此本研究计算得到的γObs与真值γ的差异可以忽略。

图5 基于Cramer-Rao下界的相干性计算误差模拟Fig.5 Simulation of coherence estimation bias based on Cramer-Rao lower bound

在体散射失相干模型中，kz将f(z)映射到γV，决定着相位和相干性对植被高度的敏感程度，同时也是SINC模型中关键的缩放因子。kz的表达式见式(12)：

(12)

式中：Bn为垂直基线；R为SAR系统与地物间的距离；hoa为模糊高，代表2π相位变化对应的垂直高度；θ0和θi分别为中心入射角和局部入射角，θi需要由参考DEM计算。

本研究分别以SRTM DEM、LiDAR DEM作为参考DEM，得到的kz分别如图6a、b所示。可以看到，二者整体趋势相同，其中LiDAR DEM生成的kz对细节刻画更清晰，而SRTM DEM生成的kz则较为粗糙。图6c为二者的统计直方图，其分布很相近。

2.3 DSM-DEM差分法

DSM-DEM差分法采用的DSM由TanDEM-X干涉相位生成，DEM为高精度的LiDAR DEM。从DSM中减去DEM，得到有效散射相位中心高度hphase，由于X-波段的微波对森林具有一定的穿透性，hphase往往低于森林高度hv，因此需要根据一部分实测数据对hphase进行校正才能获得准确的森林高度hv。校正模型采用式(13)所示的线性回归方程：

hv=B+A×hphase。

(13)

式中：B和A为回归方程系数，其值采用75对(hv，hphase)数据作为训练样本通过最小二乘法拟合得到。

训练样本从研究区中均匀选取，其空间分布见图7，红色三角形代表训练样本中心点位置。

2.4 森林高度估测精度评价方法

为了对以上2种森林高度估测模型或方法的精度进行评价，在研究区内均匀布设150个检验样本，其空间分布见图7，黄色圆点代表检验样本(与差分法训练样本独立)中心点位置，以该点为中心设一个取样窗口，窗口内所有有效像元的平均值为该样本的取值。样本的待检验值分别从SINC模型法、DSM-DEM差分法估测结果中提取，样本的实测值(参考值)从LiDAR CHM中获取。精度评价指标包括决定系数(R2)、均方根误差(root mean square error, RMSE)和估测精度(estimation accuracy, EA)，计算公式如下：

图7 训练及检验样本中心点位置Fig.7 Location of the center of the samples for training and validating

(14)

(15)

(16)

以15 m×15 m、30 m×30 m、50 m×50 m和100 m×100 m为样本取样窗口大小进行精度评价。

3 结果与分析

3.1 干涉相干性处理结果

研究区的相干系数和相干性密度图如图8a、b所示，相干系数集中分布在0.9附近，干涉数据质量较好。图8c为检验样本相干系数与LiDAR 冠层高度的散点图，随着冠层高度增加，相干性降低，表明相干系数具有一定的森林高度估测潜力。

图8 干涉相干性处理结果Fig.8 InSAR coherence processing results

3.2 森林高度估测结果

图9a为LiDAR CHM，图9b为DSM-DEM差分法得到的森林高度估测结果，图9c、d分别为以LiDAR DEM和SRTM DEM为参考DEM的SINC模型法估测结果。可以看到，3种森林高度估测结果与LiDAR CHM具有很好的一致性，其中DSM-DEM差分法的森林高度分布与LiDAR CHM在细节上最为接近，体现出该方法的有效性，但DSM-DEM差分法与LiDAR CHM相比森林高度估测结果存在低估现象，表明X-波段的微波对本研究区林分具有一定穿透性。图9c、d的森林高度分布相似，说明2种DEM对SINC模型估测结果的影响较小。

图9 森林高度参考和InSAR估测结果Fig.9 Forest height as reference and that estimated with InSAR data

3.3 森林高度估测精度评价

为了定量评价森林高度估测结果，同时分析样本尺度对估测精度的影响，利用均匀选取的样本，分别以15 m×15 m、30 m×30 m、50 m×50 m和100 m×100 m为样本取样窗口大小进行平均得到森林高度，对DSM-DEM差分法和SINC模型法的估测结果进行精度评价，分析高分辨率与中等分辨率参考DEM对SINC模型估测精度的影响。

3.3.1 DSM-DEM差分法精度评价 DSM-DEM差分法精度评价结果如图10所示，其中左侧为采用式(13)校正前的结果，右侧为校正后的结果。可见，DSM-DEM差分法森林高度估测结果与LiDAR CHM之间具有良好的相关性，且随样本尺度增大估测精度稳定提高，当样本大小从15 m×15 m变化到100 m×100 m时，R2从0.57增加到0.79，RMSE由3.38 m降至2.57 m，精度由67.62%提高到75.44%。在15 m×15 m样本尺度上，校正对DSM-DEM差分法估测精度的改善效果并不明显，甚至有些下降，从校正前的67.62%降至校正后的61.79%，原因可能是该尺度下地理编码和相位上的误差对估测结果影响较大。样本尺度逐渐增大，校正后估测结果的精度大幅度提高，当样本大小为100 m×100 m时，总精度由75.44%提高到84.41%。

图10 差分法精度评价(左侧为校正前，右侧为校正后，从上到下样本采样窗口依次增大)Fig.10 Accuracy assessment of phase difference method before (left) and after (right) calibration with increased sampling window size from up to down panels

3.3.2 SINC模型法精度评价图11左侧为采用LiDAR DEM进行kz估算和地理编码的SINC模型估测结果，右侧为采用SRTM DEM的估测结果。可以看到，二者在各样本尺度上精度均很接近，且随样本尺度增加估测精度逐渐提高。当样本大小为15 m×15 m时，估测结果的精度较差，RMSE在4 m左右，存在个别误差较大的样本，考虑到LiDAR数据和SAR数据获取时间间隔不到1个月，可能是由于地理编码误差引起的。随样本尺度增大，地理编码误差对估测结果的影响减弱，精度逐渐提高。当样本大小为100 m×100 m时，2种SINC模型估测结果的R2分别为0.54、0.51，RMSE分别为2.38、2.51 m，总精度已达到DSM-DEM差分法校正前的水平，分别为77.19%、75.99%。

图11 SINC模型法精度评价(左、右侧分别以LiDAR DEM、SRTM DEM为参考DEM)Fig.11 Accuracy assessment of SINC model with LiDAR DEM (left) and SRTM DEM (right)

3.3.3 不同分辨率DEM对SINC模型估测精度的影响分析图12为基于SINC模型采用2种参考DEM(LiDAR DEM和SRTM DEM)估测森林高度结果之间的相关性。随样本尺度增大，相关系数逐渐从15 m×15 m时的0.70增加到100 m×100 m时的0.96，说明高分辨率和中等分辨率DEM对SINC模型的影响较小，尤其当样本尺度较大时，其影响可忽略不计。该现象可从2个方面解释：首先，参考DEM在SINC模型森林高度估测中的作用是估算坡度，结合成像几何计算得到kz，而研究区坡度变化缓慢，由LiDAR DEM和SRTM DEM估算的坡度差异并不明显，因此最终结果差异较小；另外，随样本尺度增大，地理编码误差对估测精度的影响也会随平均像元增加而降低，进一步减轻2种DEM对SINC模型森林高度估测的影响。

图12 SINC模型估测结果在不同样本尺度下的相关性Fig.12 Correlation of SINC model estimation results at different sample sizes

3.4 大范围森林高度制图结果

基于相干系数的SINC模型估测森林高度可取得较高精度，且高分辨率和中等分辨率DEM对估测结果影响较小。本研究以SRTM DEM为参考DEM，采用SINC模型制作整景TanDEM-X影像覆盖范围的森林高度分布图(图13a)，可以看到，不同高度林分在影像上具有很好的空间变异性(图中黑色方框区域为上述模型评价采用的试验区范围)。

图13 采用SINC模型制作的森林高度分布图及与Google Earth多光谱遥感影像的对比Fig.13 Forest height distribution map using SINC model and comparison with Google Earth multispectral remote sensing image

图13b、c分别为从图13a中选取的典型区域与Google Earth多光谱遥感影像的对比。图13b中，一条河流分布在多光谱影像左侧，而在对应的森林高度分布图中，估测的高度则有较大误差(森林高度应为0 m)，分析其原因在于水体的相干性很低，在SINC模型中被认为是具有较强体散射的植被，进而得到错误的高度信息，因此在大范围森林高度制图时，对水体区域进行掩膜很有必要。图13c中，多光谱影像右下角的林分存在间伐现象，在SINC模型估测的森林高度分布图中也可清楚反映出来，体现了该方法的有效性。

4 讨论

本研究采用TanDEM-X单极化InSAR数据，分别以基于相位信息的DSM-DEM差分法和基于相干系数的SINC模型法估测森林高度，并分析LiDAR DEM和SRTM DEM对SINC模型估测精度和稳定性的影响。结果表明，与DSM-DEM差分法相比，尽管SINC模型法估测精度略有下降，但其既不需要实测森林高度数据进行标定，也不需要输入高分辨率的DEM，具有大范围森林高度制图的潜力和更大的实际应用价值。同时，本研究的相关内容仍存在一定的不足之处和值得继续研究的方向：

1) 使用DSM-DEM差分法估测森林高度时，需要对有效散射中心高度进行校正。本研究采用线性方程进行校正，得到估测森林高度；但实际情况下，微波对较低森林的穿透能力要强于较高森林，即有效散射中心高度与森林高度并不严格符合线性关系。因此，经线性校正后，可能会出现校正不足或过校正的现象。

2) 在SINC模型中，kz是体散射失相干反演森林高度的关键参数，其准确计算和取值范围对估测精度有一定影响。图14模拟了不同森林高度下相干性与kz的关系。一方面，kz计算准确与否直接关系到森林高度估测结果误差大小，如当相干性γVol=0.8时，不同kz可分别对应森林高度10～50 m范围。在本研究区，地形起伏缓慢，30 m分辨率的SRTM DEM可以很好刻画坡度，计算的kz较为准确，但在地形变化迅速区域，SRTM DEM对地形的刻画可能不够细致，由kz计算误差引起的森林高度误差会增大，需要更精细的DEM来计算kz。另一方面，kz取决于垂直基线长度、波长、斜距和入射角等参数，其变化范围对SINC模型估测精度影响较大。当kz过小时，相干性对森林高度变化的敏感性降低，如当kz=0.05时，对于最大高度为50 m的森林，相干性变化范围仅为0.75～1；而当kz过大时，模糊高降低，对于超过模糊高的高度将引起低估，如当kz=0.30时，对于高度超过20 m的森林，相干性降为0。因此，基于SINC模型估测森林高度时，应根据研究区森林高度变动范围选择具有合理干涉成像参数的数据，以减少由于kz不合适导致的误差。

图14 kz敏感性分析Fig.14 kz-sensitivity analysis

3) 非体散射失相干校正是SINC模型法关键的一步。在本研究中，TanDEM-X InSAR数据已进行距离向频谱滤波和精确配准，且时间基线可忽略，需要考虑的非体散射失相干只有信噪比失相干γSNR。假设γSNR=0.95，图15分析了γSNR对SINC模型估测结果的影响。图15a表示不同kz下相干性与森林高度的关系，可以看到，在森林高度较低时，γSNR校正后相干性水平明显提高，而对于较高的森林，相干性变化则不明显。图15b表示不同森林高度下由γSNR引起的高度估测误差，同样说明非体散射失相干对较低森林的影响强于较高森林。本研究在校正γSNR时，采用的γSNR为整景影像均值，实际上不同地物及不同高度森林的γSNR并不相同，因此校正后的森林高度可能在一些区域仍存在高估现象，而在另一些区域则低估。