1、全国大学生数学建模竞赛国家一等奖论文 A题 太阳影子定位模型的分析太阳影子定位模型摘要针对太阳影子定位问题，本文结合地理学和天文学的相关知识，建立了不同 数据类型下的太阳影子定位模型，实现了视频拍摄地点和日期的快速精准确定。 对于问题一，首先从地理学角度，基于地理坐标，直杆长度，时间这三个影 响影子长度的参数，计算出时角，赤纬角，太阳高度角，进而给出了影子长度与 三个参数之间的关系式。结果显示，影长对日期和时刻都呈现出先减小后增大的 趋势；对杆长呈正比关系增长；对经度呈现先急剧增长到峰值再突变为 0 ，而后突变到峰值后再急剧下降；对纬度呈缓慢上升趋势。然后，根据附件 1 中提供的 数据，画出了
2、天安门广场上直杆的太阳影子分布曲线图。 对于问题二，基于问题一中对影响影子长度因素的分析，根据地理学知识建20 21 i i 归i 归i立双目标规划模型，确立目标函数分别为：min | A - A | , min | S- S |。i=1i=1然后在约束条件下对杆子的地点坐标应用网格逼近算法优化求解，得出最符合题 目所提供数据的杆子地理位置为：(19.1E,108.71N ) 海南东方市境内，此时， 杆长为 2.03米 ，太阳方向角残差比为1.8% ，影长残差比为 0.9%，误差均很小。 对于问题三，首先建立了与问题二相似的目标规划模型，由于日期未知，模 型求解的时间复杂度较高。为提高计算速度
3、，引入了粒子群算法。分别对附件 2 和 3 中 的 数 据 进 行 分 析 ， 确 定 出 的 地 点 坐 标 分 别 为 (80.51E，32.13N ) ， (110.20E，24.83N ) 和 (81.43E，32.24N ) ，(111.56E，23.68N ) ，附件 2 为西藏阿里，日期为 8 /14 或 4 / 29 ，附件 3 为广西梧州市，日期为12 / 27或12 /14 。可以 发现，两种算法的结果极为接近，但粒子群算法计算时间要远小于网格逼近算法。 对于问题四，首先对视频数据进行采集和预处理，由于视频拍摄角度的存在， 从视频中直接得到的影长并不是实际长度，而是其投影长
4、度，这里采用基于 Hough变换和透视变换的图像矫正法，对斜视图像进行矫正，得出实际影长。然 后将得到的数据带入问题二的模型中，给出视频拍摄地点为 (110.70E，42.31N )内蒙包头市境内；在拍摄日期未知的情况下，将变化而来的实际影长代入问 题 三 的 基 于 粒 子 群 算 法 的 目 标 规 划 模 型 ， 求 解 出 视 频 拍 摄 地 点 为 (109.76E，42.66N ) 内蒙包头市境内，拍摄日期为 6 /11或 7 /13 。 对于模型的推广，根据物体采集到的太阳地理信息进行计算，可以应用到求 建筑物群合理间距问题。 关键词：双目标规划 粒子群算法 Hough 变换 透
5、视变换 一、 问题的重述1.1 问题的背景“日长影移”是生活中人人熟知的自然现象，这个词说明地面上的影子变化 与太阳活动有着密切的联系。而古代智慧的先民就利用了这个现象制作了日晷， 是最早且最精确的计时工具之一。在图像信息充斥的当代，如何通过图像数据获 得图像拍摄时的相关信息是图像分析学科的重要课题，而利用太阳光影变换获得 时间和地理信息，是非常方便可靠的。 1.2 问题的提出太阳影子定位技术，是通过分析视频中物体的太阳影子变化，确定视频拍摄 的地点和日期的一种方法。针对上述背景和应用需求，提出以下问题： 1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律， 并应用你们建立的模
6、型画出 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）3 米高的直杆的太阳影 子长度的变化曲线。 2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确 定直杆所处的地点。将你们的模型应用于附件 1 的影子顶点坐标数据，给出若干 个可能的地点。 3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型 确定直杆所处的地点和日期。将你们的模型分别应用于附件 2 和附件 3 的影子顶 点坐标数据，给出若干个可能的地点与日期。 4如果已有一根直杆在太阳下的影子变化的视
7、频，并且已通过某种方式估 计出直杆的高度为 2 米。请建立确定视频拍摄地点的数学模型，并应用你们的模 型给出若干个可能的拍摄地点。若拍摄日期未知，能否根据视频确定出拍摄地点 与日期。 二、 问题的分析2.1 问题一的分析题目要求在固定地点，给定日期和杆长的条件下，求解出直杆投影长度的变 化曲线。对于水平地面上的垂直直杆，直杆长度与影子的比值即为太阳高度角的 正切值，因此需要知道此时间段内的太阳高度角变化。查阅资料6 可得，太阳高 度角与当地地方时、经纬度密切相关，因此知道上述两个量就可确定直杆的变化 过程。 2.2 问题二的分析题目要求根据影子的变化情况和给出的日期求出直杆的位置，实际上是问题
8、一的逆求解过程。这里杆长和地点都是未知量，逆求解是非常困难的，于是将问 题二转化为双目标规划问题。当太阳方位角与影长的实际值与理论值差值的绝对 值之和达到最小时，所得经纬度即为杆子的地点坐标。 2.3 问题三的分析与问题二不同的是，该问中日期是个未知量。首先考虑沿用上一问的模型。 由于日期未知，所以要考虑一年 365 天的所有情况，这将大大增加运算时间。从 减少运算量的角度考虑，有必要改进算法。考虑引入现代优化算法之一的粒子群 算法，将所有解视为粒子所要去的位置，由于适应度与目标函数相联系，选取合 理的适应度函数，期望提高计算效率。 2.4 问题四的分析本问的数据由视频给出，那么首先要对视频进
9、行数据预处理。由于视频时间 较长，所以不考虑将视频逐帧分析，而是每隔一段时间对视频进行获取分析，对 图像所给出的影长和影子角度进行测量。考虑到视频拍摄角度的存在，从视频中 直接得到的影长并不是实际长度，而是其投影长度，因而要对斜视图像进行矫正， 得出实际影长。然后再分别利用问题二和问题三的模型进行求解。 三、 问题的假设1. 每年的太阳活动情况是相同的，均为“恒星年”。 2. 地球是一个完美的球形，不考虑海拔、地球扁率的影响。 3. 无光线衍射造成的影子减淡现象。 5. 在小尺度考虑直杆投影问题时，地表是绝对水平的。 6. 不考虑地球公转的影响。 7. 题目所给的数据是真实的，可靠的。 四、
10、符号说明符号说明 太阳高度角 L 水平地面上直杆长度 S 水平地面上直杆影子长度 h 当地地方时时角 当地太阳赤纬 当地纬度值 当地经度值 t当地地方时 t0北京时间 n当日日期序号 P当地经纬度坐标 A当地太阳高度角 T时间（包含时刻与日期） 离散化后的方位朝向 i每组数据影长与 x 坐标轴的夹角 xi ， yi题目附件提供的第 i 组直杆坐标 d日期 注：其他符号将在下文中给出具体说明。五、 模型的建立与求解5.1 问题一5.1.1 模型的建立本模型结合相关地理学知识，对影子的变化情况进行分析描述。下面将明确 一些地理学定义，以及重新定义一些本模型需要用到的参数。 太阳高度角，也称太阳高度
11、，是指某地的太阳光线与当地地平面的所交的最 小线面角，这是以太阳视盘面的几何中心和理想地平线所夹的角度。在水平地面 上，直杆长度与影长的比值即为太阳高度角的的正切值： tan = L (1) S图 1 太阳方位角示意图 其中 L 为杆长， S 为影长。 通过查阅 相关资料6 ，得知太阳高度角的计算公式为： sin = cos h cos cos + sin sin (2) 其中 为太阳高度角， h 为地方时时角， 为当时的太阳赤纬， 为当地纬度。 以一个地方太阳升到最高的地方的时间为正午 12 时，将连续两个正午 12 时之间等分为 24 个小时，所成的时间系统，称为地方时。地球上每一个地点都
12、有 其相应的地方时。由于题目只提供了当地时间的北京时间，因此在计算地方时时 角时，要先将北京时间换算为当地地方时 t ： -120t = t0 +15 (3) 其中 t0 为北京时间， 为当地经度。 根据某地地方时，可以换算出当地的地方时时角。地方时时角 h 即为当地与子午线之间相差的角度： h = 15(t -12) (4) 太阳的赤纬等于太阳入射光与地球赤道之间的角度，由于地球自转轴与公转 平面之间的角度基本不变，因此太阳的赤纬随季节不同而有周期性变化。太阳赤 纬的最高度数为 2326 ，夏至时太阳的赤纬为 +2326 ，冬至时太阳的赤纬为-2326。春分和秋分时太阳的赤纬为 0。 由于地
13、球公转轨道的偏心率非常低，可以看作是一个圆圈，太阳赤纬 可用 下面这个公式来计算： = 23.45sin 2 (284 + n) (5) 365 其中 n 为当日日期序号，1 月 1 日时 n =1，以此类推得 10 月 22 日 n = 295 。 联立式子 (2) - (5) 得到方程组： Ltan =Ssin = cos h cos cos + sin sin t = t0 + -12015 (6) h = 15(t -12) + = 23.45sin 2 (284 n) 求解上述方程组，得： 365 S = L / tan (arcsin (cos(15(t -12)cos( )cos
14、 + sin ( )sin ) (7) 可见，影子长度的变化与当地地理位置 P (, )，直杆长度 L ，时间 T (t, d ) 这三个参数有关。 5.1.2 模型的求解首先计算题目所给条件下的 h ， 与 ，再将上述参数值代入 (2) 式，得出从 9:00-15:00 的太阳高度角随时间的变化（具体值见附录）。相应地，可求得直杆 影长数据（具体见附录）。从结果中挑选出几个比较重要的时间点，将相应结果 制成下表以供参考： 表 1 影子长度分布 北京时间 太阳高度角/度 影长/m 9:00 21.18 7.74 12:00 37.88 3.85 12:14 37.99 3.84 15:00 2
15、5.35 6.33 将影长随时间变化的情况用 MATLAB 绘制成图像： 图 2 影子长度变化曲线 对于影子长度关于各个参数的变化规律，从时间，地点和杆长三个方面对其进行分析，分析某一个参数时，将其他变量看做常量，只改变一个未知量。我们 将时间参数分为日期和时刻这两种情况进行运算求解，影长关于时刻的变化在图 2 中已经给出；对于地点，将其分为经度 和纬度 这两个参数进行运算求解， 求解后的变化规律图如下： 图 3 影长与日期关系曲线 图 4 影长与杆长关系曲线 图 5 影长与经度关系曲线 图 6 影长与纬度关系曲线 5.1.3 结果分析从图 3 中可以发现，影长随着时间的增加，呈现先减小后增大
16、的趋势，影长 最小点出现在12 :14 ，这是由于北京时刻为120E 的地方时，换算到116E 附近时， 会产生时差，显然是符合常理的。 从图 4 中可以发现，影长与杆长呈正比关系，这是由于 tan = L ，从而验证S了模型的准确性。 从图 5 中可以发现，影长随时间变化出现了两个明显的“脉冲”。这是因为 当太阳逼近地平线时，影长变化的速度非常快，影子也是最长的。而太阳一旦在 地平线以下，就不存在影子，影子长度也为 0。 从图 6 中可以发现，图 5 中的曲线随纬度增加，总趋势是增加的，这是随着 纬度的增加太阳高度角减小，从而导致了影长的增加。 5.2 问题二基于网格逼近算法的双目标规划模型
17、5.2.1 数据预处理分析题目所给的数据（附件 1），可以发现这些数据不仅能计算出影长，还能 够得知影子方位角的变化，但是附件没有给出坐标系的方位朝向，故假设坐标系 x 轴的方位朝向角为 (0 360) ，取正北方向时，x 轴的方位朝向角为 0。 由此，可以得出 21 个时间点影子的太阳方位角。我们以每五个时间点为间隔，选 取部分呈现在下表中： 表 3 影子方位角数据表 北京时间 x 坐标/米 y 坐标/米 方位角/度 14:42 1.0356 0.4973 + 25.63 14:57 1.2087 0.5255 + 23.49 15:12 1.3955 0.5541 + 21.65 15:2
18、7 1.6033 0.5833 + 20.02 15:42 1.8277 0.6135 +18.55 5.2.2 模型的建立首先引入太阳方位角的定义。太阳方位角是太阳在方位上的角度，它通常被 定义为从北方沿着地平线顺时针量度的角： cos A = sin cos - cos h cos sin cos(8) 上述公式可以用来计算近似的太阳方位角，不过因为公式是使用余弦函数， 所以方位角永远是正值，因此，角度永远被解释为小于 180 度，而必须依据时角来修正。当时角为负值时 (上午)，方位角的角度小于 180 度，时角为正值时 (下午)，方位角应该大于 180 度，即要取补角的值，故作如下修正：
19、 arccos sin cos - cos h cos sin ,h 0 (21) 式中， Ai 表示实际的太阳方位角， S归i 表示归一化后的实际影长。 5.3.2 模型的求解同样，对经纬度进行离散化处理后，选取步长为 0.01，对其进行 MATLAB 编 程求解： 表 7 杆子所处地点和日期及误差分析数据表 附件二 附件三 经度 80.51E 经度 110.20E 纬度 32.13N 纬度 24.83N 杆高 2.04m 杆高 3.10m 太阳方向角残差比 2.0% 太阳方向角残差比 2.3% 影长残差比 1.1% 影长残差比 0.35% 日期 8 /14 或 4 / 29 日期 12 /
20、 27或12 /14 5.3.3 模型的改进基于粒子群算法的目标规划模型在实际求解中，发现由于加入新的未知量日期，使得模型的可行域大大增加， 增加了遍历算法的时间复杂度，因此，我们引入 粒子群算法9 对目标函数进行求 解，从而降低模型求解的时间复杂度。 假定有一个 D 维的目标搜索空间，有 n 个微粒组成了一个粒子群，其中每个 微 粒 都 用 一 个 D 维 的 向 量 描 述 ， 将 它 的 空 间 位 置 表 示 为 mi = (mi1, mi 2 , , miD ),i = 1, 2, , n ；这可看做目标优化问题中的一个解，代入适应 度函数计算出适应度值可以衡量微粒的优劣；第 i 个微粒的飞行速度也是一个 D 维的向量，记为 vi = (vi1, vi 2 , , viD )；第 i 个微粒所经历过的具有最好适应值的位 置称为个体历史最好位置，记为 pi = ( pi1, pi 2 , , piD )；整个
