CN103791885A - 面向单幅卫星影像的自获取三角元高度计算方法 - Google Patents

Abstract

本发明提供了一种面向单幅卫星影像的自获取三角元高度计算方法，利用单幅卫星影像中存在的方向信息以及太阳方位角与太阳高度角之间的关系，计算出太阳高度角、太阳方位角和卫星方位角三个角度参数，然后将卫星高度角以建筑物高度和建筑物投影长度进行等效表示，最终求解出不含卫星高度角的三角元高度。针对单幅卫星影像中的建筑物高度求取问题效果良好。

Description

面向单幅卫星影像的自获取三角元高度计算方法

技术领域

本发明涉及一种利用单幅卫星影像信息计算物体高度的方法。

背景技术

近几年来各种高分辨率卫星影像相继出现，为人们利用遥感影像获取建筑物高度提供了可能性。建筑物高度获取在城市规划、城市经济活动以及军事领域中有广泛的应用，在利用卫星影像实现真实城市环境的三维虚拟重建中，建筑物高度获取是为建筑物三维建模提供建模数据的必要途径，决定着三维模型的建模精度，对整个三维虚拟重建环境的真实程度产生直接影响。目前，利用卫星影像获取建筑物高度信息的方法主要有两个研究方向：（1）利用立体像对进行立体量测；（2）从单幅遥感影像入手，利用阴影计算建筑物高度。由于利用立体像对获取高度对两幅遥感影像有着严格的要求，因此在很难在实际中得到应用。而利用阴影从单幅遥感影像中获取高度，则因其方便、经济、操作简单受到了广大研究者的关注。早期的研究简化了三维空间关系，假设太阳、卫星和建筑物位于同一平面，提出卫星成像的二维空间模型，对卫星影像中的建筑物高度进行计算，没有考虑太阳方位角和卫星方位角对阴影的影响。后续研究扩展到三维空间成像模型，进行了多方面的改进：一种是假设遥感影像采样方向为东西、南北方向，即像元的排列是东西、南北方向的，这种情况下没有考虑卫星方位角对阴影测量的影响；一种是根据太阳、卫星与建筑物的相对位置分情况讨论建筑物高度计算模型的方法，这种方法假设“太阳、卫星对屋顶同一点在地面上投影点的连线垂直于建筑物主轴方向（即：建筑物的长边）”，但是多数卫星影像中的图像特征并不满足这一特定情况；另外还有利用垂直于建筑物主轴方向的阴影长度计算建筑物高度的算法，这种算法需要单独量测垂直于建筑物主轴方向的阴影长度和建筑物的方位角，增加了计算的工作量。另一种思路是利用建筑物高度H与阴影长度L之间的线性关系H=K*L，通过已知建筑物的高度和对应阴影长度，反推出K值，进而计算其他建筑物高度，这种方法需要预先知道影像中某一建筑物的实际高度，这一约束条件在实际应用中往往难以达到。

以上利用阴影计算建筑物高度的方法中，除方法“H=K*L”外，从原理上讲，太阳方位角、太阳高度角、卫星方位角和卫星高度角四个参数对高度的求解结果都有影响，建筑物高度计算有解的前提建立在已经提供所有四个角度参数的情况下（一些算法简化地不考虑某些参数的影响，会增大计算误差）。这里称这类方法为“阴影基四角元高度计算法”，其中“四角元”指太阳方位角、太阳高度角、卫星方位角和卫星高度角这四个角度参数。但实际中，这四个角度参数的数据由卫星发射方存放入遥感影像的元数据中，对普通研究者来说很难获取到，这一矛盾严重限制了基于阴影计算建筑物高度的方法的应用和推广。

发明内容

为了克服现有技术的不足，本发明提供一种自获取三角元高度计算法，基于摄影测量理论，首先充分利用单幅卫星影像中存在的方向信息，以及太阳方位角与太阳高度角之间的关系，并辅以星空模拟软件Stellarium计算出太阳高度角、太阳方位角和卫星方位角三个角度参数。因这三个角度参数由该方法自行推导获得，故称之为“自获取三角元”。在此基础上，针对卫星高度角仍无法得到的问题，将卫星高度角以建筑物高度和建筑物投影长度进行等效表示，进一步推导，最终求解出不含卫星高度角的三角元高度。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤一、从卫星影像中测量计算参数：

（1）卫星方位角α：测量卫星影像中建筑物侧棱l在像平面上的投影与正北方向的夹角，即为卫星方位角；以正北方向为起始位置，顺时针方向旋转为正角；

（2）太阳方位角α'：测量卫星影像中建筑物侧棱阴影在像平面上的投影与正北方向的夹角，即为太阳方位角；

（3）建筑物侧棱投影长度L_pb：在像平面上，根据卫星影像中的比例尺获得建筑物侧棱l在像平面上的投影线段长度；

（4）侧棱顶点P及其阴影点的像点间距离L_ps：在像平面上，根据卫星影像中的比例尺获得建筑物侧棱l的顶点P在像平面上的投影点p与顶点P的阴影点在像平面上的投影点s之间的距离；

步骤二、计算建筑物当地当时的太阳高度角；

步骤三、三角元高度计算关系式的推导：

将公式tanλ＝L_h/L_pb代入公式 $L_h=L_{\mathrm{ps}}/\sqrt{\frac{1}{{\tan }^2{\mathrm{\λ}}^{\′}}+\frac{1}{{\tan }^2\mathrm{\λ}}-\frac{2\·\cos (\mathrm{\α}-{\mathrm{\α}}^{\′})}{\tan {\mathrm{\λ}}^{\′}\·(L_h/L_{\mathrm{pb}})}},$ 有：

$L_h=L_{\mathrm{ps}}/\sqrt{\frac{1}{{\tan }^2{\mathrm{\λ}}^{\′}}+\frac{1}{{(L_h/L_{\mathrm{pb}})}^2}-\frac{2\·\cos (\mathrm{\α}-{\mathrm{\α}}^{\′})}{\tan {\mathrm{\λ}}^{\′}\·(L_h/L_{\mathrm{pb}})}}$

其中，λ为卫星高度角，是卫星拍摄方向与像平面的夹角，λ′为太阳高度角，是太阳光线方向与像平面的夹角，L_h为建筑物的高度，

计算得到建筑物高度 $L_h=\tan {\mathrm{\λ}}^{\′}\·(L_{\mathrm{pb}}\·\cos (\mathrm{\α}-{\mathrm{\α}}^{\′})+\sqrt{{L_{\mathrm{pb}}}^2\·{\cos }^2(\mathrm{\α}-{\mathrm{\α}}^{\′})+{L^2}_{\mathrm{ps}}-{L^2}_{\mathrm{pb}}}).$

本发明的有益效果是：针对单幅卫星影像中的建筑物高度求取问题，基于摄影测量学原理，在现有的通过阴影计算高度算法的基础上，解决了实际建筑物求高时，太阳高度角、太阳高方位角、卫星高度角、卫星方位角不易获取的问题，提出了建筑物高度计算新方法。首先通过分析影像中的方向信息，测量计算得到影像中的卫星方位角与太阳方位角大小，然后根据太阳方位角和太阳高度角之间的关系，并借助星空模拟软件Stellrium求解出高度计算中所需的太阳高度角；最后，将卫星高度角对高度计算的影响等效转化为影像中建筑物侧棱投影长度L_pb和建筑物侧棱顶点P及其阴影点的像点间距离L_ps对高度计算的影响，给出了新的高度计算公式。该方法对侧棱高度与建筑物高度一致的平顶建筑物情况，高度计算效果良好。对高杆类物体的高度计算同样适用。

附图说明

图1是太阳、卫星、建筑物阴影之间的空间关系示意图；

图2是本发明方法的流程图；

图中，1–太阳光线方向，2–卫星拍摄方向，3–正北方向，4–建筑物侧棱投影L_pb，5–建筑物侧棱顶点P及其阴影点的像点间距离L_ps，6–卫星高度角λ，7–太阳高度角λ′，8–卫星方位角α，9–太阳方位角α'，10–建筑物侧棱l，11–像平面。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

本发明在未提供太阳方位角、太阳高度角、卫星方位角和卫星高度角数据的情况下，利用单幅卫星影像中建筑物的侧棱及阴影信息计算建筑物的高度。

如附图1所示，该图表示建筑物成像时，太阳、卫星与建筑物的空间几何关系，为了更清楚地表现它们之间的空间关系，正北设置为如图1中箭头所示方向。图中直线段l表示建筑物的一个侧棱，其长度L_h等于建筑物的高度。经过透视投影变换，三维空间中侧棱l的顶点P映射到二维像平面（即投影平面）上的p点，p点称为P点的像点；棱l的下端点映射到像平面上的b点，b点称为棱l下端点的像点；同时P点的阴影映射到像平面上的s点，s点称为P点阴影的像点。卫星方位角、太阳方位角、卫星高度角、和太阳高度角分别以α、α'、λ、λ′表示，假定像平面上，以正北方向为起始位置，顺时针方向旋转的角度为正值，则各角度定义如下：

（1）卫星方位角α：指像平面上，建筑物的侧棱投影L_pb与正北方向的夹角，值域范围为[0,360°]；

（2）太阳方位角α'：指像平面上，建筑物侧棱阴影的投影L_sb与正北方向的夹角，值域范围为[0,360°]；

（3）卫星高度角λ：指三维空间中，卫星拍摄方向与像平面的夹角，值域范围为[0,90°]；

（4）太阳高度角λ′：指三维空间中，太阳光线方向与像平面的夹角，值域范围为[0,90°]

在像平面内，计算p点到b点的长度L_pb、p点到s点的长度L_ps，在△Ppb、△Pbs、△pbs中，分别有：

L_h＝tanλ·L_pb   （1）

L_h＝tanλ′·L_sb   （2）

L² _ps＝L² _pb+L² _sb-2·L_pb·L_sb·cos(α-α')   （3）

将公式（1）和公式（2）带入公式（3）中，经推导可得到建筑物的高度L_h计算关系式如下：

在像平面内计算L_ps，即可计算出高度。

为了使用公式（4）计算建筑物高度，就需要知道卫星方位角、卫星高方位角、太阳方位角和太阳高度角四个角度参数，而这四个角度参数通常封装在图像的元数据中，由卫星发射方提供，对于一般的研究者来说是很难获取到的。

星空模拟软件Stellrium软件可以精确模拟地球上任意一天、任意地点上空各星球的运行情况，包括该星球的类型、赤经/赤纬、方位角/高度角等信息。本发明利用星空模拟软件以及影像中存在的方向信息获取了卫星的方位角、太阳的方位角、太阳的高度角三个参数，并利用这三个参数给出了新的建筑物高度计算公式。具体包括以下步骤：

步骤一、从卫星影像中测量计算部分参数。

（1）卫星方位角α：在图1的像平面上，测量卫星影像中建筑物侧棱l的投影L_pb与正北方向的夹角，即为卫星方位角。以正北方向为起始位置，顺时针方向旋转为正角。

（2）太阳方位角α'：在图1的像平面上，测量卫星影像中建筑物侧棱阴影的投影L_sb与正北方向的夹角，即为太阳方位角。以正北方向为起始位置，顺时针方向旋转为正角。

（3）建筑物侧棱投影长度L_pb：本发明以建筑物侧棱高度代表建筑物高度，L_pb为建筑物侧棱l在像平面上的投影线段长度，也即建筑物投影长度。在图1的像平面上，根据卫星影像中的比例尺工具获得实际长度值。

（4）侧棱顶点P及其阴影点的像点间距离L_ps：在图1的像平面上，L_ps为建筑物侧棱l的顶点P在像平面上的投影点（即像点p）与顶点P的阴影点在像平面上的投影点（即像点s）之间的距离。根据卫星影像中的比例尺工具获得实际长度值。

步骤二、计算建筑物当地当天当时的太阳高度角。

利用步骤一中所确立太阳方位角α'，结合星空模拟软件Stellrium，可以得到当地当天当时的太阳高度角λ′。方法为：通过卫星影像中附带的经纬度、日期信息设定好地点、日期之后，利用已求得的太阳方位角，逐渐调整时间，使得软件中的虚拟太阳方位角与已知的太阳方位角一致，从而得到对应的太阳高度角。

步骤三、三角元高度计算关系式的推导。

在前面两个步骤中，仍无法得到卫星高度角。解决方法为：

由公式（1）有：

tanλ＝L_h/L_pb   （5）

代入公式（4）有：

$L_h=L_{\mathrm{ps}}/\sqrt{\frac{1}{{\tan }^2{\mathrm{\λ}}^{\′}}+\frac{1}{{(L_h/L_{\mathrm{pb}})}^2}-\frac{2\·\cos (\mathrm{\α}-{\mathrm{\α}}^{\′})}{\tan {\mathrm{\λ}}^{\′}\·(L_h/L_{\mathrm{pb}})}}---\left(6\right)$

此处，用L_pb间接反映卫星高度角对高度计算的影响。对上式进行整理得到：

L² _ps＝L² _pb+(L_h/tanλ′)²-2L_pb·(L_h/tanλ′)·cos(α-α')   （7）

通过公式（8）可以得到建筑物高度计算公式：

因此，通过所求得的α、α'、λ'三个角度和测量L_pb、L_ps的长度，利用公式（8）可以计算得到建筑物高度。

对于同一幅卫星影像中的所有建筑物，只需要计算一次共同的卫星方位角、太阳方位角和太阳高度角即可。对不同建筑物测量其L_pb、L_ps，求得各自的高度。

步骤四、误差分析

设建筑物实际高度为h_real，计算高度为h_cal，则：

绝对误差=|h_real-h_cal|，

对于标准的平顶建筑物，一般其侧棱高度与建筑物高度一致，可以得到最佳计算结果。误差的主要来源是卫星影像的分辨率不够高和手工测量产生的误差，卫星影像分辨率越高，计算所得的建筑物高度越接近于实际高度。

如图2所示，本发明包括以下四个步骤：

步骤一、从卫星影像中测量计算各参数，包括：α、α'，L_ps和L_pb。

（1）卫星方位角α

在图1的像平面上，测量卫星影像中建筑物侧棱l的投影L_pb与正北方向的夹角，即为影像拍摄时的卫星方位角。测量角度时，起始边位于正北方向，终止边则是像平面中建筑物侧棱的投影L_pb，顺时针方向旋转为正角。

（2）太阳方位角α'

在图1的像平面上，测量卫星影像中建筑物侧棱l的阴影的投影与正北方向的夹角获得太阳方位角。测量角度时，起始边位于正北方向，终止边则是像平面中建筑物侧棱l阴影的投影L_sb。

（3）建筑物侧棱投影长度L_pb

在图1的像平面上，测量建筑物侧棱l在像平面上的投影线段长度，也即建筑物投影长度。然后，根据卫星影像中的比例尺工具获得实际长度值。

（4）侧棱顶点P及其阴影点的像点间距离L_ps

在图1的像平面上，测量建筑物侧棱顶点P的像点p与顶点P阴影的像点s之间的距离。然后，根据卫星影像中的比例尺工具获得实际长度值。

步骤二、计算建筑物当地当天当时的太阳高度角。

卫星影像中附带有经纬度、日期等信息。在星空模拟软件Stellrium中，首先输入影像中的经纬度信息，设置好具体的地理位置；其次设置影像的拍摄日期，选中星空模拟软件中的虚拟太阳对象后，调节时、分、秒，使得虚拟太阳的方位角等于步骤一中所确立的太阳方位角α'，此时读取相应的虚拟太阳高度角，即为所求的太阳高度角λ′，保存即可。

步骤三、三角元高度计算关系式推导。

在像平面内，计算p点到b点的长度L_pb、p点到s点的长度L_ps，在△Ppb、△Pbs、△pbs中，分别有：

L_h＝tanλ·L_pb   （1）

L_h＝tanλ′·L_sb   （2）

L² _ps＝L² _pb+L² _sb-2·L_pb·L_sb·cos(α-α')   （3）

将公式（1）和公式（2）带入（3）中，经推导可得到建筑物的高度L_h计算关系式如下：

为了解决卫星高度角λ仍无法获取的问题，进一步做如下处理：

由公式（1）得：

tanλ＝L_h/L_pb   （5）

代入公式（4）得：

$L_h=L_{\mathrm{ps}}/\sqrt{\frac{1}{{\tan }^2{\mathrm{\λ}}^{\′}}+\frac{1}{{(L_h/L_{\mathrm{pb}})}^2}-\frac{2\·\cos (\mathrm{\α}-{\mathrm{\α}}^{\′})}{\tan {\mathrm{\λ}}^{\′}\·(L_h/L_{\mathrm{pb}})}}---\left(6\right)$

对上式进行整理得到：

L² _ps＝L² _pb+(L_h/tanλ′)²-2L_pb·(L_h/tanλ′)·cos(α-α')   （7）

通过公式（8）可以得到建筑物高度计算公式：

$L_h=\tan {\mathrm{\λ}}^{\′}\·(L_{\mathrm{pb}}\·\cos (\mathrm{\α}-{\mathrm{\α}}^{\′})+\sqrt{{L_{\mathrm{pb}}}^2\·{\cos }^2(\mathrm{\α}-{\mathrm{\α}}^{\′})+{L^2}_{\mathrm{ps}}-{L^2}_{\mathrm{pb}}})---\left(8\right)$

将步骤一、步骤二中得到的卫星方位角α、太阳方位角α'、太阳高度角λ′、建筑物投影长度L_pb、侧棱顶点及其阴影点的像点间距离L_ps分别带入公式（8）中，可计算得到建筑物高度L_h。

对于同一幅卫星影像中的所有建筑物，只需要计算一次共同的卫星方位角、太阳方位角和太阳高度角即可。对不同建筑物测量其L_pb、L_ps，求得各自的高度。

步骤四、误差分析

设建筑物实际高度为h_real，计算高度为h_cal，则：

绝对误差=|h_real-h_cal|，

误差的主要来源是卫星影像的分辨率不够高和手工测量产生的误差。

实施例1：以2009年10月4日拍摄的上海环球金融中心大厦附近区域的卫星影像为例，计算环球金融中心大厦的高度。

测量影像中建筑物侧棱与正北方向的夹角获得影像拍摄时的卫星方位角α=148.50°。

测量影像中建筑物侧棱阴影的方向与正北方向的夹角获得太阳方位角α'=154.23°。

卫星遥感影像自身提供的比例尺信息为0.498m/像素：

测量影像中建筑物侧棱长度为223个像素，故建筑物投影长度L_pb为：

L_pb=0.498×223=111.054m

测量影像中侧棱顶点P及其阴影点的像点间距离为576个像素，故：

L_ps=0.498×576=286.848m

打开Stellrium软件，设置地点：北纬31°14′16.10″，东经121°30′06.03″；设置日期：2009-10-04；选中虚拟太阳后调整时间，当虚拟太阳的方位角等于前面测量的太阳方位角α'=154.23°时，得到相应的虚拟太阳高度角，此即为所求的太阳高度角λ′，读取得数为λ′=51.33°。

利用公式（8）计算得到环球金融中心大厦高度为496.2349m。结果如表1所示：

表1：环球金融中心大厦数据

实施例2：以2009年6月23日拍摄的北京新保利大厦附近区域的卫星影像为例，计算新保利大厦的高度。

测量影像中建筑物侧棱与正北方向的夹角获得影像拍摄时的卫星方位角α=263.93°。

测量影像中建筑物侧棱阴影的方向与正北方向的夹角获得太阳方位角α'=142.03°。

卫星遥感影像自身提供的比例尺信息为0.222m/像素：

测量影像中建筑物侧棱长度为196个像素，故建筑物投影长度L_pb为：

L_pb=0.222×196=43.512m

测量影像中侧棱顶点P及其阴影点的像点间距离为326个像素，故：

L_ps=0.222×326=72.372m

打开Stellrium软件，设置地点：北纬39°55'52.37"，东经116°25'35.83"；设置日期：2009-6-23；选中虚拟太阳后调整时间，当虚拟太阳的方位角等于前面测量的太阳方位角α'=142.03°时，得到相应的虚拟太阳高度角，此即为所求的太阳高度角λ′，读取得数为λ′=70.08°。

利用公式（8）计算得到新保利大厦高度为108.2837m。结果如如表2所示：

表2：新保利大厦数据

实施例3：以2010年10月27日拍摄的北京水立方游泳馆附近区域卫星影像为例，计算水立方游泳馆高度。

测量影像中建筑物侧棱与正北方向的夹角获得影像拍摄时的卫星方位角α=10.11°。

测量影像中建筑物侧棱阴影的方向与正北方向的夹角获得太阳方位角α'=164.58°。

卫星遥感影像自身提供的比例尺信息为0.234m/像素：

测量影像中建筑物侧棱长度为73个像素，故建筑物投影长度L_pb为：

L_pb=0.234×73=17.082m

测量影像中侧棱顶点P及其阴影点的像点间距离为252个像素，故：

L_ps=0.234×252=58.968

打开Stellrium软件，设置地点：北纬39°59'19.36"，东经116°22'50.85"；设置时间：2010-10-27；选中虚拟太阳后调整时间，当虚拟太阳的方位角等于前面测量的太阳方位角α'=164.58°时，得到相应的虚拟太阳高度角，此即为所求的太阳高度角λ′，读取得数为λ′=35.96°。

利用公式（8）计算得到水立方游泳馆高度为31.2626m。结果如表3所示：

表3：水立方游泳馆数据

Claims (1)

1.一种面向单幅卫星影像的自获取三角元高度计算方法，其特征在于包括下述步骤：

步骤一、从卫星影像中测量计算参数：

（1）卫星方位角α：测量卫星影像中建筑物侧棱l在像平面上的投影与正北方向的夹角，即为卫星方位角；以正北方向为起始位置，顺时针方向旋转为正角；

（2）太阳方位角α'：测量卫星影像中建筑物侧棱阴影在像平面上的投影与正北方向的夹角，即为太阳方位角；

（3）建筑物侧棱投影长度L_pb：在像平面上，根据卫星影像中的比例尺获得建筑物侧棱l在像平面上的投影线段长度；

（4）侧棱顶点P及其阴影点的像点间距离L_ps：在像平面上，根据卫星影像中的比例尺获得建筑物侧棱l的顶点P在像平面上的投影点p与顶点P的阴影点在像平面上的投影点s之间的距离；

步骤二、计算建筑物当地当时的太阳高度角；

步骤三、三角元高度计算关系式的推导：

将公式tanλ＝L_h/L_pb代入公式 $L_h=L_{\mathrm{ps}}/\sqrt{\frac{1}{{\tan }^2{\mathrm{\λ}}^{\′}}+\frac{1}{{\tan }^2\mathrm{\λ}}-\frac{2\·\cos (\mathrm{\α}-{\mathrm{\α}}^{\′})}{\tan {\mathrm{\λ}}^{\′}\·\tan \mathrm{\λ}}},$ 有：

$L_h=L_{\mathrm{ps}}/\sqrt{\frac{1}{{\tan }^2{\mathrm{\λ}}^{\′}}+\frac{1}{{(L_h/L_{\mathrm{pb}})}^2}-\frac{2\·\cos (\mathrm{\α}-{\mathrm{\α}}^{\′})}{\tan {\mathrm{\λ}}^{\′}\·(L_h/L_{\mathrm{pb}})}}$

其中，λ为卫星高度角，是卫星拍摄方向与像平面的夹角，λ′为太阳高度角，是太阳光线方向与像平面的夹角，L_h为建筑物的高度，

计算得到建筑物高度 $L_h=\tan {\mathrm{\λ}}^{\′}\·(L_{\mathrm{pb}}\·\cos (\mathrm{\α}-{\mathrm{\α}}^{\′})+\sqrt{{L_{\mathrm{pb}}}^2\·{\cos }^2(\mathrm{\α}-{\mathrm{\α}}^{\′})+{L^2}_{\mathrm{ps}}-{L^2}_{\mathrm{pb}}}).$

Images