面向森林高度提取的光学多角度立体观测影像模拟

姚远， 倪文俭， 张志玉

1. 中国科学院空天信息创新研究院 遥感科学国家重点实验室, 北京 100101;2. 中国科学院大学 资源与环境学院, 北京 100049

1 引 言

森林生态系统作为陆地生态系统的主体，在地球系统碳循环过程中发挥着重要的作用。森林植被碳储量直接取决于森林冠层高度、密度和类型。因此，森林冠层高度作为3 个决定变量之一，是森林植被碳储量估算的重要参数（Lei 等，2021）。高精度森林冠层高度制图对于陆地生态系统碳收支核算研究具有重要意义（Drake等，2002；Potapov等，2021）。

借助于各种不同的测距技术，通过对森林垂直结构的直接探测，实现高精度的森林冠层高度制图，是国际学术研究前沿（Duncanson 等，2010；Khati等，2019；Qi等，2019），如激光雷达LiDAR（Light Detection and Ranging）和雷达干涉InSAR（Interferometric Synthetic Aperture Radar）分别使用的是激光测距技术与微波测距技术。激光雷达能够实现森林垂直结构的直接测量，但现有激光雷达数据（GEDI、ICESat-2）通常为点/线采样，缺乏空间连续性。可用于森林高度提取的微波测距技术包括雷达干涉，雷达极化干涉和雷达干涉层析（InSAR，PolInSAR 和TomoSAR），但受时间去相干影响，且能够进行大区域森林高度制图的雷达干涉数据源有限（Lei 等，2021）。光学多角度立体观测（即摄影测量）依靠视觉测距技术，通过对不同观测角度引起的视差的测量提取森林冠层顶部的高程，结合林下地形，也具备对森林冠层垂直结构进行直接探测的潜力（Ni 等，2014；梁顺林 等，2020）。近十几年来，星载摄影测量系统不断发展，光学立体数据源不断增加，目前已有多个星载摄影测量系统具备全球覆盖能力，如ASTER 于1999 年开始立体观测，ALOS/PRISM 在2006 年至2011 年获取了全球尺度的立体影像，随后陆续发射了ZY-3，WorldView 系列等。通过多星多时相数据融合，具备区域高精度森林冠层高度制图及森林高度变化监测的能力。国内外学者利用星载光学多角度立体观测影像针对森林冠层高度提取开展了大量的研究工作（Hobi 和Ginzler，2012；Ni 等，2012；St-Onge 等，2008；Tian 等，2017）。根据所用影像获取模式的不同，这些研究可以分为两类，即沿轨前后视立体观测和侧视立体观测。

沿轨前后视立体观测影像通常由星载摄影测量系统同时搭载多个相机来获取，如日本的ALOS/PRISM 和中国的资源三号（ZY-3）等。受专业化数据处理技术的限制，早期星载摄影测量数据的使用主要集中在测绘领域。随着计算机视觉和数字摄影测量的发展，数据处理自动化程度不断提高，利用星载摄影测量数据进行森林冠层高度提取研究日益受到重视。如Avtar 和Sawada（2013）的研究区位于柬埔寨约1340 km2的森林地区，研究结果表明，ALOS/PRISM 立体像对得到的数字表面模型DSM（Digital Surface Model）与SRTM-DEM的差值和与ASTER-GDEM 的差值所估算的森林高度变化范围分别为-29 m—-20 m 和-29 m—-10 m；Ni 等（2014） 的研究区位于美国缅因州约35×35 km2的森林地区，研究表明，ALOS/PRISM 三线阵影像多角度融合的点云能够增强对森林冠层结构的刻画能力，所提取的森林高度与机载大光斑激光雷达（LVIS）获取的中心高（RH50）的相关性（r）可达0.86，RMSE为2.57 m。

侧视立体观测影像通常由高分辨率卫星通过相机侧摆来获取，如IKONOS、WorldView-1/2 和中国的高分二号（GF-2）等。Neigh 等（2014）基于IKONOS 立体数据生成的DSM，结合美国国家高程数据集NED（National Elevation Dataset）获取冠层高度模型CHM（Canopy Height Model），利用机载激光雷达G-LiHT 获取的CHM 进行精度评价，结果表明，两者具有较好的一致性（R2＞0.7，RMSE＜3 m）。倪文俭等（2018）的研究区位于内蒙古大兴安岭面积约19.1 km×14.5 km 的林区，利用交会角约20°的GF-2立体数据提取的的森林高度，在空间分布格局上与机载激光雷达获取的森林高度具有较好一致性，两者相关性（r）为0.71，RMSE为3.6 m。Persson和Perko（2016）利用WorldView-2 立体像对生成瑞典南部约1200 ha 森林的DSM，结合机载激光扫描系统ALS（Airborne Scanning System）获取的数字地形模型DTM（Digital Terrain Model）估算森林高度，结果表明，估算的森林高度与激光雷达参考数据之间相关系数（r）＞0.9，表明WorldView-2 图像的立体匹配适用于森林高度制图。

上述研究证明了利用光学多角度立体观测影像进行森林高度提取的潜力。光学遥感影像容易受到云、雨等天气因素的影响，多星多时相数据融合是实现区域高精度森林高度制图的关键，但不同卫星获取的影像分辨率和观测角度差异较大，如ASTER 后视观测角度为27.6°，影像分辨率为15 m；ALOS/PRISM 的前、后视观测角度为±24°，影像分辨率仅2.5 m；ZY-3 的01 星前、后视角度为±22°，正视影像分辨率为2.1 m，前、后视影像分辨率为3.5 m；WorldView 系列卫星通过快速改变姿态，可进行多角度斜侧视观测，并提供亚米级影像；GF-2影像分辨率达0.8 m，且具备±35°的侧摆能力。利用多角度立体观测影像进行视差测量的基础是基于影像纹理特征所识别的同名点，而同名点密度和视差的测量精度受立体像对间影像纹理特征的影响较大（Montesano 等，2017），因为对于给定的森林场景，影像纹理特征主要受空间分辨率、观测角度等影响。因此，系统研究影像分辨率和观测角度变化对多角度立体观测点云垂直分布的影响规律，是实现基于多星多角度立体观测影像区域森林冠层高度制图的基础。森林光学多角度立体观测模拟模型是开展该研究的有效理论工具。针对这一需求，Ni 等（2019）研发了面向森林高度提取的“森林光学多角度立体观测模型（LandStereo）”。LandStereo 模型能够在景观尺度上精细刻画森林冠层沿轨前后视多角度立体特征，但该模型尚不具备侧视立体观测影像的模拟能力。因此，本文在原有模型的基础上，发展了侧视立体成像模拟方法，使LandStereo 模型具备在任意指定观测角度的影像模拟能力。

2 立体观测模型的改进

LandStereo 模型采用计算机图像仿真工具POV-Ray（Persistence of Vision Raytracer）进行图像的渲染，通过自身的场景描述语言SDL（Scene Description Language）进行复杂场景的定义和渲染。模型根据输入的地形和单木结构信息，首先利用SDL 语言进行地理三维场景的构建，然后根据给定的光照与观测几何关系和相机成像参数，严密控制相机运动，借助于光线追踪算法进行线阵推扫影像模拟。为便于使用常规摄影测量处理软件对模拟图像进行分析，模型也建立了对应的有理多项式函数模型RFM （Rational Function Model）（Tao 和Hu，2001），依据输入场景的地理范围和成像几何模型，建立空间格网，采用与地形无关的方式迭代求解有理多项式系数RPC（Rational Polynomial Coefficient）。

原模型仅能模拟沿轨前后视线阵推扫影像，为使其具备任意指定观测角度影像的模拟能力，本文在沿用LandStereo 模型的原有结构的基础上，一方面改进了线阵摄影中心坐标（即推扫过程中的相机位置）的计算方式，由原来只考虑观测高度角的变化，改进为同时考虑观测方位角与高度角变化，以形成任意指定角度的图像模拟能力；另一方面，为了获取用于计算RPC 的地面控制点GCP（Ground Control Point），改进了严密成像几何模型，由原来只考虑垂直于轨道窄视场成像情况，改为考虑垂直于轨道宽视场成像情况，以形成任意指定角度图像RPC的计算能力。

2.1 线阵摄影中心坐标的计算

设卫星轨道高度为H，观测方位角为α，其以正北方向作为起始方向（即0°），沿顺时针逐渐增大，观测高度角为β，模拟区域起始扫描行中心A坐标为（X0，Y0，Z0）。图1展示了太阳—传感器—观测目标相对位置及传感器推扫过程。其中，以推扫方向为y轴，以垂轨方向为x轴，以铅垂方向为z轴。传感器和太阳的位置是相对于模拟区域扫描行中心定义的，αs和βs表示太阳方位角和高度角，Si表示传感器推扫至第i行时的摄影中心坐标。

图1 太阳—传感器—观测目标几何关系Fig.1 Sun-sensor-target geometry

则 第i行 影 像 的 摄 影 中 心 坐 标Si（Xs，i，Ys，i，Zs，i）计算公式如下：

式中，Δx和Δy分别为x和y方向的地面分辨率。

因此，通过对观测方位角和高度角的定义，即可实现不同的立体观测方式，如：

（1）当观测高度角β为90°时，为正视观测。

（2）当观测方位角α为0°或180°时，为前后视观测。

（3）当观测方位角α为90°或270°时，为正侧视观测。

（4）当观测方位角α与高度角β为其他角度组合时，为斜侧视观测。

2.2 严密成像几何模型

设相机焦距为f0，相机探测单元（CCD）尺寸为Δc，相机CCD数目为mc。图2展示了成像范围内任意点的地理坐标（Xi，Yi，Zi）与其对应的图像坐标（xi，yi）的严密成像几何关系，其中，以沿轨方向为y轴，以垂轨方向为x轴，以铅垂方向为z轴。图2（a）表示y向坐标求解，A（X0，Y0，Z0）为初始行中心坐标，S0和Si为传感器在初始行和推扫至第i行的摄影中心坐标，θ1和θ2分别为视线沿轨 （俯仰角）和垂轨方向的角度分量（翻滚角），yi表示观测目标点（Xi，Yi，Zi）在地面的投影点的y坐标。图2（b）表示x向坐标求解，z′轴表示在θ2平面内与铅垂方向相差θ1的方向，摄影中心Si的坐标（Xs，i，Ys，i，Zs，i），f为相机沿轨方向的等效焦距，O为地主点，x0表示扫描行中心所对应的图像坐标，其对应CCD阵列的中心位置，xi表示观测目标点（Xi，Yi，Zi） 在地面的投影点的x坐标。

图2 严密成像几何模型Fig.2 Rigorous imaging geometric model

按照上述构建的严密成像几何模型，（Xi，Yi，Zi）与（xi，yi）的转换关系可表示为

式中，θ1和θ2与观测方位角和高度角间的转换关系为

由于模拟图像为像素坐标系（以行、列号定义像素的位置），假设任意点的图像坐标（xi，yi）对应像素坐标（ri，ci），则两者转换关系如式（7）所示。

式中，Δx和Δy与前述定义一致，分别为x和y方向的地面分辨率，x0为模拟区域扫描行中心根据式（4）求得的x坐标，f为相机沿轨方向的等效焦距，分别表达为

利用构建的严密成像几何模型，采用与地形无关的方式生成空间格网，迭代求解任意指定角度图像RPC。在构建的有理多项式函数模型中，地理坐标（Xi，Yi，Zi）与像点坐标（ri，ci）的数学关系为（Ni等，2019）：

式中，式（11）中的（Xi，Yi，Zi）和（ri，ci）分别为归一化后的地理坐标与像点坐标。式（12）中的Pn（n=1，2，3，4）表示从（Xi，Yi，Zi）到（ri，ci）转换的多项式，ai，n（i=1，2，3，…，20）为RPC系数。

3 多角度图像模拟与处理分析

改进后的LandStereo 模型理论上具备了对任意指定角度观测图像的模拟能力，通过模拟实例对此进行分析验证。

3.1 输入参数与模拟图像

LandStereo 模型所需的输入包括相机参数、观测几何参数和三维场景构建参数，如表1所示。相机参数包括焦距、CCD 单元尺寸、CCD 单元数目；观测几何参数包括太阳方位角和高度角、观测方位角和高度角、卫星轨道高度、模拟场景的起始扫描行中心坐标和模拟的图像行数；三维场景构建参数包括用于构建三维地形的数字地形模型DTM 相关参数和地物颜色信息。输入DTM（图3）位于内蒙古自治区呼伦贝尔市根河市大兴安岭林区（121°25′E-121°32′E，50°53′N-50°58′N），空间参考为WGS_1984_UTM_zone_51N，图3 中红框所示区域为模拟实例的最小模拟范围，具体参数如表1所示。目前常见的森林生长模型（如ZELIG，JABOWA-FORET，SIBBORK 等）能够模拟复杂森林场景并提供单木结构信息，本文以白桦为例，共获取230386 条单木数据，每条数据具体给出每木的位置、胸径、树高、冠长和冠幅等信息，部分示例如表2所示。

表1 多角度图像模拟参数Table 1 Parameters of simulation of multi-view images

表2 单木结构参数示例Table 2 Examples of Structural parameters of individual tree

图3 输入的数字地形模型Fig.3 Input DTM

为了充分展示模型对任意指定角度图像模拟能力，除了表1给出的热点方向外，另外模拟了观测方位角为160°，高度角分别为30°、45°、75°、90°的4景影像，用以分析相同观测方位角不同高度角的差异；模拟了观测高度角为60°，方位角分别为70°、250°、340°的3景影像，用以分析相同观测高度角不同方位角的差异。以上角度组合如图4（a）所示，其中三角形号表示观测目标位置，五角星号表示太阳照射角度，黑点表示不同传感器的观测角度。图4（b）展示了根据表1的输入参数模拟的热点方向图像。图4（c）—（j）为图4（b）中红点（1）稀疏林区不同观测角度模拟图像的局部放大图像，所对应的观测角度已在图4（a）中用相应的字母标出。为了充分分析不同角度观测模拟图像的差异，将各局部图像的分辨率进一步提高。图4（c）—（g）展示了观测方位角为160°，高度角分别为30°，45°，60°，75°，90°的模拟影像，即相同观测方位角，不同观测高度角；图4（e）的观测方向与太阳照射方向相同，即为热点方向，影像中只有光照树冠和光照背景两个分量；图4（c），（d）观测高度角低于热点方向，影像中含有光照树冠、光照背景和阴影背景3个分量，且由于图4（c）的观测高度角低于4（d），因此，能看到更多的阴影背景；图4（f），（g）观测高度角高于热点方向，影像中可以看到光照树冠、光照背景、阴影背景和阴影树冠四个分量，由于图4（g）的观测高度角高于4（f），且都高于热点方向，因此，图4（g）中能看到更多的阴影背景；此外可以看到，随着观测高度角的变化，光照树冠与阴影的相对位置关系在发生变化，且变化只发生在视线方向。

图4（h）—（j）展示了观测高度角为60°，方位角分别为70°，250°，340°的模拟影像，与图4（e）共同组成了相同观测高度角，不同方位角的情况。图4（k）—（n）分别展示了图4（h），（e），（i），（j）中红框所示区域的进一步放大的图像。由于太阳照射方向没有发生变化，因此，阴影背景位置不会发生变化，而光照树冠的投影方向会随观测方位角发生变化。红线十字交叉处标注的是相同的地理位置，对比图4（k）—（n）可以看到，阴影背景与红线十字的相对位置没有发生变化，但光照树冠和阴影背景的相对位置随着观测方位角的变化而在方向上发生变化。特别是图4（k）和（m），图4（l）和（n）的观测方位角之间相差180°，因此，图4（k）和（m），图4（l）和（n）的光照树冠投影方向相反。图4（o）—（q）为图4（b）中红点（2）浓密森林区域的局部放大图像，所对应的观测角度分别与图4（j），（g），（e）相同。可以看出，浓密林区光照树冠与阴影的相对位置关系随观测角度的变化规律与稀疏林区一致。这些模拟结果与理论预期相一致，表明实际模拟的任意指定观测角度图像能够精确表达所定义的模型输入。

3.2 模型的验证与分析

为了对改进的模型进行验证，模拟了正视观测和观测高度角为75°，方位角分别为0°（后视），90°（正侧视），225°（斜侧视）的4 景裸露地表（无单木覆盖）影像。根据表1 的图像模拟输入参数及成像几何模型，4景模拟影像的空间分辨率及覆盖范围如表3所示，其中正视模拟图像覆盖范围最小为3.84 km×6 km，如图3 中的红框所示。其他角度模拟影像的覆盖范围在正视基础上进行扩展。

表3 多角度模拟图像分辨率及覆盖范围Table 3 Spatial resolution and coverage of multi-view images

以正视观测影像为左片，分别以其他角度为右片构成立体像对，同时需结合各观测角度模拟影像的RPC。在RPC 求解过程中，利用构建的严密成像几何模型，在模拟覆盖范围内分别生成3000（30 行×20 列×5 层）个控制点和24000（60行×40列×10层）个检查点，经最小二乘多次迭代，最终达到收敛。表4展示了以均方根误差和最大误差两项指标依次评价4 景模拟图像RPC 的解算精度。从表4可以看出，利用控制点和检查点拟合的4景图像RPC精度很高，x、y向的最大误差和均方根误差均小于0.1 个像素，表明求解的RPC 精度满足立体匹配要求，也验证了严密成像几何模型的准确性。因此，任意角度组合构成的立体像对可使用ENVI、PCI 等通用软件进行立体匹配处理，可得到数字表面模型（DSM）。

表4 RPC求解精度Table 4 Accuracy of RPC solution

针对裸露地表，通过将提取的地表高程与输入DTM 对比来对模型进行验证。以正视与斜侧视（α=225°，β=75°） 观测图像的立体匹配为例，图5展示了裸露地表提取精度验证结果，图5（a）展示了地表高程提取结果，图5（b）表示提取的地表高程与输入DTM 的散点分析结果，其他角度模拟影像组合提取的裸地高程精度如表5所示。

通过图5（a）与输入DTM（图3）的对比可以看出，提取的裸露地表DSM 的地形起伏特征与输入DTM 相同。从图5（b）和表5 可以看出，不同立体像对提取的DSM 与输入DTM 具有明显的线性相关关系，两者相关性（r）达0.99，均方根误差均小于2 m。结果表明，不同观测角度裸露地表构成的立体像对提取的地表高程值与输入DTM 具有很好的一致性，证明了改进的LandStereo 模型在几何上的准确性。

图5 裸露地表高程提取精度验证Fig.5 Accuracy verification of extracted elevation of bare ground

表5 多角度图像裸露地表高程提取精度Table 5 Extraction accuracy of bare ground elevation in multi-view images

3.3 基于不同角度模拟影像组合的森林冠层高度提取

在改进的模型能够提取高精度的地表高程的基础上，为了分析不同观测方向对森林高度提取精度的影响，采用与裸露地表相同的观测方式和模拟输入参数，进行山地森林场景的图像模拟。图6 展示了模拟结果（以斜侧视为例）。通过图6与图4（b）的对比可以看出，两景模拟图像的整体纹理特征有所差异，这是由于模拟图像间不同的观测方向导致的。

图6 斜侧视观测模拟结果Fig.6 Results of Simulated image of （α=225°， β=75°）

由于山地森林模拟图像的观测方式和成像几何与裸露地表相同，基于模拟的多角度山地森林图像，采用与裸露地表观测相同的RPC 和立体像对组合方式进行匹配处理，将提取的DSM 与输入DTM 做差值得到CHM。图7 展示了森林高度提取对比结果。图7（a）为根据输入的单木列表得到的CHM 制图，图7（b）展示了提取的森林冠层高度制图。从图7（a）、（b）可以看出，提取的森林冠层高度整体的空间分布格局与参考CHM较为相似。

图7 森林高度提取结果与CHM的对比Fig.7 Comparison between extracted forest height and CHM

根据输入的单木列表得到的CHM 作为参考数据，通过在30 m×30 m 的林分尺度上取平均值，将提取的森林冠层高度与参考CHM 进行精度评价。图8分别展示了不同角度模拟影像组合提取的森林高度的精度。图8（a）—（c）分别表示正视与正侧视（α=90°，β=75°），后视（α=0°，β=75°）和斜侧视（α=225°，β=75°）观测模拟的山地森林影像提取的森林高度精度， 它们的观测高度角相同，但方位角不同，即展示了不同观测方位角对森林高度提取精度的影响。

从图8（a）—（c）可以看出，不同角度模拟影像组合能够有效提取森林冠层高度。森林高度的提取精度随着观测方位角的变化而发生变化，如图8（b）所示，正视与后视提取的森林高度与参考值之间线性相关，其相关系数（r）为0.94，RMSE 为2.4 m。而图8（a），（c）所示的正视与正侧视，正视与斜侧视提取的森林高度与参考值之间具有较大的散点分布差异。由表3可知，构成不同像对组合的模拟影像分辨率近似相等。对于给定的森林场景，由不同观测角度引起的图像空间分辨率的变化对森林高度提取精度影响很小，森林高度提取精度可能与不同的观测角度引起的森林结构差异与纹理变化有关。具体的原因和影响规律需要做进一步的深入分析。现有结果可以表明，观测角度组合的变化是利用立体像对提取森林冠层高度精度的重要影响因素。

图8 提取的森林高度与参考CHM的散点结果Fig.8 Scattering plots of extracted forest height with reference CHM

4 结 论

本文在原有“森林光学多角度立体观测模型（LandStereo）” 仅具备前后视多角度立体观测模拟能力的基础上，通过对相机观测几何严密控制算法的扩展，发展了侧视立体成像模拟方法，使LandStereo 模型具备林区任意指定观测角度的图像模拟能力。

利用改进的模型对相同观测方位角、不同高度角和相同观测高度角、不同方位角的多角度影像进行了模拟分析，结果表明，随着观测角度的变化，光照树冠、光照背景、阴影背景和阴影树冠四个分量的比例变化与理论预期相一致。表明改进的LandStereo 模型模拟的任意指定观测角度图像能够精确表达所定义的模型输入。进而，利用改进的模型模拟了正视观测和观测高度角为75°，方位角分别为0°（后视），90°（正侧视），225°（斜侧视）的裸露地表影像，用以验证模拟图像的精度，结果表明，改进的LandStereo 模型提取的地表高程与输入DTM 具有很好的一致性，证明了改进的LandStereo 模型在几何上的准确性。在此基础上，采用与裸露地表相同的观测方式和模拟参数，进行了山地森林图像的模拟与森林冠层高度的提取，并结合参考数据进行了分析。结果表明，对于给定的森林场景，不同角度组合的立体像对对森林冠层高度的提取精度存在一定的差异，表明观测角度组合的变化是利用立体像对提取森林冠层高度精度的重要影响因素。

在后续的研究中，将模拟各种观测几何条件下的多角度立体影像，系统研究不同成像条件（观测角度、分辨率等）对森林高度提取精度的影响规律。