太阳影子定位 2015 数学建模 国赛 A题.pdf

1 太阳影子定位方法的研究与应用太阳影子定位方法的研究与应用 摘摘 要要 太阳影子定位技术是确定视频的拍摄地点和拍摄日期是视频数据分析的重要方法。本文对太阳高度角、赤纬角、方向角、时角进行了研究，建立了影子长度变化的数学模型，并分析了影子长度和模型各参数之间的变化规律，采用全局最优搜索算法和广度遍历策略对题目中的四个问题逐一求解，对所建立的四个模型采用残差分析法等进行了准确性及稳定性检验。 问题一：本文研究了影子长度与拍摄日期、时间、地点和杆高之间的关系，通过太阳高度角、方向角、赤纬角和时角公式联合建立了影子长度变化的数学模型： 同时分析了日期序号n、 北京时间t与当前地区经度、 纬度四个参数与影子长度的变化规律。应用建立的模型求解了 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）3 米高的直杆的太阳影子长度变化的二次曲线方程，参见公式(9)，画出了曲线图，参见图 11。 问题二：本文利用正午日影计算当地经度，根据问题一的数学模型，采用最小二乘法将影长和时间的关系进行二次拟合得到影长随时间变化的拟合二次曲线方程，求得真太阳时为 12 点时最短影长，继而确定直杆所在地的经度；同时在确定最短影长后使用全局最优搜索法和广度遍历在一段纬度范围内求得不同纬度相对应的杆高，再应用问题一的模型求出不同纬度所对应的影长，最后利用误差分析的原理，确定了直杆高度和当地纬度。根据附件 1 中给出的影子顶点坐标数据，求得北半球直杆高度为 1.9747m、东经 111.015 度，北纬 24.5799 度或高度为 1.7571 米、东经 111.015 度，南纬 5.145 度。 问题三： 采用问题二的方法先确定了地点， 随后将一年中的日期进行全局最优搜索，循环应用问题二中确定直杆所处地点的数学模型，最后根据误差分析得到最优解。依据建立的数学模型，根据附件 2 和 3 中给出的影子顶点坐标数据，分析了真实影长与模型计算影长之间的误差，确定了日期，详见表 11 和 12。 问题四：本文采用一种半自动的方法对视频进行阴影检测。首先对视频预处理，按4 秒间隔提取 611 张图像，然后构建一张背景图片，通过背景剪切技术检测出阴影，最后通过区域生长算法确定图像中直杆阴影的左右端点，利用图像坐标系的比例关系求出了不同时刻直杆的影子长度。由问题二中的模型确定了在已知日期的情况下的若干个可能的拍摄地点。由问题三中的模型确定在拍摄日期未知的情况下的拍摄地点与日期。 关键词关键词： 最小二乘拟合 全局最优搜索 背景剪切技术 区域生长算法 误差分析 2 一一 问题的提出问题的提出 如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中直杆的太阳影子变化，确定视频拍摄的地点和日期的一种方法。 1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场 （北纬 39度 54 分 26 秒,东经 116 度 23 分 29 秒）3 米高的直杆的太阳影子长度的变化曲线。 2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。 将你们的模型应用于附件 1 的影子顶点坐标数据， 给出若干个可能的地点。 3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件 2 和附件 3 的影子顶点坐标数据，给出若干个可能的地点与日期。 4附件 4 为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为 2 米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。 如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？ 二二 模型假设模型假设 1.由于太阳距离地球较远，假设地球接受的太阳光线为平行光； 2.物体本身不弯曲； 3.不考虑特殊的天气情况； 4.忽略日序与年的关系； 5.当地地面与当地水平面平行； 6.忽略附件 4 视频中摄像机的抖动； 赤纬角：指观测时刻的太阳倾角，即太阳直射光线和赤道平面之间的夹角； 时角：指太阳所在的时圈和通过正南方的时圈之间的夹角； 太阳高度角：地平面和太阳光的入射方向之间的夹角； 方位角：指太阳直射光线在地平面上的投影与地平面正南方之间的夹角； 当前地区纬度（度数） 当前地区经度（问题一中的单位：弧度制） st 真太阳时：真太阳视圆面中心的时角加 12 小时 三三 符号说明符号说明 3 n 一年中的日期序号，如 1 月 2 日， 2n ；10 月 22 日，295n ； H 固定直杆的高度； L 直杆被太阳光照射后的影子在地表的长度； 程序所求影长与问题二附件数据求得影长的方差； 程序所求影长与问题二附件数据求得影长的误差精度； t 北京时间； 四四 模型的建立与求解模型的建立与求解 （一）问题一（一）问题一 4.1.1 问问题分析题分析 题目要求建立影子长度变化的数学模型，分析出影子长度关于模型各参数的变化规律，并根据已给的日期、时间、地点、杆高数据，应用所建立的模型画出影子长度变化曲线。从题目要求可分析出，影子长度与拍摄日期、时间、地点和杆高是相关的，可通过研究天文知识和太阳位置相关理论，通过太阳高度角、方向角、赤纬角和时角公式联合建立影子长度变化的数学模型，在此基础上对各参数进行定量和定性分析。 4.1.2 模型的建立与求解模型的建立与求解 1.太阳位置基本概念太阳位置基本概念 在实际生活当中，地球绕着太阳运转，但人们为了描述太阳在天空的移动位置，假想存在一个以观察者为中心，以无限长为半径的球体，天空中所有星体，都在该球体上绕地轴转动，这个假想的球即为天球。太阳在天球上的位置时刻都在变化，可以用赤道坐标系和地平坐标系从不同角度来确定其位置1。 天文学中，赤道坐标系是将地球上的经纬度坐标系扩展到天球，与赤道面平行的纬度圈在天球上则为赤纬圈，通过南北极的经度圈在天球上称为时圈。太阳的位置可用赤纬角和时角表示，如图 1 所示。所谓赤纬角，是指观测时刻的太阳倾角，即太阳直射光线和赤道平面之间的夹角。而时角则是指太阳所在的时圈和通过正南方的时圈之间的夹角。从天球北极看，逆时针方向为负1。 天文学中，地平坐标系是以子午圈为主圈，地平圈为基圈，南点为主点的坐标系。可以用太阳高度角和方位角来描述太阳的位置，如图 2 所示。对于地球上的某个地点，太阳高度角是指地平面和太阳光的入射方向之间的夹角。太阳方位角即太阳所在的方位，指太阳直射光线在地平面上的投影与地平面正南方之间的夹角，以正南方向为零度，向西为正值逐渐变大，向东为负值逐渐变小，直至在正北方合为1801。 4 图 3 表示太阳在不同时刻的高度变化以及形成的高度角。图 4 中，表示高度角，表示方位角，Z表示天顶角，即光束与竖直方向的夹角。 图 1 高度角与方位角 图 2 赤纬角与时角 图 3 太阳的运动过程与形成的高度角 图 4 高度角、方位角与天顶角 2.建立影子长度变化的数学模型建立影子长度变化的数学模型 太阳在天球坐标系的位置可以用高度角和方位角来确定，计算公式如下2： 其中为当地纬度，为太阳赤纬角，为时角。 太阳赤纬角公式如下3： 5 其中2 (1)/365bn，n为一年中的日期序号，如当日期为 1 月 1 日时，1n ；当日期为 3 月 22 日，81n 。 一年中春分秋分时0，夏至时23.5，冬至时23.5； 为时角，即用角度表示时间，每15为1h，且在正午时0 ，上午0 ，下午0 ，st为真太阳时，公式如下： 其中st 为真太阳时,t为北京时间，为当前地区经度。 从图 5 可知，当太阳在运动的过程中，影子的长短会随着太阳位置的变化而变化；通过图 6 的直角三角形中， 可以直观的得到影长L、 杆高H和太阳高度角之间的关系： 根据以上公式从而建立出影子长度变化的模型： 其中： 且：2 (1)/365bn 图 5 太阳位置移动后影长的变化图 图 6 杆高、高度角和影长的关系图 (7) 6 3.分析影子长度关于各参数的变化规律分析影子长度关于各参数的变化规律 在上一节中，我们已得出影子长度L变化相关的四个参数：一年中的日期序号n、北京时间t与当前地区经度、纬度。为了分析影子长度L关于某一个参数的变化规律，我们固定其他三个参数值。 (1)影子长度L与北京时间t的变化关系 以北京天安门广场9：0015：00之间时刻为例，用Matlab运行 Test_4.m程序拟合出影子长度与时间的关系图，如下图7所示。 随着时间t的增长，影子长度L逐渐减小，在到达最低点后再增大。其中最低点为(11.9533,3.6633),即当北京时间t为 11 点 57 分时，影子长度L最短，约为3.66m。由此可知，时间决定太阳的位置，位置决定影子的长短，影子长度与最低点成左右对称关系。 图 7 影长与北京时间的关系 图 8 影长与日期序号的关系 (2)影子长度L与日期n的变化关系 以北京天安门为例，在一年当中n的取值范围在1,365，规定以当天影长的均值作为计算第n天影子长度L的取值，关系如图 8 所示。随着日期n的增大，影子长度L先逐渐减短，降至最低点后逐渐增长。基本符合抛物线。最低点为(173,1.5445)，最高点为（355,7.7369）。即在一非闰年的 6 月 22 日，影子长度为全年最低1.54m；12 月 22 日，影子长度为全年最高7.7369m。从公式(3)可知，日期的变化影响着太阳赤纬角的变化，因为赤纬角是由于地球绕太阳运行造成的现象，所以它随时间而变。赤纬角以年为周期，在+23 26与-23 26的范围内移动， 7 成为季节的标志。每年 6 月 21 日或 22 日赤纬达到最大值+23 26称为夏至，即影子最短的时候。至12月21日或22日赤纬减至最小值-23 26为冬至，此时阳光斜射北半球，此时影子最长。 (3)影子长度L与当前地区经度的变化关系 同样， 第n天影子长度L以当天影长的均值进行计算。中国领土位于北半球，东起黑龙江与乌苏里江的汇合处(东经 135 05),西到帕米尔高原(东经 73 40)，我们取考虑中国境内(东经90度东经130度)一部分区域，关系如图9所示。 影子长度L随着经度的增大而逐渐减短，至最低点后又逐渐增长。从公式(5)、(6)可知，经度与地方时间相关。不同经度的地方，日出日落有先有后。由于我国的标准时间北京时间是东经120度的地方太阳时，因此图中影长最短在东经 120 度。 图 9 影长与经度关系 图 10 影长与纬度关系 （4）影子长度L与当前地区纬度的变化关系 在北半球，纬度的范围在(0,)2 。直杆高度H、经度固定的情况下，求解出影子长度L与纬度关系如图 10 所示。随着纬度的增加，第n天的影长最大值先增大，在北纬 80.21 度时突然骤减，影子长度L随着纬度的增大，反而递减。在太阳赤纬角、时角一定时，太阳的高度随着纬度变化，纬度高，太阳高度小，纬度低，太阳高度大，因此纬度高的地方影子长，纬度低的地方影子短。 8 4.求解影子长度的变化求解影子长度的变化 如图 11，在不同时刻中，太阳位置发生变化，照射一直杆PO的顶点P也在不同的位置，即：1p、2p和p处。若研究更多的太阳位置，得到的投影顶点所形成的图形应为一条曲线。因此猜想直杆的太阳影子长度随北京时间 9：00-15：00 的变化曲线应为一条抛物线。 图 11 太阳位置移动后影长的变化图 通过 Matlab 计算，从时刻 9 点到 15 点，每增加 0.5 小时计一次当地时间，已知当地经度为东经 116 度 23 分 29 秒，取116.3度，t为北京时间，根据公式(5)得当地时间，即真太阳时st为表 1 中的数据： 表 1 真太阳时与时间的关系表 t 9 9.5 10 10.5 11 11.5 12 st 8.753 9.253 9.753 10.253 10.753 11.253 11.753 t 12.5 13 13.5 14 14.5 15 st 12.253 12.753 13.253 13.753 14.253 14.753 已知真太阳时st后，可以根据公式(4)，求得时角w的数据（转换成弧度制），如表 2： 表 2 时角与真太阳时的关系表 st 8.753 9.253 9.753 10.253 10.753 11.253 11.753 w -0.8500 -0.7191 -0.5882 -0.4573 -0.3264 -0.1955 -0.0646 st 12.253 12.753 13.253 13.753 14.253 14.753 w 0.0663 0.1972 0.3281 0.4590 0.5899 0.7208 求得时角w后，进而可以计算赤纬角的度数，已知日期为 10 月 22 日，因此295n ，根据公式(3)可求得： 最后，根据求出来的赤纬角、时角与已知的转换成弧度制的当地纬度 0.6965，并且根据公式(1): 9 求得高度角与时角的关系，如表 3： 表 3 高度角与视角的关系表 -0.8500 -0.7191 -0.5882 -0.4573 -0.3264 -0.1955 -0.0646 0.3870 0.4634 0.5317 0.5898 0.6358 0.6678 0.6842 0.0663 0.1972 0.3281 0.4590 0.5899 0.7208 0.6841 0.6675 0.6353 0.5891 0.5309 0.4625 最后，根据计算结果与公式(6),得到影子长度L与高度角的关系为： 表 4 影子长度和高度角的关系表 0.3870 0.4634 0.5317 0.5898 0.6358 0.6678 0.6842 L 7.3607 6.0025 5.1002 4.4823 4.0645 3.8036 3.6780 0.6841 0.6675 0.6353 0.5891 0.5309 0.4625 L 3.6788 3.8061 4.0689 4.4891 5.1101 6.0170 即可以得到影子长度L与北京时间t的关系为： 表 5 影子长度和北京时间的关系表 t 9 9.5 10 10.5 11 11.5 12 L 7.3607 6.0025 5.1002 4.4823 4.0645 3.8036 3.6780 t 12.5 13 13.5 14 14.5 15 L 3.6788 3.8061 4.0689 4.4891 5.1101 6.0170 因此经过数据曲线拟合，得到下图影子长度L和高度角的关系： 图 12 直杆影长岁时间的变化曲线 10 将图 12 中的曲线进行二次拟合，可以得到影长随时间变化的拟合二次曲线方程： 通过 Matlab 运行 Test_1.m 程序可以求得， 影长L最短时刻：12.2h,即，此时的最短影长为：3.6633m 4.1.3 结果分析结果分析 上述结果中，表示出了理想状态下影子长度变化的数学模型，其中影子最短的时刻为 12 点 12 分，影子最长的时刻为 9 点，较符合标准。而在整个地球上空，太阳光束射进地球的大气层，由于大气层密度一定比真空中的密度大，因此光束会发生折射现象，因此结果有一定的误差性。如图 13 与图 14，太阳近似看成光源 B，当光束照进大气层后，由光束折射定理可以得到，光束偏转，与法线的夹角偏小，折射后光束射在地表 C处。而反向延长 C 处光束，我们将点 A 近似看成了太阳所在位置，而其真正的位置在点 B 处，因此，对于太阳位置来说，得到的结果相较真正太阳位置来是偏高的。 如果不受大气层影响，太阳从 B 射入，应该直射在地表 D 处，如果发生折射现象，高度角会变大，如果当地云层较厚，大气层密度较大，所形成的高度角会继续偏大；因此经过折射后，我们计算的高度角相较没有折射情况的高度角，大小是偏大的，并且大气层密度越大，高度角越大。 因此如果有一竖立直杆，没有大气层的情况下形成的影子长度为Ll，而在大气层的影响下影子长度会逐渐变短，如图 14 的折射情况，影子长度就为l，而如果当地云层较厚，大气层密度较大，相同的直杆所形成的影子长度会继续偏小；因此经过折射后，我们计算的影子长度L相较没有折射情况的影子长度Ll，长度是偏短的，并且大气层密度越大，影子长度越短。 图 13 折射现象 图 14 折射现象引起的误差 11 （二）问题二（二）问题二 4.2.1 问题分析问题分析 题目中要求根据某固定直杆在水平地面投影出的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点，即需要确定该直杆所在地的经纬度。从问题一的求解过程中可知，经度可以通过真太阳时计算关系确定，而通过问题一中影子长度变化与时间的关系可知在真太阳时为 12 点时影子长度最短，因此只要找到最短影长所对应的北京时间便可确定直杆所在的经度；纬度的确定涉及到影长、杆高、赤纬角和高度角的关系。确定最短影长后使用全局最优搜索法和广度遍历策略在一段纬度范围内求得不同纬度相对应的杆高，再运用问题一中的模型求出不同纬度所对应的影长。最后利用误差分析的原理，根据附件 1 中给出的影子顶点坐标数据，分析真实影长与模型计算影长之间的误差，找出方差最小的便可以得到所对应的直杆高度H和纬度。 4.2.2 模型的建立与求解模型的建立与求解 1.影长变化的曲线影长变化的曲线 图 15 太阳位置变化引起的影长与影子顶点形成的曲线 在空间建立坐标系，如图 15，以地平面为XOY平面，竖直位置为Z方向，设直杆顶端一点P的坐标为( , )x y，直杆被太阳光照射投影在地面的影子顶端为点P，设其坐标为00(,)xy，则直杆的影子长度2200Lxy。 根据问题一所求出的数学模型，采用最小二乘法将影长和时间的关系曲线图进行二次拟合得到影长随时间变化的拟合二次曲线方程，推断出该曲线图也有一个最低点，因此拟合出该曲线图后，设曲线中的最低点横坐标即北京时间为mt，纵坐标即相对应的最短影长为mL。 2. 求直杆所在地的经度求直杆所在地的经度 当mtt时，根据真太阳时的定义，此时的st即为 12。根据公式(5)： P Z X Y 12 2()/312stt 可以求出mt时刻，直杆所在地区的经度： 3.求解高度角与杆高求解高度角与杆高 (1)求正午时刻的高度角 . 根据已知的日期序号n与公式(3)求出太阳赤纬角： 在太阳真时12st 即正午时， 正午时刻的高度角、 赤纬角与当地纬度的关系如下： 首先考虑北半球，南半球的分析将同理得出。所有区域的纬度范围在,02 之内。利用广度遍历策略，将此范围内的纬度每 0.001 弧度计一次，根据公式(11)与全局最优搜索的纬度1、2.1571可以表示出与纬度相对应的正午时刻的高度角(1,2.1571)ii (2)求杆高 根据在(1)中求出的正午时刻高度角与最小二乘拟合得出的最短影长mL，并结合公式(6)，可以将每个纬度所对应的杆高表示出来，即(, .)1 2 1571iH i (3)求任意高度角 根据公式(4)与公式(5)可以将时角用北京时间与已求出的经度表示出来，即： 因此在公式(1)中： 13 赤纬角、时角已求出，再次利用广度遍历策略，将,02 范围内的纬度每 0.001 弧度计一次，根据公式(3)、公式(12)与公式(1)与此 1571 个全局最优搜索的纬度1、2.1571可以表示出当时角取(1,2.21)ii纬度为(1,2.1571)jj时，相对应的高度角为： 根据公式(3)、 公式(12)与公式(1)与此 1571 个全局最优搜索的纬度1、2.1571可以表示出与纬度相对应的高度角(1,2.21,1,2.1571)ijij 图 16 纬度变化范围 4.利用误差分析求解最优纬度利用误差分析求解最优纬度 根据杆高(1,2.1571)iH i 与 (3)求出的任意高度角，并且运用公式(6)，求出影子长度(1,2.1571)iL i ； 根据附件 2 中所给的坐标数据计算出影长2200Lxy，与iL求两者方差，即： 2121()ijiijiLL 14 取方差最优值min，并且使得min ， 为误差精度。 当方差取最优值min时，根据对应的影长iL求解出所对应的杆高和高度角，并根据公式(1)求解出最优纬度0，即为此模型中所求直杆所在地区的纬度。 5.根据附根据附件一数据求解模型件一数据求解模型 设直杆被太阳光照射投影在地面的影子顶端为点P其坐标为00(,)xy， 数据为附件 1中的 21 个数值。则直杆的影子长度2200Lxy 表 6 附件一数据对应的影子长度 1.1496 1.1822 1.2153 1.2491 1.2832 1.3180 1.3534 1.3894 1.4262 1.4634 1.5015 1.5402 1.5799 1.6201 1.6613 1.7033 1.7462 1.7901 1.8350 1.8809 1.9279 将影长和时间的关系曲线图进行二次拟合，得到影长随时间变化的拟合二次曲线图： 图 17 影长随时间变化的拟合二次曲线图 从拟合出的二次曲线图，得出曲线中的最低点横坐标即北京时间为mt， 12.599mt 纵坐标即相对应的最短影长为mL， 0.4936mL 当mtt时，根据真太阳时的定义，此时的st即为 12。根据公式(5)可以求出mt时刻，直杆所在地区的经度： 1 1 1 . 0 1 5 15 已知附件 1 中数据日期为 2015 年 4 月 18 日，因此日期序号108n ，根据赤纬角公式： 求出赤纬角： 0.1841 .在12st 即正午时，由地理知识可以得到正午时刻的高度角、赤纬角与当地纬度的关系： 首先考虑北半球，所有区域的纬度范围在,02 之内。利用广度遍历策略，将此范围内的纬度每 0.001 弧度计一次，根据公式(11)与此 1571 个全局最优搜索的纬度1、2.1571可以表示出与纬度相对应的正午时刻的高度角： (1,2.1571)ii 表 7 全局最优搜素出的纬度与相对应的正午高度角关系图 1 2 3 4 5 6 . 1571 0.001 0.002 0.003 0.004 0.005 0.006 . 1.571 1 2 3 4 5 6 . 1571 1.1671 1.1661 1.1651 1.1641 1.1631 1.1621 . -0.4029 由于数据较多，只列举一部分。 分析求出的正午时刻高度角(1,2.1571)ii与最小二乘拟合得出的最短影长mL，并且根据公式(6)可以将每个纬度所对应的杆高表示出来，即： (1,2.1571)jHj 根据公式(4)与公式(5)以将时角用北京时间与已求出的经度表示出来，即： 5123t 16 表 8 附件一中的北京时间与时角的关系表 1t 2t 3t 4t 5t 6t . 21t 14.70 14.75 14.80 14.85 14.90 14.95 15.00 15.05 1 2 3 4 5 6 . 21 0.5500 0.5631 0.5762 0.5893 0.6024 0.6155 0.6286 0.6417 在公式(1)中：赤纬角已求出、时角已表示出，再次利用广度遍历策略，将,02 范围内的纬度每 0.001 弧度计一次，根据公式(3)、公式(12)与公式(1)与此 1571 个全局最优搜索的纬度1、2.1571可以表示出当时角取(1,2.21)ii纬度为(1,2.1571)jj时，相对应的高度角为： 考虑求出的杆高(1,2.1571)iH i 与高度角(1,2.21,1,2.1571)ijij， 并且仍然运用公式(6)，求出影子长度(1,2.21,1,2.1571)ijLij 比较所得结果iL与题目中所给的坐标计算出的影长2200Lxy， 求出两者方差， 即：2121()ijiijiLL 取方差最优值min，使得min，为误差经度 根据 Matlab 运行 Test_2.m 程序得到： 61.2511 10 当方差取最优值时， 根据对应的影长ijL求解出所对应的杆高H=1.9747m和高度角，并根据公式(1)求解出最优纬度，即为此模型中所求杆高的纬度： 24.5799 因此求解出若直杆在北半球，则直杆所处地点和杆高为： 纬度纬度 经度经度 杆高杆高H 北纬北纬 24.5799 东经东经 111.015111.015 杆高杆高 1.9747m1.9747m。 表 9 北半球直杆所在地点与杆高 17 同理对应南半球，在同一经度时，运用上述模型，求得： （三）问题三（三）问题三 4.3.1 问题分析问题分析 问题三要求根据某固定直杆在水平地面上的太阳影子顶点坐标数据建立数学模型，以确定直杆所处的地点和日期。地点在问题二中已经解决，只需确定日期。嵌套利用广度遍历策略，将一年中的日期进行搜索，循环运用问题二中确定直杆所处地点的数学模型，最后根据误差计算可得到最优解。依据建立的数学模型，根据附件 2 和 3 中给出的影子顶点坐标数据，分析真实影长与模型计算影长之间的误差，可得到所对应的日期n和纬度。 4.3.2 模型的建立与求解模型的建立与求解 1. 根据日期根据日期求求解直杆的所在地解直杆的所在地 在问题二的基础上，以日期n作为一个变量（1,365n），搜索在1,365范围内的日期， 每一个日期序号计一次， 每一个日期根据问题二建立的模型确定直杆的经纬度。 根据公式可以确定日期序号为(1,2.365)in i 时所对应的 365 个直杆所在地区的经度与纬度。 确定日期序号为(1,2.365)in i 时所对应的 365 个最小方差值min(1,2.365)ii与日期序号的关系图像，找到最小方差值的最小值。 因此可以确定直杆所在地区的经度为： (20)12imt， 纬度为： (1,2.365,1,2.21,1,2.1571)ijkijk 2.2.根据数据求解模型根据数据求解模型 将附件二的坐标与时间数据分别带入上述模型中，得到以下数据与关系图。 纬度纬度 经度经度 杆高杆高H 南纬南纬 5.145 东经东经 111.015111.015 杆高杆高 1.7571m1.7571m。 表 10 南半球直杆所在地点与杆高 18 图 18 附件二中的数据得出纬度与日期关系图 图 19 附件二中的数据得出最小方差值与日期关系图 从图中可以观察出，最小方差值中有两点较低，取到极小值点。通过 Matlab 运行Test_3.m 程序得出这两点的日期与纬度为： 表 11 附件二中数据得出的直杆所在经纬度 n（日期） 102（4月月 12日）日） 246（9月月 3日日） （东经） 71.8650 71.8650 （北纬） 32.7519 32.3418 从表格中可以观察出， 在一定误差范围内， 两点的经纬度较为相似。 根据地理常识，我们知道太阳在一年中的直射点在北回归线-赤道-南回归线-赤道-北回归线范围内运动，考虑北半球上的区域，在一年之中的确有两天使得该区域上的直杆影子的长度变化 19 规律相同，据此判断，求出的两组经纬度表示的是同一区域。 同样将附件三的坐标与时间数据分别带入上述模型中，得到以下数据与关系图： 图 20 附件三中的数据得出纬度与日期关系图 图 21 附件三中的数据得出最小方差值与日期关系图 从图中可以观察出，最小方差值中有两点较低，取到极小值点。通过 Matlab 得出这两点的日期与纬度为： 表 12 附件三中数据得出的直杆所在经纬度 n（日期） 97（4月月 7日）日） 248（9月月 5日）日） （东经） 108.9600 108.9600 （北纬） 43.6021 43.9459 同理南半球的直杆所在地点与日期也同样用问题三中的模型进行求解，不再赘述。 （四四）问题）问题四四 4.4.1 问题分析问题分析 由于题目没有给出摄像机的内参和外参，因此要想根据附件四中一根直杆在太阳下 20 的影子变化视频与已知的直杆长度建立确定视频拍摄地点的数学模型，首先需要从视频中提取图像，然后构建一张背景图片，通过背景剪切技术检测出阴影，通过区域生长算法确定图像中直杆阴影的左右端点和坐标，最后利用图像坐标系的比例关系可求出不同时刻直杆的影子长度。根据坐标，由问题二中的模型可以确定在已知日期的情况下，确定若干个可能的拍摄地点。由问题三中的模型可以确定在拍摄日期未知的情况下，确定出拍摄地点与日期。 4.4.2 模型的建立与求解模型的建立与求解 1.求解实际影子长度求解实际影子长度 用Matlab读入视频，分析可得帧数为61020帧，帧速率为：25帧/秒，宽度为1920像素，高度为1080像素。首先我们从视频每隔4秒提取一张图片，共提取得到611张图片，然后构建一张背景图片，通过背景剪切技术检测出阴影，通过区域生长算法确定了图像中直杆阴影的左右端点，利用图像坐标系的比例关系求出不同时刻直杆的影子长度，如图19-21。从处理后的像素坐标中得到固定直杆的底端坐标00(,)xy，直杆顶端坐标( ,)x y以及直杆投影在地面的影子顶端坐标( ,)iix y。其中影子顶端坐标根据递归函数，向右衍生，能到达的最远点，即为影子在上述坐标中的顶点坐标。 得到坐标之后，我们就可以求得杆高长度，即为22000()()Hxxyy，与影子长度2200()()iiiLxxyy，其中1,2.61i 。 根据公式iiHLHL，得到： iiHLLH 其中L为实际的影子长度，H为实际固定直杆的高度(2m)。 图22 视频中原始图像 21 图23 背景图像 图24 背景剪切后的图像 2.确定拍摄地点与日期确定拍摄地点与日期 (1)根据附件四中的视频，在已知拍摄日期的情况下，通过问题二中的模型得出若干个可能的拍摄地点： 表13 求解出的拍摄地点的经纬度 纬度（N） 22.92 51.57 45.84 51.57 经度（E） 113.70 170.30 128.5 142.5 (2)根据附件四中的视频，在拍摄日期未知的情况下，通过问题三中的模型求解出拍摄地点与日期。 在建模过程中，我们确实得到了若干个可能的拍摄地点及对应的日期，经过分析，所得结果回代检验后，符合题中要求。在此并未给出具体的数据！ 22 五五 模型的检验模型的检验 问题三中的模型是关于日期n、纬度、经度、北京时间t、影子长度L和直杆高度H的模型。现在作以下模型检验，验证模型的合理性。 问题一中，根据北京时间t、经纬度、和杆高H求出了影子长度L。现在把影子长度L当作已知条件，由模型三反求经纬度，进行模型一的检验。结果如下表所示: 表15 检验出的实际值与计算值的相对误差表 实际值 计算值 相对误差 纬度（N） 39.9072 41.5012 0.0399 经度（E） 116.3914 114.2554 0.0184 问题二中，根据顶点坐标数据（即影子长度L）、日期n、北京时间t求出了经纬度和直杆高度。反过来，现在把经纬度、直杆高度H、日期n和北京时间t当作已知条件，由模型三反求影子长度L，进行模型二的检验。作出残差图23如下所示。散点图上的点均在0附近，可以知道计算影子长度与给定影子长度高度吻合。 图 23 实际值与计算值的残差图 由于问题三中的模型是根据问题一、二中的模型得出的，当上述两种检验方式效果较好， 说明问题三所建立的模型是合理的。 此时， 可以说明三个模型都是稳定和可靠的。而问题四的求解模型由模型二和模型三推出，那么模型四也是可靠的。据此，证明了本论文中提出的四个模型具有一致性、合理性。 23 六六 模型模型的的评价评价和推广和推广 （一）（一）模型优点模型优点 1.本文根据地球纬度、经度、太阳高度角、日期、时间、赤纬角、时角、固定直杠的高度等信息求解出了具体时间和地点的影子长度，利用多种信息更接近实际情况； 2.文章中把太阳影子顶点坐标都统一成影子长度，计算更为方便； 3.本文的一大特色是对四个模型均进行了模型的检验，得出了经纬度的相对误差和影长的残差图。事实证明所建立的模型效果很好。 （二）（二）模型模型不足不足 1.模型忽略了特殊天气的影响，只考虑正常天气； 2.模型四中不考虑相机拍摄角度对影字长度的影响，这样会存在一定的误差。 （三三）模型推广模型推广 1.第一问中推导处理影子长度的计算公式，使得影长在后面的表示求解有了根据。 2.题目中多问要求确定一些参数求解另外一个或多个参数，我们建立了多个求解模型。事实上可以建立一个整体模型，输入若干参数，得到若干输出。这样模型的应用就更广泛了。 3.考虑到模型通过提取视频建立模型，求解拍摄地点和日期，我们可以建立一个能够处理图像和视频的数学模型。 七七 参考文献参考文献 1张闯，吕东辉，太阳实时位置计算及在图像光照方向中的应用，200072，TN911.7； 2单黎明，太阳跟踪定位技术及其应用研究，北京控制工程研究所，北京 100190； 3lishucheng_71，太阳赤纬，/view/862819.html，2015 年 9 月 13 日； 4朱志红，陈为，毛行奎，自动太阳跟踪器综述C，中国电源学会第 18 届全国电源技术大会，厦门，2009； 5谷林，胡泽东，陈良益，等.用 CPLD 替代单片机实现线阵 CCD 自动变扫描控制J，光电子激光，2003，14（8），830-834 6吴翊、吴孟达、成礼智,数学建模的理论与实践，长沙：国防科技大学出版社，1999； 7彭肖约，何平安，袁炳夏，基于 CPLD 的线阵 CCD 驱动电路设计与实现J，光电子激光，2007，18（7），803-807 24 附录一附录一：各程序代码：各程序代码 （一）（一）、问题一的求解程序问题一的求解程序 求解影子长度模型的代码求解影子长度模型的代码 %问题一的求解程序 Test_1.m clear,clc; n = datenum(2015,10,22)-datenum(2015,0,0);%日期转换成天数 b=2*pi*(n-1)/365; delta =deg2rad(0.006918-0.399912*cos(b)+0.070257*sin(b)-0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+0.00148*sin(3*b)*(180/pi);%赤纬角 fai = deg2rad(39+54/60+26/3600);%纬度 sc = (120-116.3)/15;%时差 st = 9-sc:0.1:15-sc; omega = deg2rad(st-12)*15);%时角 alpha = asin(sin(delta)*sin(fai)+cos(delta)*cos(fai)*cos(omega);%高度角 H = 3;%杆长 L = H*cot(alpha);%影长 %曲线拟合 s = 9:0.1:15; plot(s,L); title(天安门广场直杆影子长度变化曲线); xlabel(时间轴); ylabel(影长); grid on; hold on; %多项式拟合 A=polyfit(s,L,2); %n 是给定的多项式的次数，拟合出来的结果 A 是系数向量 L1=polyval(A,s); %计算出拟合的 L 值 plot(s,L1,r-); %红线是拟合曲线 legend(所求影长曲线,二次拟合曲线); (二二) 问题二的求解程序问题二的求解程序 %问题二的求解程序 Test_2 clear,clc; %将Excel里的时间数据用函数 =HOUR(A4)+MINUTE(A4)/60 转换成数值形式,然后存储到txt文件中 data = load(data_2.txt);%加载数据 time = data(:,1); time = time; X = data(:,2); X = X; Y = data(:,3); 25 Y = Y; L = sqrt(power(X,2)+power(Y,2); plot(time,L,*); hold on %最小二乘法拟合 x = 8:0.001:15.7; p = polyfit(time,L,2); y = polyval(p,x); plot(x,y); hold off ml,mt = min(y); minTime = 8+mt*0.001; E = 120-(minTime-12)*15 n = datenum(2015,4,18)-datenum(2015,0,0); %日期转换成天数 b=2*pi*(n-1)/365; delta =deg2rad(0.006918-0.399912*cos(b)+0.070257*sin(b)-0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+0.00148*sin(3*b)*(180/pi);%赤纬角 fai = 0:0.001:pi/2; %北半球 %fai = -pi/2:0.1:0; %南半球 sc = (120-E)/15; %时差 st = time-sc; omega = deg2rad(st-12)*15); %时角 for i=1:length(fai) alpha(i) = pi/2-fai(i)+delta;%北半球正午高度角 %alpha(i) = pi/2+fai(i)-delta;%南半球正午高度角 H(i) = ml*tan(alpha(i);%杆长 alpha1 = asin(sin(delta)*sin(fai(i)+cos(delta)*cos(fai(i)*cos(omega);%高度角 l = H(i)*cot(alpha1);%影长 error(i) = var(l-L); end minErr = 100; index = 0; for k=1:length(error) if all(minErr error(k) minErr = error(k); index = k; end end 26 fai(index)*180/pi H(index) (三三) 问题问题三三的求解程序的求解程序 %问题三的求解程序 Test_3 clear,clc; %将Excel里的时间数据用函数 =HOUR(A4)+MINUTE(A4)/60 转换成数值形式,然后存储到txt文件中 data = load(data_3_1.txt);%加载数据 time = data(:,1); time = time; X = data(:,2); X = X; Y = data(:,3); Y = Y; L = sqrt(power(X,2)+power(Y,2); plot(time,L,*); hold on %最小二乘法拟合 x = 8:0.001:18; p = polyfit(time,L,2); y = polyval(p,x); plot(x,y); hold off ml,mt = min(y); minTime = 8+mt*0.001; E = 120-(minTime-12)*15 N = zeros(365,1); for n=1:365 b=2*pi*(n-1)/365; delta =deg2rad(0.006918-0.399912*cos(b)+0.070257*sin(b)-0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+0.00148*sin(3*b)*(180/pi);%赤纬角 fai = 0:0.001:pi/2; %北半球 %fai = -pi/2:0.1:0; %南半球 sc = (120-E)/15; %时差 st = time-sc; omega = deg2rad(st-12)*15); %时角 for i=1:length(fai) alpha(i) = pi/2-fai(i)+delta;%北半球正午高度角 %alpha(i) = pi/2+fai(i)-delta;%南半球正午高度角 H(i) = ml*tan(alpha(i);%杆长 alpha1 = asin(sin(delta)*sin(fai(i)+cos(delta)*cos(fai(i)*cos(omega);%高度角 27 l = H(i)*cot(alpha1);%影长 error(i) = var(l-L); end minErr = 100; index = 0; for k=1:length(error) if all(minErr error(k) minErr = error(k); index = k; end end N(n) = fai(index)*180/pi; high(n) = H(index); minErrs(n) = minErr; end minErr1 = 100; minErr2 = 100; index1 = 0; index2 = 0; for k=1:length( minErrs) if (minErr1 minErrs(k) minErr1 = minErrs(k); index1 = k; end end for k=1:length( minErrs) if (minErr2 minErrs(k)&(minErrs(k)=minErr1) minErr2 = minErrs(k); index2 = k; end end index1 FAI1 = N(index1) index2 FaI2 = N(index2) t = 1:365; figure(1); plot(t,N,r.); title(纬度随日期变化图); xlabel(日期轴); ylabel(纬度); legend(所求纬度坐标); figure(2); plot(t,minErrs,*); 28 title(所求方差随日期变化图); xlabel(日期轴); ylabel(方差值); (四四) 问问题题四四的求解程序的求解程序 %提取视频中的数据 %提取图片 mov = VideoReader(Appendix4.avi); for i=1:100:mov.numberofframes b=read(mov,i); imwrite(b,strcat(m,int2str(i),.jpg),bmp); end % 获取背景图片 Files=dir(m*.jpg); z_num=length(Files); Sum=zeros(1080,1920); for i =1:z_num It=imread(Files(i).name); Xt=rgb2gray(It); %取灰度图 Sum=Sum+double(Xt); end Ver=Sum/z_num; %求平均 jpg=uint8(Ver); imshow(jpg); imwrite(jpg,ver.jpg,jpg) %得到所有图片以后求平均 Files=dir(j_m*.jpg); z_num=length(F