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

面向森林高度提取的光学多角度立体观测影像模拟 Simulation of multi-view stereo optical imagery for extraction of forest canopy height 1 2 first-author 姚远，研究方向为森林遥感。E-mail： yaoyuan19@mails.ucas.ac.cn 姚远，研究方向为森林遥感。E-mail： yaoyuan19@mails.ucas.ac.cn yaoyuan19@mails.ucas.ac.cn 1 2 corresp 倪文俭，研究方向为森林遥感机理模型发展及基于多源遥感数据的森林植被碳储量及其动态变化监测。E-mail：niwj@aircas.ac.cn 倪文俭，研究方向为森林遥感机理模型发展及基于多源遥感数据的森林植被碳储量及其动态变化监测。E-mail：niwj@aircas.ac.cn niwj@aircas.ac.cn 1 中国科学院空天信息创新研究院 遥感科学国家重点实验室, 北京 100101 中国科学院空天信息创新研究院 遥感科学国家重点实验室, 北京 100101 State Key Laboratory of Remote Sensing Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100101, China State Key Laboratory of Remote Sensing Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100101, China 中国科学院大学 资源与环境学院, 北京 100049 中国科学院大学 资源与环境学院, 北京 100049 College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China 森林冠层高度是森林生态系统植被碳储量的重要指示因子之一。已有研究表明光学多角度立体观测（即摄影测量）数据具备森林冠层高度提取的潜力，然而光学影像容易受到云、雨等天气因素的影响，区域高精度森林冠层高度制图必需依靠多星多时相数据融合才能实现。现有星载摄影测量系统获取的影像分辨率和观测角度差异较大，系统研究影像分辨率和观测角度变化对多角度立体观测点云垂直分布的影响规律，是实现多星多时相数据融合的前提。森林光学多角度立体观测模型是开展该研究的有效理论工具，但现有“森林光学多角度立体观测模型（LandStereo）”仅具备沿轨前后视立体观测影像模拟能力，无法用于WorldView等侧视立体观测影像的模拟分析。因此，本文对LandStereo模型进行了改进，使模型具备林区任意指定观测角度的影像模拟能力，进而利用改进后的模型对正视、正侧视和斜侧视等典型观测角度影像进行了模拟，并对不同观测角度组合构成的立体像对森林高度的提取精度进行了分析。结果表明，改进后的模型能够用于刻画森林冠层多角度立体特征，观测角度的变化是利用立体像对提取森林冠层高度精度的重要影响因素。本文为系统研究林区任意角度组合的多角度立体观测点云随影像分辨率和观测角度的变化规律提供了理论工具。 Canopy height is one of important indicators of carbon storage in forest ecosystems. Previous studies have demonstrated that optical stereo imagery held great potential for deriving forest canopy height. However, optical images are easily affected by cloud coverage and rains. The regional mapping of forest canopy height has to be achieved by the synthesis of multi-sensor and multi-temporal imagery. The spatial resolutions and viewing angles of existing spaceborne stereoscopic systems are quite different. It is essential to make a systematic investigation about the impact of image spatial resolutions and viewing angles on the vertical distribution of stereoscopic point clouds, which is the basis for the synergy of images acquired at different seasons even by different cameras. The theoretical model for the simulation of optical stereo imagery over forested areas is necessary for such studies. The LandStereo is such kind of model, which can simulate the along-track(only In-track viewing) stereoscopic imagery like ALOS/PRISM and ZY-3. However, the current version of LandStereo has no mode to simulate the along-track(including In-track and Cross-track) stereoscopic imagery like Worldview-1/2/3. Therefore, this study reported the modification of the LandStereo model to have it being able to simulate images acquired by any possible viewing directions in forest areas. Firstly, the calculation method of linear array projection center coordinates is improved, from originally only considering the change of observation altitude angle to considering the change of azimuth angle and altitude angle. Secondly, Rigorous imaging geometric model is improved to obtain RPC of images acquired by any possible viewing directions. Based on the improved LandStereo model, the bare and mountainous forest images with altitude angle of 75° and azimuth angle of 0°, 90° , and 225° are simulated to verify the accuracy and extract forest canopy height. The surface elevation extracted by the improved LandStereo model is consistent with the input DTM with high accuracy ( r =0.99, RMSE=1.507), which proves the geometric accuracy of the improved LandStereo model. There are certain differences in the extraction accuracy of forest canopy height between stereo images of different angles. The results showed that the modified model could correctly simulate the stereoscopic features of forest canopy with given view direction, also initially demonstrated that view angle was an important factor affecting the estimation accuracy of forest canopy height by stereoscopic images. 森林高度 多角度 立体观测 影像模拟 LandStereo POV-Ray 植被光学模型 forest height multi-view stereoscopic observations imagery simulation LandStereo POV-Ray vegetation optical model 1　引 言 1 森林生态系统作为陆地生态系统的主体，在地球系统碳循环过程中发挥着重要的作用。森林植被碳储量直接取决于森林冠层高度、密度和类型。因此，森林冠层高度作为3个决定变量之一，是森林植被碳储量估算的重要参数（ Lei等，2021 Lei等，2021 ）。高精度森林冠层高度制图对于陆地生态系统碳收支核算研究具有重要意义（ Drake等，2002 Drake等，2002 ； Potapov等，2021 Potapov等，2021 ）。 借助于各种不同的测距技术，通过对森林垂直结构的直接探测，实现高精度的森林冠层高度制图，是国际学术研究前沿（ Duncanson等，2010 Duncanson等，2010 ； Khati等，2019 Khati等，2019 ； Qi等，2019 Qi等，2019 ），如激光雷达LiDAR（Light Detection and Ranging）和雷达干涉InSAR（Interferometric Synthetic Aperture Radar）分别使用的是激光测距技术与微波测距技术。激光雷达能够实现森林垂直结构的直接测量，但现有激光雷达数据（GEDI、ICESat-2）通常为点/线采样，缺乏空间连续性。可用于森林高度提取的微波测距技术包括雷达干涉，雷达极化干涉和雷达干涉层析（InSAR，PolInSAR和TomoSAR），但受时间去相干影响，且能够进行大区域森林高度制图的雷达干涉数据源有限（ Lei等，2021 Lei等，2021 ）。光学多角度立体观测（即摄影测量）依靠视觉测距技术，通过对不同观测角度引起的视差的测量提取森林冠层顶部的高程，结合林下地形，也具备对森林冠层垂直结构进行直接探测的潜力（ Ni等，2014 Ni等，2014 ； 梁顺林 等，2020 梁顺林 等，2020 ）。近十几年来，星载摄影测量系统不断发展，光学立体数据源不断增加，目前已有多个星载摄影测量系统具备全球覆盖能力，如ASTER于1999年开始立体观测，ALOS/PRISM在2006年至2011年获取了全球尺度的立体影像，随后陆续发射了ZY-3，WorldView系列等。通过多星多时相数据融合，具备区域高精度森林冠层高度制图及森林高度变化监测的能力。国内外学者利用星载光学多角度立体观测影像针对森林冠层高度提取开展了大量的研究工作（ Hobi和Ginzler，2012 Hobi和Ginzler，2012 ； Ni等，2012 Ni等，2012 ； St-Onge等，2008 St-Onge等，2008 ； Tian等，2017 Tian等，2017 ）。根据所用影像获取模式的不同，这些研究可以分为两类，即沿轨前后视立体观测和侧视立体观测。 沿轨前后视立体观测影像通常由星载摄影测量系统同时搭载多个相机来获取，如日本的ALOS/PRISM和中国的资源三号（ZY-3）等。受专业化数据处理技术的限制，早期星载摄影测量数据的使用主要集中在测绘领域。随着计算机视觉和数字摄影测量的发展，数据处理自动化程度不断提高，利用星载摄影测量数据进行森林冠层高度提取研究日益受到重视。如 Avtar和Sawada（2013） Avtar和Sawada（2013） 的研究区位于柬埔寨约1340 km 2 的森林地区，研究结果表明，ALOS/PRISM立体像对得到的数字表面模型DSM（Digital Surface Model）与SRTM-DEM的差值和与ASTER-GDEM的差值所估算的森林高度变化范围分别为-29 m—-20 m和-29 m—-10 m； Ni等（2014） Ni等（2014） 的研究区位于美国缅因州约35×35 km 2 的森林地区，研究表明，ALOS/PRISM三线阵影像多角度融合的点云能够增强对森林冠层结构的刻画能力，所提取的森林高度与机载大光斑激光雷达（LVIS）获取的中心高（RH50）的相关性（ r ）可达0.86，RMSE为2.57 m。 侧视立体观测影像通常由高分辨率卫星通过相机侧摆来获取，如IKONOS、WorldView-1/2和中国的高分二号（GF-2）等。 Neigh等（2014） Neigh等（2014） 基于IKONOS立体数据生成的DSM，结合美国国家高程数据集NED（National Elevation Dataset）获取冠层高度模型CHM（Canopy Height Model），利用机载激光雷达G-LiHT获取的CHM进行精度评价，结果表明，两者具有较好的一致性（ R 2 > 0.7 ， RMSE < 3 m）。 倪文俭等（2018） 倪文俭等（2018） 的研究区位于内蒙古大兴安岭面积约19.1 km×14.5 km的林区，利用交会角约20°的GF-2立体数据提取的的森林高度，在空间分布格局上与机载激光雷达获取的森林高度具有较好一致性，两者相关性（ r ）为0.71，RMSE为3.6 m。 Persson和Perko（2016） 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 Montesano等，2017 ），因为对于给定的森林场景，影像纹理特征主要受空间分辨率、观测角度等影响。因此，系统研究影像分辨率和观测角度变化对多角度立体观测点云垂直分布的影响规律，是实现基于多星多角度立体观测影像区域森林冠层高度制图的基础。森林光学多角度立体观测模拟模型是开展该研究的有效理论工具。针对这一需求， Ni等（2019） Ni等（2019） 研发了面向森林高度提取的“森林光学多角度立体观测模型（LandStereo）”。LandStereo模型能够在景观尺度上精细刻画森林冠层沿轨前后视多角度立体特征，但该模型尚不具备侧视立体观测影像的模拟能力。因此，本文在原有模型的基础上，发展了侧视立体成像模拟方法，使LandStereo模型具备在任意指定观测角度的影像模拟能力。 2　立体观测模型的改进 1 LandStereo模型采用计算机图像仿真工具POV-Ray（Persistence of Vision Raytracer）进行图像的渲染，通过自身的场景描述语言SDL（Scene Description Language）进行复杂场景的定义和渲染。模型根据输入的地形和单木结构信息，首先利用SDL语言进行地理三维场景的构建，然后根据给定的光照与观测几何关系和相机成像参数，严密控制相机运动，借助于光线追踪算法进行线阵推扫影像模拟。为便于使用常规摄影测量处理软件对模拟图像进行分析，模型也建立了对应的有理多项式函数模型RFM（Rational Function Model）（ Tao和Hu，2001 Tao和Hu，2001 ），依据输入场景的地理范围和成像几何模型，建立空间格网，采用与地形无关的方式迭代求解有理多项式系数RPC（Rational Polynomial Coefficient）。 原模型仅能模拟沿轨前后视线阵推扫影像，为使其具备任意指定观测角度影像的模拟能力，本文在沿用LandStereo模型的原有结构的基础上，一方面改进了线阵摄影中心坐标（即推扫过程中的相机位置）的计算方式，由原来只考虑观测高度角的变化，改进为同时考虑观测方位角与高度角变化，以形成任意指定角度的图像模拟能力；另一方面，为了获取用于计算RPC的地面控制点GCP（Ground Control Point），改进了严密成像几何模型，由原来只考虑垂直于轨道窄视场成像情况，改为考虑垂直于轨道宽视场成像情况，以形成任意指定角度图像RPC的计算能力。 2.1　线阵摄影中心坐标的计算 2 设卫星轨道高度为 H ，观测方位角为 α， 其以正北方向作为起始方向（即0°），沿顺时针逐渐增大，观测高度角为 β ，模拟区域起始扫描行中心 A 坐标为（ X 0 ，Y 0 ，Z 0 ）。 图1 图1 展示了太阳—传感器—观测目标相对位置及传感器推扫过程。其中，以推扫方向为 y 轴，以垂轨方向为 x 轴，以铅垂方向为 z 轴。传感器和太阳的位置是相对于模拟区域扫描行中心定义的 ，α s 和 β s 表示太阳方位角和高度角， S i 表示传感器推扫至第 i 行时的摄影中心坐标。 则第 i 行影像的摄影中心坐标 S i （ X s ，i ，Y s ，i ，Z s ，i ）计算公式如下： <math display="block" id="M1"><msub><mrow><mi>X</mi></mrow><mrow><mi mathvariant="normal">s</mi><mo>,</mo><mi mathvariant="normal">i</mi></mrow></msub><mo>=</mo><msub><mrow><mi>X</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>+</mo><mfenced separators="|"><mrow><mi>H</mi><mo>-</mo><msub><mrow><mi>Z</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow></mfenced><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">t</mi><mi>β</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi>α</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331506&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331500&type= 43.51866913 4.40266657 <math display="block" id="M2"><msub><mrow><mi>Y</mi></mrow><mrow><mi mathvariant="normal">s</mi><mo>,</mo><mi mathvariant="normal">i</mi></mrow></msub><mo>=</mo><msub><mrow><mi>Y</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo>+</mo><mfenced separators="|"><mrow><mi>H</mi><mo>-</mo><msub><mrow><mi>Z</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow></mfenced><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">t</mi><mi>β</mi><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mi>α</mi><mo>+</mo><mfenced separators="|"><mrow><mi>i</mi><mo>-</mo><mn mathvariant="normal">1</mn></mrow></mfenced><msub><mrow><mo>∆</mo></mrow><mrow><mi>y</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331518&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331502&type= 61.55266953 4.40266657 <math display="block" id="M3"><msub><mrow><mi>Z</mi></mrow><mrow><mi mathvariant="normal">s</mi><mo>,</mo><mi>i</mi></mrow></msub><mo>=</mo><mi>H</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331510&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331509&type= 11.43000031 3.72533321 式中， Δ x 和 Δ y 分别为 x 和 y 方向的地面分辨率。 因此，通过对观测方位角和高度角的定义，即可实现不同的立体观测方式，如： （1）当观测高度角 β 为90°时，为正视观测。 （2）当观测方位角 α 为0°或180°时，为前后视观测。 （3）当观测方位角 α 为90°或270°时，为正侧视观测。 （4）当观测方位角 α 与高度角 β 为其他角度组合时，为斜侧视观测。 2.2　严密成像几何模型 2 设相机焦距为 f 0 ，相机探测单元（CCD）尺寸为 Δ c ，相机CCD数目为 m c 。 图2 图2 展示了成像范围内任意点的地理坐标（ X i ，Y i ，Z i ）与其对应的图像坐标（ x i ，y i ）的严密成像几何关系，其中，以沿轨方向为 y 轴，以垂轨方向为 x 轴，以铅垂方向为 z 轴。 图2（a） 图2（a） 表示 y 向坐标求解， A （ X 0 ，Y 0 ，Z 0 ）为初始行中心坐标， S 0 和 S i 为传感器在初始行和推扫至第 i 行的摄影中心坐标， θ 1 和 θ 2 分别为视线沿轨 （俯仰角）和垂轨方向的角度分量（翻滚角）， y i 表示观测目标点（ X i ，Y i ，Z i ）在地面的投影点的 y 坐标。 图2（b） 图2（b） 表示 x 向坐标求解， z′ 轴表示在 θ 2 平面内与铅垂方向相差 θ 1 的方向，摄影中心 S i 的坐标（ X s ，i ，Y s ，i ，Z s ，i ）， f 为相机沿轨方向的等效焦距， O 为地主点， x 0 表示扫描行中心所对应的图像坐标，其对应CCD阵列的中心位置， x i 表示观测目标点（ X i ，Y i ，Z i ） 在地面的投影点的 x 坐标。 1,2; 按照上述构建的严密成像几何模型，（ X i ，Y i ，Z i ）与（ x i ，y i ）的转换关系可表示为 <math display="block" id="M4"><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mfenced separators="|"><mrow><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>-</mo><msub><mrow><mi>X</mi></mrow><mrow><mi>s</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow></mfenced><mo>-</mo><mfrac><mrow><mo stretchy="false">(</mo><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>-</mo><msub><mrow><mi>Z</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo stretchy="false">)</mo><mi mathvariant="normal">t</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">n</mi><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mrow><mrow><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mrow><mi>y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mfenced separators="|"><mrow><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>-</mo><msub><mrow><mi>Y</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow></mfenced><mo>+</mo><mfenced separators="|"><mrow><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>-</mo><msub><mrow><mi>Z</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow></mfenced><mi mathvariant="normal">t</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">n</mi><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331549&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331534&type= 50.20733261 15.49400043 式中， θ 1 和 θ 2 与观测方位角和高度角间的转换关系为 <math display="block" id="M5"><mi mathvariant="normal">t</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">n</mi><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>=</mo><mo>-</mo><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">t</mi><mi>β</mi><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mi>α</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331553&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331551&type= 27.00866890 4.48733330 <math display="block" id="M6"><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>=</mo><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><mi>β</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi>α</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331557&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331555&type= 24.72266769 4.48733330 由于模拟图像为像素坐标系（以行、列号定义像素的位置），假设任意点的图像坐标（ x i ，y i ）对应像素坐标（ r i ， c i ），则两者转换关系如 式（7） 式（7） 所示。 <math display="block" id="M7"><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><msub><mrow><mi>c</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mrow><mi>m</mi></mrow><mrow><mi>c</mi></mrow></msub></mrow><mrow><mn mathvariant="normal">2</mn></mrow></mfrac><mo>+</mo><mfrac><mrow><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>-</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow><mrow><msub><mrow><mo>∆</mo></mrow><mrow><mi>x</mi></mrow></msub></mrow></mfrac></mtd></mtr><mtr><mtd><msub><mrow><mi>r</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mrow><mi>y</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow><mrow><msub><mrow><mo>∆</mo></mrow><mrow><mi>y</mi></mrow></msub></mrow></mfrac></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331538&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331537&type= 27.85533524 19.38866615 式中， Δ x 和 Δ y 与前述定义一致，分别为 x 和 y 方向的地面分辨率， x 0 为模拟区域扫描行中心根据 式（4） 式（4） 求得的 x 坐标， f 为相机沿轨方向的等效焦距，分别表达为 <math display="block" id="M8"><msub><mrow><mo>∆</mo></mrow><mrow><mi>x</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mrow><mo>∆</mo></mrow><mrow><mi>c</mi></mrow></msub><mfenced separators="|"><mrow><mi>H</mi><mo>-</mo><msub><mrow><mi>Z</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub></mrow></mfenced></mrow><mrow><mi>f</mi><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331540&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331561&type= 25.31533241 9.05933285 <math display="block" id="M9"><msub><mrow><mo>∆</mo></mrow><mrow><mi>y</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mrow><mo>∆</mo></mrow><mrow><mi>c</mi></mrow></msub><mo stretchy="false">(</mo><mi>H</mi><mo>-</mo><msub><mrow><mi>Z</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mo stretchy="false">)</mo></mrow><mrow><mi>f</mi><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><msup><mrow><mi mathvariant="normal">s</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331544&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331542&type= 25.56933212 8.63599968 <math display="block" id="M10"><mi>f</mi><mo>=</mo><mfrac><mrow><msub><mrow><mi>f</mi></mrow><mrow><mn mathvariant="normal">0</mn></mrow></msub><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi>β</mi></mrow><mrow><mi mathvariant="normal">c</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">s</mi><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331564&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331547&type= 15.57866764 7.95866632 利用构建的严密成像几何模型，采用与地形无关的方式生成空间格网，迭代求解任意指定角度图像 RPC 。在构建的有理多项式函数模型中，地理坐标（ X i ，Y i ，Z i ）与像点坐标（ r i ，c i ）的数学关系为（ Ni等，2019 Ni等，2019 ）： <math display="block" id="M11"><msub><mrow><mi>r</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mrow><mi>p</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mfenced separators="|"><mrow><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mfenced></mrow><mrow><msub><mrow><mi>p</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mfenced separators="|"><mrow><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mfenced></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331579&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331566&type= 24.46866608 9.90600014 , <math display="block" id="M12"><msub><mrow><mi>c</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mrow><mi>p</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub><mfenced separators="|"><mrow><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mfenced></mrow><mrow><msub><mrow><mi>p</mi></mrow><mrow><mn mathvariant="normal">4</mn></mrow></msub><mfenced separators="|"><mrow><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mfenced></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331581&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331567&type= 24.72266769 9.90600014 <math display="block" id="M13"><msub><mrow><mi>p</mi></mrow><mrow><mi>n</mi></mrow></msub><mo>=</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">1</mn><mo>,</mo><mi>n</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">2</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">3</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">4</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">5</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">6</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">7</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">8</mn><mo>,</mo><mi>n</mi></mrow></msub><msubsup><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">9</mn><mo>,</mo><mi>n</mi></mrow></msub><msubsup><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">10</mn><mo>,</mo><mi>n</mi></mrow></msub><msubsup><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">11</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">12</mn><mo>,</mo><mi>n</mi></mrow></msub><msubsup><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msubsup><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">13</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><msubsup><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">14</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow></msub><msubsup><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">15</mn><mo>,</mo><mi>n</mi></mrow></msub><msubsup><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">16</mn><mo>,</mo><mi>n</mi></mrow></msub><msubsup><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msubsup><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">17</mn><mo>,</mo><mi>n</mi></mrow></msub><msub><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow></msub><msubsup><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">18</mn><mo>,</mo><mi>n</mi></mrow></msub><msubsup><mrow><mi>X</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">19</mn><mo>,</mo><mi>n</mi></mrow></msub><msubsup><mrow><mi>Y</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup><msub><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mi>a</mi></mrow><mrow><mn mathvariant="normal">20</mn><mo>,</mo><mi>n</mi></mrow></msub><msubsup><mrow><mi>Z</mi></mrow><mrow><mi>i</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331584&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=43331569&type= 157.90333557 9.31333351 式中， 式（11） 式（11） 中的（ X i ，Y i ，Z i ）和（ r i ， c i ）分别为归一化后的地理坐标与像点坐标。 式（12） 式（12） 中的 P n （ n =1，2，3，4）表示从（ X i ，Y i ，Z i ）到（ r i ， c i ）转换的多项式， a i，n （ i =1，2，3，…，20）为 RPC 系数。 3　多角度图像模拟与处理分析 1 改进后的LandStereo模型理论上具备了对任意指定角度观测图像的模拟能力，通过模拟实例对此进行分析验证。 3.1　输入参数与模拟图像 2 LandStereo模型所需的输入包括相机参数、观测几何参数和三维场景构建参数，如 表1 表1 所示。相机参数包括焦距、CCD单元尺寸、CCD单元数目；观测几何参数包括太阳方位角和高度角、观测方位角和高度角、卫星轨道高度、模拟场景的起始扫描行中心坐标和模拟的图像行数；三维场景构建参数包括用于构建三维地形的数字地形模型DTM相关参数和地物颜色信息。输入DTM（ 图3 图3 ）位于内蒙古自治区呼伦贝尔市根河市大兴安岭林区（121°25'E-121°32'E，50°53'N-50°58'N），空间参考为WGS_1984_UTM_zone_51N， 图3 图3 中红框所示区域为模拟实例的最小模拟范围，具体参数如 表1 表1 所示。目前常见的森林生长模型（如ZELIG，JABOWA-FORET，SIBBORK等）能够模拟复杂森林场景并提供单木结构信息，本文以白桦为例，共获取230386条单木数据，每条数据具体给出每木的位置、胸径、树高、冠长和冠幅等信息，部分示例如 表2 表2 所示。 序号 参数 参数值 1 DTM列数 8814 2 DTM行数 9560 3 DTM- x 向分辨率/m 1 4 DTM- y 向分辨率/m 1 5 DTM最大值/m 1163.0 6 DTM最小值/m 760.0 7 DTM左上 x 坐标 389052 8 DTM左上 y 坐标 5648265 9 焦距/mm 998.432 10 卫星轨道高度/km 500 11 CCD单元数目 3840 12 CCD单元尺寸/μm 2.0 13 模拟起始点 x 坐标 395190 14 模拟起始点 y 坐标 5640550 15 模拟起始点 z 坐标 784 16 模拟行数 6000 17 太阳高度角/（°） 60 18 太阳方位角/（°） 160 19 观测高度角/（°） 60 20 观测方位角/（°） 160 21 树冠RGB值 （0.358，0.451，0.262） 22 土壤RGB值 （0.545，0.466，0.348） 23 树干RGB值 （0.346，0.329，0.330） 序号 x y z 胸径/cm 树高/m 冠长/m 冠幅/m 1 390893.243 5639716.789 869.718 11.2 13.045 3.780 2.587 2 390895.233 5639429.976 839.466 9.4 11.054 3.165 3.071 3 390895.667 5638714.744 790.480 6.6 8.796 2.933 2.452 4 390902.842 5648070.270 968.068 16.8 17.969 5.341 3.496 5 390900.118 5647907.498 951.397 15.9 18.268 4.547 5.436 6 389495.046 5641788.931 1045.328 23.8 25.374 6.924 6.232 7 395077.012 5646967.111 956.958 16.2 17.375 4.460 4.152 8 395789.378 5642531.206 883.766 12.0 13.395 3.222 3.610 9 395787.670 5642191.759 800.150 7.2 9.340 1.756 2.856 10 396108.745 5644761.545 857.865 10.5 12.514 2.680 3.661 11 396112.980 5642888.376 906.344 13.3 14.881 4.029 4.015 12 394918.374 5646023.969 1034.196 20.6 22.929 5.307 6.159 … 为了充分展示模型对任意指定角度图像模拟能力，除了 表1 表1 给出的热点方向外，另外模拟了观测方位角为160°，高度角分别为30°、45°、75°、90°的4景影像，用以分析相同观测方位角不同高度角的差异；模拟了观测高度角为60°，方位角分别为70°、250°、340°的3景影像，用以分析相同观测高度角不同方位角的差异。以上角度组合如 图4（a） 图4（a） 所示，其中三角形号表示观测目标位置，五角星号表示太阳照射角度，黑点表示不同传感器的观测角度。 图4（b） 图4（b） 展示了根据 表1 表1 的输入参数模拟的热点方向图像。 图4（c） 图4（c） —（j）为 图4（b） 图4（b） 中红点（1）稀疏林区不同观测角度模拟图像的局部放大图像，所对应的观测角度已在 图4（a） 图4（a） 中用相应的字母标出。为了充分分析不同角度观测模拟图像的差异，将各局部图像的分辨率进一步提高。 图4（c） 图4（c） —（g）展示了观测方位角为160°，高度角分别为30°，45°，60°，75°，90°的模拟影像，即相同观测方位角，不同观测高度角； 图4（e） 图4（e） 的观测方向与太阳照射方向相同，即为热点方向，影像中只有光照树冠和光照背景两个分量； 图4（c） 图4（c） ，（d）观测高度角低于热点方向，影像中含有光照树冠、光照背景和阴影背景3个分量，且由于 图4（c） 图4（c） 的观测高度角低于4（d），因此，能看到更多的阴影背景； 图4（f） 图4（f） ，（g）观测高度角高于热点方向，影像中可以看到光照树冠、光照背景、阴影背景和阴影树冠四个分量，由于 图4（g） 图4（g） 的观测高度角高于4（f），且都高于热点方向，因此， 图4（g） 图4（g） 中能看到更多的阴影背景；此外可以看到，随着观测高度角的变化，光照树冠与阴影的相对位置关系在发生变化，且变化只发生在视线方向。 1,2;3,4,5,6;7,8,9,10;11,12,13,14;15,16,17; 图4（h） 图4（h） —（j）展示了观测高度角为60°，方位角分别为70°，250°，340°的模拟影像，与 图4（e） 图4（e） 共同组成了相同观测高度角，不同方位角的情况。 图4（k） 图4（k） —（n）分别展示了 图4（h） 图4（h） ，（e），（i），（j）中红框所示区域的进一步放大的图像。由于太阳照射方向没有发生变化，因此，阴影背景位置不会发生变化，而光照树冠的投影方向会随观测方位角发生变化。红线十字交叉处标注的是相同的地理位置，对比 图4（k） 图4（k） —（n）可以看到，阴影背景与红线十字的相对位置没有发生变化，但光照树冠和阴影背景的相对位置随着观测方位角的变化而在方向上发生变化。特别是 图4（k） 图4（k） 和（m）， 图4（l） 图4（l） 和（n）的观测方位角之间相差180°，因此， 图4（k） 图4（k） 和（m）， 图4（l） 图4（l） 和（n）的光照树冠投影方向相反。 图4（o） 图4（o） —（q）为 图4（b） 图4（b） 中红点（2）浓密森林区域的局部放大图像，所对应的观测角度分别与 图4（j） 图4（j） ，（g），（e）相同。可以看出，浓密林区光照树冠与阴影的相对位置关系随观测角度的变化规律与稀疏林区一致。这些模拟结果与理论预期相一致，表明实际模拟的任意指定观测角度图像能够精确表达所定义的模型输入。 3.2　模型的验证与分析 2 为了对改进的模型进行验证，模拟了正视观测和观测高度角为75°，方位角分别为0°（后视），90°（正侧视），225°（斜侧视）的4景裸露地表（无单木覆盖）影像。根据 表1 表1 的图像模拟输入参数及成像几何模型，4景模拟影像的空间分辨率及覆盖范围如 表3 表3 所示，其中正视模拟图像覆盖范围最小为3.84 km×6 km，如 图3 图3 中的红框所示。其他角度模拟影像的覆盖范围在正视基础上进行扩展。 2 观测角度 2 空间分辨率/m 2 覆盖范围/m x 向 y 向 x 向 y 向 正视 1 1 3840 6000 （ α= 90° ， β= 75°） 1.035 1.035 3974.4 6209.9 （ α= 0° ， β= 75°） 1.035 1.072 3974.4 6432 （ α= 225° ， β= 75°） 1.035 1.054 3974.4 6324 以正视观测影像为左片，分别以其他角度为右片构成立体像对，同时需结合各观测角度模拟影像的RPC。在RPC求解过程中，利用构建的严密成像几何模型，在模拟覆盖范围内分别生成3000（30行×20列×5层）个控制点和24000（60行×40列×10层）个检查点，经最小二乘多次迭代，最终达到收敛。 表4 表4 展示了以均方根误差和最大误差两项指标依次评价4景模拟图像RPC的解算精度。从 表4 表4 可以看出，利用控制点和检查点拟合的4景图像RPC精度很高， x 、 y 向的最大误差和均方根误差均小于0.1个像素，表明求解的RPC精度满足立体匹配要求，也验证了严密成像几何模型的准确性。因此，任意角度组合构成的立体像对可使用ENVI、PCI等通用软件进行立体匹配处理，可得到数字表面模型（DSM）。 3 观测角度 4 最大误差/像素 4 RMSE/像素 2 控制点 2 检查点 2 控制点 2 检查点 x y x y x y x y 正视 0.039 0.072 0.042 0.073 0.021 0.033 0.021 0.033 （ α= 0°， β= 75°） 0.044 0.067 0.045 0.083 0.019 0.029 0.019 0.034 （ α =90°， β =75°） 0.0475 0.06 0.045 0.075 0.0195 0.0305 0.0195 0.0315 （ α =225°， β =75°） 0.046 0.0545 0.046 0.06 0.0195 0.03 0.0205 0.0305 针对裸露地表，通过将提取的地表高程与输入DTM对比来对模型进行验证。以正视与斜侧视（ α =225°， β =75°） 观测图像的立体匹配为例， 图5 图5 展示了裸露地表提取精度验证结果， 图5（a） 图5（a） 展示了地表高程提取结果， 图5（b） 图5（b） 表示提取的地表高程与输入DTM的散点分析结果，其他角度模拟影像组合提取的裸地高程精度如 表5 表5 所示。 1,2; 立体像对 线性拟合公式 r RMSE/m 正视 & （ α= 0°， β= 75°） y =1.004 x -3.38869 0.99956 1.79 正视 & （ α= 90°， β= 75°） y =1.0038 x -3.5932 0.99913 1.62 正视 & （ α= 225°， β= 75°） y =1.0042 x -2.1821 0.99925 1.51 通过 图5（a） 图5（a） 与输入DTM（ 图3 图3 ）的对比可以看出，提取的裸露地表DSM的地形起伏特征与输入DTM相同。从 图5（b） 图5（b） 和 表5 表5 可以看出，不同立体像对提取的DSM与输入DTM具有明显的线性相关关系，两者相关性（ r ）达0.99，均方根误差均小于2 m。结果表明，不同观测角度裸露地表构成的立体像对提取的地表高程值与输入DTM具有很好的一致性，证明了改进的LandStereo模型在几何上的准确性。 3.3　基于不同角度模拟影像组合的森林冠层高度提取 2 在改进的模型能够提取高精度的地表高程的基础上，为了分析不同观测方向对森林高度提取精度的影响，采用与裸露地表相同的观测方式和模拟输入参数，进行山地森林场景的图像模拟。 图6 图6 展示了模拟结果（以斜侧视为例）。通过 图6 图6 与 图4（b） 图4（b） 的对比可以看出，两景模拟图像的整体纹理特征有所差异，这是由于模拟图像间不同的观测方向导致的。 由于山地森林模拟图像的观测方式和成像几何与裸露地表相同，基于模拟的多角度山地森林图像，采用与裸露地表观测相同的RPC和立体像对组合方式进行匹配处理，将提取的DSM与输入DTM做差值得到CHM。 图7 图7 展示了森林高度提取对比结果。 图7（a） 图7（a） 为根据输入的单木列表得到的CHM制图， 图7（b） 图7（b） 展示了提取的森林冠层高度制图。从 图7（a） 图7（a） 、（b）可以看出，提取的森林冠层高度整体的空间分布格局与参考CHM较为相似。 1,2;3; 根据输入的单木列表得到的CHM作为参考数据，通过在30 m×30 m的林分尺度上取平均值，将提取的森林冠层高度与参考CHM进行精度评价。 图8 图8 分别展示了不同角度模拟影像组合提取的森林高度的精度。 图8（a） 图8（a） —（c）分别表示正视与正侧视（ α =90°， β =75°），后视（ α =0°， β =75°）和斜侧视（ α =225°， β =75°）观测模拟的山地森林影像提取的森林高度精度， 它们的观测高度角相同，但方位角不同，即展示了不同观测方位角对森林高度提取精度的影响。 1,2,3; 从 图8（a） 图8（a） —（c）可以看出，不同角度模拟影像组合能够有效提取森林冠层高度。森林高度的提取精度随着观测方位角的变化而发生变化，如 图8（b） 图8（b） 所示，正视与后视提取的森林高度与参考值之间线性相关，其相关系数（ r ）为0.94，RMSE为2.4 m。而 图8（a） 图8（a） ，（c）所示的正视与正侧视，正视与斜侧视提取的森林高度与参考值之间具有较大的散点分布差异。由 表3 表3 可知，构成不同像对组合的模拟影像分辨率近似相等。对于给定的森林场景，由不同观测角度引起的图像空间分辨率的变化对森林高度提取精度影响很小，森林高度提取精度可能与不同的观测角度引起的森林结构差异与纹理变化有关。具体的原因和影响规律需要做进一步的深入分析。现有结果可以表明，观测角度组合的变化是利用立体像对提取森林冠层高度精度的重要影响因素。 4　结论 1 本文在原有“森林光学多角度立体观测模型（LandStereo）” 仅具备前后视多角度立体观测模拟能力的基础上，通过对相机观测几何严密控制算法的扩展，发展了侧视立体成像模拟方法，使LandStereo模型具备林区任意指定观测角度的图像模拟能力。 利用改进的模型对相同观测方位角、不同高度角和相同观测高度角、不同方位角的多角度影像进行了模拟分析，结果表明，随着观测角度的变化，光照树冠、光照背景、阴影背景和阴影树冠四个分量的比例变化与理论预期相一致。表明改进的LandStereo模型模拟的任意指定观测角度图像能够精确表达所定义的模型输入。进而，利用改进的模型模拟了正视观测和观测高度角为75°，方位角分别为0°（后视），90°（正侧视），225°（斜侧视）的裸露地表影像，用以验证模拟图像的精度，结果表明，改进的LandStereo模型提取的地表高程与输入DTM具有很好的一致性，证明了改进的LandStereo模型在几何上的准确性。在此基础上，采用与裸露地表相同的观测方式和模拟参数，进行了山地森林图像的模拟与森林冠层高度的提取，并结合参考数据进行了分析。结果表明，对于给定的森林场景，不同角度组合的立体像对对森林冠层高度的提取精度存在一定的差异，表明观测角度组合的变化是利用立体像对提取森林冠层高度精度的重要影响因素。 在后续的研究中，将模拟各种观测几何条件下的多角度立体影像，系统研究不同成像条件（观测角度、分辨率等）对森林高度提取精度的影响规律。 参考文献（References） Avtar R and Sawada H . 2013 . Use of DEM data to monitor height changes due to deforestation . Arabian Journal of Geosciences , 6 ( 12 ): 4859 - 4871 [ DOI: 10.1007/s12517-012-0768-2 http://dx.doi.org/10.1007/s12517-012-0768-2 ] Use of DEM data to monitor height changes due to deforestation Drake J B , Dubayah R O , Clark D B , Knox R G , Blair J B , Hofton M A , Chazdon R L , Weishampel J F and Prince S . 2002 . Estimation of tropical forest structural characteristics using large-footprint lidar . Remote Sensing of Environment , 79 ( 2/3 ): 305 - 319 [ DOI: 10.1016/s0034-4257(01)00281-4 http://dx.doi.org/10.1016/s0034-4257(01)00281-4 ] Estimation of tropical forest structural characteristics using large-footprint lidar Duncanson L I , Niemann K O and Wulder M A . 2010 . Estimating forest canopy height and terrain relief from GLAS waveform metrics . Remote Sensing of Environment , 114 ( 1 ): 138 - 154 [ DOI: 10.1016/j.rse.2009.08.018 http://dx.doi.org/10.1016/j.rse.2009.08.018 ] Estimating forest canopy height and terrain relief from GLAS waveform metrics Hobi M L and Ginzler C . 2012 . Accuracy assessment of digital surface models based on WorldView-2 and ADS80 stereo remote sensing data . Sensors , 12 ( 5 ): 6347 - 6368 [ DOI: 10.3390/s120506347 http://dx.doi.org/10.3390/s120506347 ] Accuracy assessment of digital surface models based on WorldView-2 and ADS80 stereo remote sensing data Khati U , Lavalle M and Singh G . 2019 . Spaceborne tomography of multi-species Indian tropical forests . Remote Sensing of Environment , 229 : 193 - 212 [ DOI: 10.1016/j.rse.2019.04.017 http://dx.doi.org/10.1016/j.rse.2019.04.017 ] Spaceborne tomography of multi-species Indian tropical forests Lei Y , Treuhaft R and Gonçalves F . 2021 . Automated estimation of forest height and underlying topography over a Brazilian tropical forest with single-baseline single-polarization TanDEM-X SAR interferometry . Remote Sensing of Environment , 252 : 112132 [ DOI: 10.1016/j.rse.2020.112132 http://dx.doi.org/10.1016/j.rse.2020.112132 ] Automated estimation of forest height and underlying topography over a Brazilian tropical forest with single-baseline single-polarization TanDEM-X SAR interferometry Liang S L , Bai R , Chen X N , Cheng J , Fan W J , He T , Jia K , Jiang B , Jiang L M , Jiao Z T , Liu Y B , Ni W J , Qiu F , Song L L , Sun L , Tang B H , Wen J G , Wu G P , Xie D H , Yao Y J , Yuan W P , Zhang Y G , Zhang Y Z , Zhang Y T , Zhang X T , Zhao T J and Zhao X . 2020 . Review of China's land surface quantitative remote sensing development in 2019 . Journal of Remote Sensing , 24 ( 6 ): 618 - 671 Review of China's land surface quantitative remote sensing development in 2019 梁顺林 , 白瑞 , 陈晓娜 , 程洁 , 范闻捷 , 何涛 , 贾坤 , 江波 , 蒋玲梅 , 焦子锑 , 刘元波 , 倪文俭 , 邱凤 , 宋柳霖 , 孙林 , 唐伯惠 , 闻建光 , 吴桂平 , 谢东辉 , 姚云军 , 袁文平 , 张永光 , 张玉珍 , 张云腾 , 张晓通 , 赵天杰 , 赵祥 . 2020 . 2019年中国陆表定量遥感发展综述 . 遥感学报 , 24 ( 6 ): 618 - 671 [ DOI: 10.11834/jrs.20209476 http://dx.doi.org/10.11834/jrs.20209476 ] 2019年中国陆表定量遥感发展综述 Montesano P M , Neigh C , Sun G Q , Duncanson L , Van Den Hoek J and Ranson K J . 2017 . The use of sun elevation angle for stereogrammetric boreal forest height in open canopies . Remote Sensing of Environment , 196 : 76 - 88 [ DOI: 10.1016/j.rse.2017.04.024 http://dx.doi.org/10.1016/j.rse.2017.04.024 ] The use of sun elevation angle for stereogrammetric boreal forest height in open canopies Neigh C S R , Masek J G , Bourget P , Cook B , Huang C Q , Rishmawi K and Zhao F . 2014 . Deciphering the precision of stereo IKONOS canopy height models for US forests with G-LiHT Airborne LiDAR . Remote Sensing , 6 ( 3 ): 1762 - 1782 [ DOI: 10.3390/rs6031762 http://dx.doi.org/10.3390/rs6031762 ] Deciphering the precision of stereo IKONOS canopy height models for US forests with G-LiHT Airborne LiDAR Ni W J , Guo Z F , Zhang Z Y , Sun G Q and Huang W L . 2012 . Semi-automatic extraction of digital surface model using ALOS/PRISM data with ENVI // 2012 IEEE International Geoscience and Remote Sensing Symposium . Munich : IEEE: 6557 - 6560 [ DOI: 10.1109/IGARSS.2012.6352097 http://dx.doi.org/10.1109/IGARSS.2012.6352097 ] Semi-automatic extraction of digital surface model using ALOS/PRISM data with ENVI Ni W J , Ranson K J , Zhang Z Y and Sun G Q . 2014 . Features of point clouds synthesized from multi-view ALOS/PRISM data and comparisons with LiDAR data in forested areas . Remote Sensing of Environment , 149 : 47 - 57 [ DOI: 10.1016/j.rse.2014.04.001 http://dx.doi.org/10.1016/j.rse.2014.04.001 ] Features of point clouds synthesized from multi-view ALOS/PRISM data and comparisons with LiDAR data in forested areas Ni W J , Zhang D F , Wang Y , Pang Y , Zhang Z Y , Liu J L , He Y T and Guo W . 2018 . Extraction of forest height by using GF-2 cross-track stereo images . Journal of Remote Sensing , 22 ( 3 ): 392 - 399 Extraction of forest height by using GF-2 cross-track stereo images 倪文俭 , 张大凤 , 汪垚 , 庞勇 , 张志玉 , 刘见礼 , 何亚婷 , 郭伟 . 2018 . 高分二号异轨立体数据的森林高度提取 . 遥感学报 , 22 ( 3 ): 392 - 399 [ DOI: 10.11834/jrs.20187088 http://dx.doi.org/10.11834/jrs.20187088 ] 高分二号异轨立体数据的森林高度提取 Ni W J , Zhang Z Y , Sun G Q and Liu Q H . 2019 . Modeling the stereoscopic features of mountainous forest landscapes for the extraction of forest heights from stereo imagery . Remote Sensing , 11 ( 10 ): 1222 [ DOI: 10.3390/rs11101222 http://dx.doi.org/10.3390/rs11101222 ] Modeling the stereoscopic features of mountainous forest landscapes for the extraction of forest heights from stereo imagery Persson H J and Perko R . 2016 . Assessment of boreal forest height from WorldView-2 satellite stereo images . Remote Sensing Letters , 7 ( 12 ): 1150 - 1159 [ DOI: 10.1080/2150704x.2016.1219424 http://dx.doi.org/10.1080/2150704x.2016.1219424 ] Assessment of boreal forest height from WorldView-2 satellite stereo images Potapov P , Li X Y , Hernandez-Serna A , Tyukavina A , Hansen M C , Kommareddy A , Pickens A , Turubanova S , Tang H , Silva C E , Armston J , Dubayah R , Blair J B and Hofton M . 2021 . Mapping global forest canopy height through integration of GEDI and Landsat data . Remote Sensing of Environment , 253 : 112165 [ DOI: 10.1016/j.rse.2020.112165 http://dx.doi.org/10.1016/j.rse.2020.112165 ] Mapping global forest canopy height through integration of GEDI and Landsat data Qi W L , Lee S K , Hancock S , Luthcke S , Tang H , Armston J and Dubayah R . 2019 . Improved forest height estimation by fusion of simulated GEDI Lidar data and TanDEM-X InSAR data . Remote Sensing of Environment , 221 : 621 - 634 [ DOI: 10.1016/j.rse.2018.11.035 http://dx.doi.org/10.1016/j.rse.2018.11.035 ] Improved forest height estimation by fusion of simulated GEDI Lidar data and TanDEM-X InSAR data St-Onge B , Hu Y and Vega C . 2008 . Mapping the height and above-ground biomass of a mixed forest using lidar and stereo Ikonos images . International Journal of Remote Sensing , 29 ( 5 ): 1277 - 1294 [ DOI: 10.1080/01431160701736505 http://dx.doi.org/10.1080/01431160701736505 ] Mapping the height and above-ground biomass of a mixed forest using lidar and stereo Ikonos images Tao C V and Hu Y . 2001 . A comprehensive study of the rational function model for photogrammetric processing . Photogrammetric Engineering and Remote Sensing , 67 ( 12 ): 1347 - 1357 A comprehensive study of the rational function model for photogrammetric processing Tian J J , Schneider T , Straub C , Kugler F and Reinartz P . 2017 . Exploring digital surface models from nine different sensors for forest monitoring and change detection . Remote Sensing , 9 ( 3 ): 287 [ DOI: 10.3390/rs9030287 http://dx.doi.org/10.3390/rs9030287 ] Exploring digital surface models from nine different sensors for forest monitoring and change detection 10.11834/jrs.20221439 P2 1007-4619（2023）04-1034-11 姚远，倪文俭，张志玉.2023.面向森林高度提取的光学多角度立体观测影像模拟.遥感学报，27（4）： 1034-1044 姚远，倪文俭，张志玉.2023.面向森林高度提取的光学多角度立体观测影像模拟.遥感学报，27（4）： 1034-1044 Yao Y， Ni W J and Zhang Z Y. 2023. Simulation of multi-view stereo optical imagery for extraction of forest canopy height. National Remote Sensing Bulletin， 27（4）：1034-1044 Yao Y， Ni W J and Zhang Z Y. 2023. Simulation of multi-view stereo optical imagery for extraction of forest canopy height. National Remote Sensing Bulletin， 27（4）：1034-1044 战略性国际科技创新合作重点专项（编号：2020YFE0200800）；国家自然科学基金（编号：42090013，42022009） Supported by National Key Research and Development Program of China (No.2020YFE0200800); National Natural Science Foundation of China (No. 42090013, 42022009) 2021-07-01 2023-04-07 2023-05-10 2022-01-12 36985750 遥感学报 4 27