CN105116470A

CN105116470A - 基于地物tin模型的日照时数计算

Info

Publication number: CN105116470A
Application number: CN201510375777.0A
Authority: CN
Inventors: 张富; 毛大鹏; 张丽娟; 邱本志
Original assignee: North China University of Water Resources and Electric Power
Current assignee: North China University of Water Resources and Electric Power
Priority date: 2015-06-29
Filing date: 2015-06-29
Publication date: 2015-12-02

Abstract

本发明提供一种基于地物TIN模型的日照时数计算，包括以下步骤：建立地物TIN模型的数据结构，数据结构中单个点的内存结构包括：大地坐标系，站心坐标系，空间直角坐标系，它们之间的转换标志Flag；计算全天日照时间，根据设定的步长对全天日常时间进行离散，形成离散的时刻序列；计算时刻序列中某一时刻点的太阳高度角和方位角，然后循环地物TIN列表中的地物数据，计算出模型中某一点在受影高度平面上形成的投影，判断需要计算日照时数的点是否被遮挡，一旦发生遮挡立即退出循环；最终累计出需要计算日照时数的点被遮挡时间，计算出日照时间段。本发明减少了经纬坐标导致的误差，改进了日照时数的计算方法，显著地提高了日照时数的计算精度。

Description

基于地物TIN模型的日照时数计算

技术领域

本发明涉及复杂地形的日常时数计算，尤其涉及一种基于地物TIN模型的日照时数计算。

背景技术

地球形体以及遮挡与被遮挡地物所处的空间位置对日照时数具有重要的影响，但是在传统建筑日照分析研究中，把整个分析区域作为同一经纬度下的理想平面，不仅忽略了地球的形体特征，同时也忽略了遮挡与被遮挡地物所处经纬坐标的差异。根据《椭球大地测量学》理论，当地球表面距离大于10公里时，法截线的计算必须考虑地球的形体特征。以南京(东经118°48′，北纬32°04′)为例，2008年12月21日(冬至日)上午7：30，按照直角三角形边角关系，产生10公里长的阴影只需700米高的地物。

反之，地物之间的经纬坐标差异将呈现在高度角和方位角的计算上，从而影响地物的阴影范围计算，最终导致日照时数计算出现较大误差。

发明内容

本发明要解决的技术问题在于提供一种日照时数计算方式，以解决现有技术存在的问题。

为解决上述技术问题，本发明采用以下技术方案：

基于地物TIN模型的日照时数计算，包括以下步骤：

建立地物TIN模型的数据结构[tin₁，tin₂，……，tin_m]，所述数据结构中单个点的内存结构包括：由经度L、纬度B、高程H构成的大地坐标系，由站心经度L′、站心纬度B′、和高程H′构成的站心坐标系，以地心为中心的X、Y、Y空间直角坐标系X，Y和Z，还包括它们之间的转换标志Flag；

计算出日出日落时间，得出全天日照时间，根据设定的步长对全天日常时间进行离散，形成离散的时刻序列[T₁，T₂，……，T_n]；

计算时刻序列中某一时刻点T_i的太阳高度角和方位角，然后循环地物TIN列表中的地物数据，计算出tin_j在受影高度平面上形成的投影S_j，然后判断需要计算日照时数的点是否被S_j遮挡，一旦发生遮挡立即退出循环；

然后继续循环T_i+1…T_n,计算出得出需要计算日照时数的点被连续遮挡的时间段序列{[T_x...T_x+n],[T_x',...,T_x'+m],...[T_x”,...,T_x”_+q]}，并累计出需要计算日照时数的点被遮挡时间T_N，然后计算出日照时间段T_M＝T_T-T_N。

判断需要计算日照时数的点是否被S_j遮挡时，如果需要计算日照时数的点在S_j内部，所述内部包含边界，则视为tin_j对需要计算日照时数的点产生遮挡。

所述需要计算日照时数的点被遮挡的判断方法为：

当受影面高度大于0米时，认为遮挡需要计算日照时数的点的地物在受影面上的投影为“平面”TIN，此时，如果平面上需要计算日照时数的点在平面的外部，那么地物TIN不对需要计算日照时数的点不产生遮挡，如果平面上需要计算日照时数的点在平面的内部，那么地物TIN不对需要计算日照时数的点产生遮挡；

当受影面高度等于0米时，认为地物在地球的椭球面上投影，此时，需要计算日照时数的点位于地球的椭球面上，认为椭球面上的TIN为球面，通过面积比较法判断需要计算日照时数的点与球面三角形的拓扑关系，如果平面上需要计算日照时数的点在球形三角形的外部，那么球面三角形对需要计算日照时数的点不产生遮挡，如果平面上需要计算日照时数的点在平面的内部，那么球面三角形对需要计算日照时数的点产生遮挡。

如果受影面与地物TIN相交，则地物TIN在受影面上投影产生四边形或者三角形或者直线或弧线。对于直线或弧线，则需要判断需要计算日照时数的点是否在线上来判断是否产生遮挡；而对于四边形，则需要分割为两个三角形，然后对每个三角形判断需要计算日照时数的点在三角形的内部还是外部，来判断是否产生遮挡。

本发明的有益效果：本发明减少了经纬坐标导致的误差，改进了日照时数的计算方法，显著地提高了日照时数的计算精度。

附图说明

图1为地物TIN模型的内存数据结构。

图2为地面点的地心空间直角坐标与日照二维平面坐标对应关系。

图3为地物独立性遮挡原理示意图。

图4为基于地物TIN投影的日照时数计算流程。

图5为点与球面三角形的关系。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细说明。

因为复杂的地理场景中包括建筑、地形、植被和标志物等地物要素，所以选择TIN数据模型表达地物数据不仅能有效实现多源数据的融合和转换，而且简化了地物投影与测试点拓扑关系的判断过程。

本发明使用天文算法即《AstronomicalAlgorithms》(Meeus,1998)公布的方法(TEP-SPA的代表算法)作为计算太阳位置的SPA算法。

如图1所示，本发明中，建筑、地形和其他地物类型统一采用TIN结构表达，单个TIN则由无序的三点(_3DPoint自定义结构)组成。为避免浮点计算过程中忽略极小量，提高计算精度，需要对地心空间直角坐标统一转换为站心坐标系进行相关计算。因此，单个点的内存结构需要包含：大地坐标系(经度L、纬度B和高程H)、站心坐标系(站心经度L′、站心纬度B′、和高程H′)、地心空间直角坐标系(X，Y和Z)以及相互之间的转换标志Flag。

因此，上述中可以实现地物原始数据类型到TIN数据结构的转换。由于建筑、地形和其他地物类型的表面空间数据都可以简化为点、面和体要素类型，因此最终都可转换为TIN结构的三角面片单体数据类型。

在天文经纬与大地经纬中存在一定差异，垂线偏差是导致其差异的主要因素之一，而坐标系的也会导致垂线偏差的出现。为达到最好的效果，本发明采用WGS84大地(质心)坐标系作为本研究的坐标系参照，下述中经纬度即指WGS84坐标系下经纬度。为区分天文坐标系和大地坐标系，下述将采用“地球坐标系统”表示适用于日照分析研究的最佳大地(质心)坐标系统，本文即指WGS84大地质心坐标系。

本发明首先构建一个二维平面坐标系，即由遮挡地物的遮挡点、地面点和太阳位置构建的平面二维坐标系统(以下简称“日照坐标系”)。

然后，为了描述与计算方便，需对同一个点的不同坐标系表达之间进行转换。本发明中采用WGS84大地(质心)坐标系作为坐标系参照，下述中经纬度即指WGS84坐标系下经纬度。同时为区分天文坐标系和大地坐标系，下述将采用“地球坐标系统”表示适用于日照分析研究的最佳大地(质心)坐标系统，本发明中即指WGS84大地质心坐标系。同时设定由遮挡地物的遮挡点、地面点和太阳位置构建的平面二维坐标系统(以下简称“日照坐标系”)。

如图2所示，首先构建地面点P在地心空间直角坐标系下(X,Y,Z)与日照坐标系下(x,y,z)的对应关系，实现步骤如下：

(1)构建空间坐标系与P点站心坐标系(x，y，z方向不变)关系式

把大地空间直角坐标系O-XYZ平移至P点，变换成以P点为中心的坐标系P-XYZ，由于P在O-XYZ坐标系下坐标为：N为P点卯酋圈曲率半径，e为椭球第一偏心率，L、B为P点经纬度，H_p为P点大地高(如果P点在球面上则H_p＝0)。

此时，空间任一点(不与P重合，假设为M)在P-XYZ坐标系中坐标(X_p-M,Y_p-M,Z_p-M)与O-XYZ坐标系下坐标(X_M,Y_M,Z_M)的对应关系为：

(2)P点站心坐标系转换为P站心地平直角坐标系(x向地平北，z向天顶)

针对坐标系P-XYZ，首先绕Z轴顺时针旋转经度180°-L(L为P点经度)，然后Y轴反向变换为Y′，再绕Y′轴顺时针旋转90°-B(B为P点纬度)，则形成以P点为中心，经纬切线方向分别为X′、Y′的新直角坐标系P-X′Y′Z′。则点M在P-X′Y′Z′坐标系中坐标变为：

其中，

因此可得：

(3)P站心地平直角坐标系转换为日照二维坐标系(x为日照投影方向)

假设某一时刻，太阳方位角为A，如果坐标系P-X′Y′Z′绕Z′轴顺时针旋转A(方位角上午为负，下午为正)，形成新坐标系P-xyz,在P-xyz坐标系中，日照光线与遮挡物在一平面内。而此时点M在P-xyz坐标系中坐标变为：

由公式(1)和公式(2)可得到P-X′Y′Z′和O-XYZ关系。又由于R为旋转矩阵(也为正交矩阵)，因此存在：R^-1(θ)＝R^T(θ)和

由此可以得出P-xyz与O-XYZ关系为：

由于日照光线与测试点和地面点组成的直线在同一平面内，因此，相对某一时刻而言y≡0，即：

(4)确定日照坐标系与大地经纬坐标系关系

地心空间直角坐标与经纬坐标之间的关系式：

采用公式(7)可构建出O-XYZ、P-xyz和(L，B)之间的对应关系，即：由O-XYZ变换为P-xyz采用公式(5)，由P-xyz变换为O-XYZ采用公式(6)，根据经纬度和方位角计算O-XYZ采用公式(7)。上述公式中，L、B为P₀和P₁纬度和经度，N为(L，B)位置对应的卯酋圈曲率半径，e为椭球第一偏心率；H_p为计算点大地高(如果P处于地面上，H_p＝0)。

接下来，进行日照投影公式的求取。

在空间直角坐标系下，M₀为遮挡体，M₂和M₃为被遮挡物，M₀大地高为H₀，M₃大地高为H₁。P₀为M₀地面点，P₁为M₀某一点(图示为最高点)；P₂为被遮挡体M₃地面点，P₃为M₃顶点。假设P₃为测试点，遮挡P₃的地物为M₀，那么M₀完全不遮挡P₃点日照光线的充分必要条件为：以M₀上任一点P_i和太阳位置所构建出的日照光线，在与距地面高为H₁的阴影面的交点，不等于P₂点经纬度。于是日照遮挡问题的求解转化为：计算某一时刻P_i点的日照光线，与一定大地高的阴影面交点问题。

P₀为遮挡地物M₀地面点，P₁为遮挡地物上产生阴影的某一点，P₁距P₀距离为H₀。在某一时刻，P₁处太阳高度角为h，G2为过P₁点的日照光线，P₃为测试点，也即G2与受影面交点。P₂与P₃具有相同大地坐标但高度为0的地面点。假设P₀处经纬度分别为L和B，P₂处经纬度即为待求经纬度Lx和Bx，则P₁和P₃大地经纬坐标分别为(L，B，H₀)和(L_x，B_x，H₁)。

在P-xyz坐标系下构建P₀、P₁、P₂与P₃的几何关系，G1为日照立面与地球椭球面相切而形成的法截弧。如果让P-xyz坐标系的Z轴与P₀P₁重合，指向天顶，X轴则定义为P₀点在椭球面上的切线，指向太阳射线G2在地面投影方向；P₄为G2与x轴交点，某一时刻P₁处太阳高度角为h是已知量，现在只需要根据B、L、H₀、H₁以及某一时刻的时间变量即可计算出P₂和P₃经纬坐标L_x、B_x，再通过判断L_x、B_x与实际测试点位置关系即可判断出P₀P₁是否遮挡P₃。

求取L_x、B_x的过程如下：

①构建P-xyz坐标系下的G1法截弧曲线方程

G₁(x,z)＝K₁x²+K₂z²+K₃z+K₄xz

其中：a为椭球长半轴，b为短半轴；A为P₁点某一时刻对应的太阳方位角。

②构建求解P₂点在P-xyz坐标系下坐标z₂、x₂方程组

假设P₂、P₃在P-xyz坐标系下坐标分别为(z₂，x₂)和(z₃，x₃)，同时根据几何关系可以确定P₀(0，0)，P₁(H₀，0)，P₃(0，H₀coth)，其中h为某一时刻对应的太阳高度角。

条件一：P₃满足直线方程G2，其中G2为：z+xtanh-H₀＝0，即：

f₁＝z₃+x₃tanh-H₀＝0(8)

条件二：P₂满足曲线方程G1，即：

条件三：O-XYZ坐标系下，P₂、P₃在XOY平面投影点P₂′、P₃′和原点O共线。

由于P₂和P₃具有相同纬度和经度，因此O-XYZ坐标系下，P₂、P₃在XOY平面投影点P₂′、P₃′和原点O共线。根据公式(6)得出：

f₃＝(C₁x₃+C₂z₃+C₃)(D₁x₂+D₂z₂+D₃)-(D₁x₃+D₂z₃+D₃)(C₁x₂+C₂z₂+C₃)＝0(10)

其中：

条件四：P₂至P₃直线距离为H₁，即：

由以上四条件得出的f₁、f₃和f₄可消去z₃和x₃，得到z₂，x₂二元方程组：

其中：

如果采用牛顿逼近法求解x₂和z₂，收敛8-12次即可得到10^-12数据精度。

③求解P₁投影点P₃的大地经纬(L_x，B_x)

把P₂点P-xyz坐标系下坐标值(x₂，0，z₂)带入(6)得出P₂点在O-XYZ坐标系下X₂、Y₂和Z₂值，再根据公式(7)反向计算(H_p＝0时)出L_x、B_x，即为P₂和P₃点经纬度(L_x，B_x)。

由日照分析原理可知，日照时数的计算方法主要是通过判断地物的投影面与测试点的空间拓扑关系，从而判断出某一时刻地物是否遮挡测试点，最终累计出整个分析时段的日照时间和被遮挡时间。根据图形几何原理可知：地物阴影无论是否受到其他地物的“阻挡”或者变形，都不影响此地物对测试点的遮挡性能，即“地物独立性遮挡”原理。

如图3所示，如果地物A能对O点产生遮挡，即使地物B完全“阻挡”A的阴影导致A的阴影无法到达O点，那么B肯定具有遮挡O点“能力”，否则并不能完全“阻挡”住能遮挡O点的A阴影。反之，如果地物A不能对O点产生遮挡，则无论B是否阻挡A的阴影都不影响A对O的遮挡关系。同理，由于D无法完全“阻挡”C的阴影，因此C对Oˊ的遮挡状态与D无关。

基于上述原理，把地物TIN逐个投影在受影面上，针对某一TIN的遮挡判断时，不需考虑周边其他遮挡物的影响，仅计算此TIN在某时刻的投影与测试点的空间拓扑关系，从而循环判断出某一时刻是否存在地物TIN对测试点产生遮挡。最终累计出测试点O被遮挡的总时间数，再计算出O点的日照时数。

其实现流程图4所示：

首先，根据设定的地理、日期和大气参数，根据太阳位置算法计算日出日落时间，从而得出全天日照时间作为日照分析时间段T_T(也可根据需要设定)，并由设定时间步长对时间进行离散，形成后续计算中的时刻序列[T₁,T₂,...T_n]。

循环时间序列[T₁,T₂,...T_n]，把某一时间段作为时刻点T_i并计算出对应太阳高度角和方位角。循环TIN列表[tin₁,tin₂,...tin_m]中地物数据，计算出tin_j在受影高度平面上形成的投影S_j，通过测试点O与受影面上投影S的拓扑关系判断O是否被tin_j遮挡。即：如果测试点O在S内部(含边界)，则视为tin_j对O产生遮挡。一旦发生遮挡即时退出循环。

继续循环T_i+1…T_n,最后得出O点被连续遮挡的时间段序列{[T_x...T_x+n],[T_x',...,T_x'+m],...[T_x”,...,T_x”_+q]}，并累计出O点的被遮挡时间T_N，则日照时间段T_M＝T_T-T_N。

在计算过程中，需要判断测试点是否被遮挡。由于地球为椭球形体，因此，根据测试点高度构建出的受影面决定了地物TIN的投影面形状。地物TIN在H高度的受影面上投影形状可分为以下几种：

(1)受影面高度H＞0米

当受影面高度H＞0米时，可以认为遮挡测试点的地物在受影面上的投影仍为“平面”TIN。此时，只需判断出平面上点(测试点)是否在平面三角形(地物TIN投影)的“内外”，既可判断出地物TIN是否对测试点产生遮挡。如果平面上需要计算日照时数的点(测试点)在平面的外部，那么地物TIN不对需要计算日照时数的点不产生遮挡，如果平面上需要计算日照时数的点在平面的内部，那么地物TIN不对需要计算日照时数的点产生遮挡。

(2)受影面高度H＝0米

当投影面高度H＝0米时，可以认为地物是在地球椭球面上投影。“椭球面”上TIN可以认为“球面”TIN。因此，问题转化为球面三角形与球面上某一点的空间拓扑关系判断。通过面积比较法判断需要计算日照时数的点与球面三角形的拓扑关系，如果平面上需要计算日照时数的点在球形三角形的外部，那么球面三角形对需要计算日照时数的点不产生遮挡，如果平面上需要计算日照时数的点在平面的内部，那么球面三角形对需要计算日照时数的点产生遮挡。

点在球面多边形的位置关系判断相对比较复杂，不仅与多边形的类型有关，而且与投影中心的位置有关。本发明中比较面积法进行判断点与球面三角形的拓扑关系。

如图5所示，A、B、C分别为某地物TIN在地球球面上三个投影点，P为测试点。采用比较面积法判断P点是否在球面三角形ABC中：连接大弧PA、PB、PC，分别计算出球面三角形ABC、APB、APC、BPC的面积S、S₁、S₂和S₃，如果满足S＝S₁+S₂+S₃则P在球面三角形ABC内部，否则P在球面三角形ABC外部。

此外，在极端情况，如果假想受影面与地物TIN相交，则地物TIN在受影面上投影可能产生四边形、三角形、直线或弧线。对于直线或弧线，则判断测试点是否在线上；而对于四边形，则需要分割为两个三角形，然后再根据上述方法进行判断。

Claims

1.基于地物TIN模型的日照时数计算，其特征在于，包括以下步骤：

2.根据权利要求1所述的基于地物TIN模型的日照时数计算，其特征在于：判断需要计算日照时数的点是否被S_j遮挡时，如果需要计算日照时数的点在S_j内部，所述内部包含边界，则tin_j对需要计算日照时数的点产生遮挡。

3.根据权利要求1所述的基于地物TIN模型的日照时数计算，其特征在于：所述需要计算日照时数的点被遮挡的判断方法为：

4.根据权利要求3所述的基于地物TIN模型的日照时数计算，其特征在于：如果受影面与地物TIN相交，则地物TIN在受影面上投影产生四边形或者三角形或者直线或弧线。对于直线或弧线，则需要判断需要计算日照时数的点是否在线上来判断是否产生遮挡；而对于四边形，则需要分割为两个三角形，然后对每个三角形判断需要计算日照时数的点在三角形的内部还是外部，来判断是否产生遮挡。