数学建模太阳影子定位的模型建立.docx
《数学建模太阳影子定位的模型建立.docx》由会员分享，可在线阅读，更多相关《数学建模太阳影子定位的模型建立.docx（31页珍藏版）》请在冰豆网上搜索。
数学建模太阳影子定位的模型建立
太阳影子定位的模型建立
摘要
太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。
为了准确地确定视频拍摄的地点和日期，从而为视屏数据的分析作铺垫，本文基于太阳影子定位技术，运用天体物理学理论和数学几何知识，建立了影子长度变化的数学模型，及其参数的变化规律，并应用此模型画出了2015年10月22日北京时间9:
00-15:
00之间天安门广场3米高的直杆的太阳影子长度的变化曲线是二次函数；利用遗传算法、三次样条插值和图像调整坐标系，建立了任一固定直杆位置的数学模型，并得到了直杆所处的地点；投影角度与太阳角速度之间的关系，以及参数估计算法，建立了任一固定直杆位置的数学模型，并得到了直杆所处的地点与日期；利用Matlabr2014b，确定了视频拍摄地点的数学模型和可能的拍摄地点。
问题一，基于天体物理学理论和数学几何知识，利用太阳高度角公式、时角公式和赤纬角公式，计算得到了影子长度变化的数学模型及其参数的变化规律，通过Matlab软件模拟仿真得到了2015年10月22日北京时间9:
00-15:
00之间天安门广场3米高的直杆的太阳影子长度的变化曲线图。
从图中可以看到，太阳影长随时间的变化先缩短后伸长，呈现出一个较为对称的图像，并可清晰看出在中午12点左右，图像有一个最低点，此时影长大致为3.8米左右。
问题二，基于遗传算法，利用三次样条插值和图像调整坐标系，建立了直杆所处地点轨迹的数学模型。
再利用二次多项式拟合影子轨迹求出当地正午时间，即可根据北京时间求出大致经度。
最后建立太阳影子定位的优化模型，从而得出若干可能的地点。
问题三，基于某固定直杆在水平地面上的太阳影子顶点坐标数据，利用投影角度和太阳角速度之间的关系，参数估计算法，建立模型直杆所处的观测日期和经纬度，从而求出直杆所处地点的参数方程。
运用Matlab软件编程数据拟合了附件1的影子顶点的坐标数据，得到了直杆影子的轨迹图和其长度的变化图。
问题四，利用Matlabr2014b工具箱中的视频查看器播放视频获得截图，使用图片查看器中的测距工具测量视频中的直杆长度为680，从而得到视频中尺度和实际尺度的比例尺为340：
1。
从北京时间8:
54:
20每间隔3051帧（2分钟）截取图片，测量出对应时刻的视频影长l，利用Excel记录数据并计算出实际影长L。
建立模型时采用最小二乘法拟合出方程，最终得出视频拍摄点的若干可能位置。
关键词：
太阳影子定位遗传算法三次样条插值法图像调整坐标系最小二乘法Matlab软件
一、问题重述
如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。
1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:
00-15:
00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。
2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。
将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。
3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。
将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。
4．附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。
请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。
如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？
二、符号与变量说明
符号
说明
L
表示直杆影子长度
d
表示直杆高度
h
表示太阳高度角
x
表示直杆所在纬度
O
表示赤纬角
t
表示时角
S
表示北京时间
Y0
表示太阳直射点时的经度
A
表示太阳方位角
X
表示所求的纬度
Y
表示所求的经度
δ
表示当时的太阳赤纬角
φ
表示当地纬度
三、模型假设
对于问题一：
1.以直杆为主要参考对象，视太阳为东边升起，西边落下；
2.忽略地球绕太阳的公转；
3.地球为一个体积均匀的球体；
4.假设10月22日当天9时至15时，太阳直射点的纬度是一个固定值。
对于问题二：
假定太阳光束为平行光。
对于问题三：
假定地球不动，以地球为中心，以任意长为半径作一假想球面。
对于问题四：
视频拍摄视角对影子长短的影响忽略不计。
四、模型的建立与求解
4.1问题一的求解
4.1.1影子长度变化数学模型的建立
根据天体物理学理论和数学几何关系可知，任一物体的影子长度与物体所处的纬度、经度、太阳直射点纬度以及太阳的位置等有关。
通过分析物体影子长度变化随太阳位置变化来建立影子长度变化模型，并对其和各参数进行分析。
将一直杆垂直立于某水平地面上，其杆的太阳影子随太阳位置的变化而变化，得到影子长度变化的模型演示图，如图1所示。
太阳光线通过直杆顶点Q，在地面上形成一个影点Q’，设影子的长度QQ’为L，太阳光线与地面的夹角为∠QQ’O=h。
图1影子长度变化模型演示图
由图1可得
（1.1）
.（1.2）
联立公式（1.1）（1.2），并运用三角函数得影子长度的变化规律：
.（1.3）
其中，L为直杆影子长度，d为直杆高度，h为太阳高度角，X为直杆所在纬度，O为赤纬角。
4.1.2影子长度变化的数学模型参数
（1）赤纬角O
赤纬角指的是太阳赤纬，是地球赤道平面与太阳和地球中心的连线之间的夹角。
赤纬角以年为周期，在+23°26′与-23°26′的范围内移动，成为季节的标志。
每年6月21日或22日赤纬达到最大值+23°26′称为夏至，该日中午太阳位于地球北回归线正上空，是北半球日照时间最长、南半球日照时间最短的一天。
随后赤纬角逐渐减少至9月21日或22日等于零时全球的昼夜时间均相等为秋分。
至12月21日或22日赤纬减至最小值-23°26′为冬至，此时阳光斜射北半球，昼短夜长而南半球则相反。
当赤纬角又回到零度时为春分即3月21日或22日，如此周而复始形成四季。
因赤纬值日变化很小，一年内任何一天的赤纬角O可用下式计算：
（1.4）
或
.（1.5）
其中N表示自1月1日起到观测当天为第N天，t为时角，s为北京时间，则可以得到：
.（1.6）
由上式可以得到当地时间正午12：
00时的时角为0，以小时来计量（1HA=150）。
（2）太阳高度角
对于地球上的某个地点，太阳高度角是指太阳光的入射方向和地平面之间的夹角，专业上讲太阳高度角是指某地太阳光线与通过该地与地心相连的地表切线的夹角。
太阳高度角简称高度角。
当太阳高度角为90°时，此时太阳辐射强度最大；当太阳斜射地面时，太阳辐射强度就小。
.（1.7）
（3）时角
时角指的是某天体的赤经（RA）与当地的恒星时（LST）的差值，即
.（1.8）
一个天体的时角表示该天体是否通过了当地的子午圈（中天）。
其数值表示了该天体与当地子午圈的角距离，并借用时间的单位，以小时来计量（1HA=150）。
当地时间12点时的时角为零，前后每隔一小时，增加360/24=150。
从正午算起,按顺时针方向为正，逆时针方向为负，也就是上午为负，下午为正。
它的数值等于离正午的时间（小时）乘以15，故日出、日落时的时角最大，正午时为零（指地平线以上）[1]。
（4）年积日
年积日是仅在一年中使用的连续计算日期的方法，是从当年1月1日起开始计算的天数。
例如：
每年的1月1日为第1日，2月1日为第32日，以此类推。
平年的12月31日为第365日，闰年的12月31日为第366日。
4.1.3模型分析处理
假设北京时间就是天安门当地时间。
同时，由于六个小时的时间，太阳直射点的纬度变化不大，因此可以假设10月22日当天9时至15时，太阳直射点的纬度是一个固定值。
由于10月22日基本上位于区间第三部分的左起1/3处，因此由函数关系可以得到：
（
），（1.9）
即该时段太阳直射点的纬度。
同时，在北京时间中午12时之前n小时这一时刻，太阳直射点的经度比天安门所在经度靠东15n0，即
.（1.10）
在北京时间中午12时之后n小时这一时刻，太阳直射点的经度比天安门所在经度靠西15n0，即
.（1.11）
知道了直射点，就可以知道直射点与天安门的距离，继而知道直射点与地球球心的连线和天安门与地球球心的连线的夹角来了。
4.1.4直杆的太阳影长的变化曲线
利用上述建立的影子长度变化的数学模型，结合MATLAB软件编程，得到了2015年10月22日北京时间9:
00-15:
00之间天安门广场（北纬39054’26’’,东经116023’29’’）3米高直杆的太阳影子长度的变化曲线。
根据题设条件可知，观测时间为2015年10月22日，为一年中的第295天，即N=295，观测点经度为116023’29’’，且观察点得当地时间就是北京时间，故北京时间12：
00时对应时角为0，9:
00-15:
00对应的时角范围为-450-450。
利用太阳高度角公式、时角公式、赤纬角公式以及三角函数关系，通过Matlab软件运行，得到了直杆的太阳影长的变化曲线，如图2所示。
图2直杆的太阳影长的变化曲线
由图2可知，太阳影长随时间的变化先缩短后增大，呈现出一个较为对称的图像，并可清晰看出在中午12点左右，图像有一个最低点，此时影长大致为3.8米左右。
4.2问题二的求解
4.2.1问题分析
问题二要求根据顶点数据确定影子顶点坐标，打算基于三次样条插值做出影子轨迹，根据图像调整坐标系，再利用二次多项式拟合影子轨迹求出当地正午时间，即可根据北京时间求出大致经度。
考虑到题中仅给出影长的坐标，相对于问题一缺少杆长条件，于是本文打算利用拟合求出的经度作为初始经度，并对地图网格化，以杆影和实际影长的相似度作为适应度，基于遗传算法对网格中的经纬度进行搜索，建立太阳影子定位的优化模型，从而得出若干可能的地点。
4.2.2模型建立
本文通过对影子轨迹拟合估计出当地正午时间，再根据北京时间求出初始经度，以此作为初值解，建立基于遗传算法的太阳影子定位模型。
4.2.3拟合确定初始经度
对于问题二，利用二次曲线拟合来估计2015年4月18日的当地正午时间，依据附件1，对杆影长度进行拟合去逼近此动态过程，最终达到估计的目的。
将杆影附件中数据加以处理，杆影长度随时间变化，如表5.1所示。
表1附件1中影长随时间变化
时间（2015年4月18日）
14:
42
14:
45
14:
48
14:
51
14:
54
14:
57
15:
00
杆影长度
1.1496
1.1822
1.2153
1.2491
1.2832
1.3180
1.3534
再利用相应时间与表5.1中数据按二次曲线进行拟合，使用最小二乘法确定拟合方程为
（2.1）
其中，z0为附件1杆影长度，t0为附件1当地时间。
因为杆影最短的时刻即为当地正午时刻，故
，所以初始经度为
.（2.2）
其中，ϕ0为附件1初始经度，t0min为附件1杆影最短时刻。
4.2.4基于启发式算法的太阳影子定位模型
问题二要求根据直杆在水平面上的影子顶点数据推断可能的地点，建立直杆在水平地面上形成的太阳影子坐标示意图，如图3所示。
图3直杆在水平地面上形成的太阳影子坐标示意图
本文通过确定约束条件，以求解的竿影长和实际竿影长的相似度作为目标函数，建立基于启发式算法的太阳影子定位模型，具体步骤如下：
（1）初始化
根据拟合方法确定的经度为111.020，考虑到误差影响，取经度范围为[105,115]，纬度范围取[0,90]，选定步长为0.10，接下来将经纬度范围网格化，记下该网格为矩阵D，通过经纬度确定可能的地点。
（2）确定目标函数
知道太阳高度角公式
，不妨设杆长为s米，可以计算出21个时刻所对应的影长L，为消去杆长的影响，按照如下公式进行处理：
.（2.3）
可以得出影长之比M，记附件中的影长为
，可以得出实际影长之比M'：
.（2.4）
n为可能的地点，优化目标记为
.（2.5）
（3）首先以步长为0.10对经度进行循环，对于每一确定的经度
，对纬度
进行循环，利用得到的经纬度算出目标函数值
，记
，再继续搜索下一对经纬度，进而计算出目标函数值
，与
进行比较：
（2.6）
其中D为网格化的矩阵，d为优化的经纬度坐标。
（4）遍历所有经纬度，确定最优的目标值，从而确定可能的拍摄地点。
4.2.5基于遗传算法的太阳影子定位模型
由于以上的模型收敛速度慢，不适用数据量大的搜索，于是本文利用遗传算法的全局搜索能力对网格中的经纬度筛选，以影长和实际影长的比例作为适应度，建立基于遗传算法的太阳影子定位模型，模型建立流程图如图4所示。
图4模型建立流程图
4.2.6模型的应用
根据直杆在水平面上的影子顶点数据推断可能的地点，本文通过确定约束条件，以求解的竿影长和实际竿影长的相似度作为目标函数，建立基于启发式算法的太阳影子定位模型，具体步骤如下：
（1）编码
采用二进制编码方法，对于n个地点的搜索问题，染色体分为n段，其中0表示该地点未被选择，1表示选择该地点。
考虑到可能有若干个可能的拍摄地点，故将染色体中出现1的概率记为5/n，表示可能有5个拍摄地点。
（2）种群初始化
在完成染色体编码以后，必须产生一个初始种群作为初始解，本文确定初始化种群的数目为50个。
（3）适应度函数
将可能的经纬度记为矩阵
其中
表示i点的经纬度（K具体值见附件2），二进制的染色体记为
，记下M中为1的位置，对应经纬度即为可能的拍摄地点。
再利用求解的经纬度按照问题一的公式（3）得出太阳高度角
，不妨设杆长为s米，可以计算出21个时刻所对应的影长L，为消去杆长的影响，按照如下公式进行处理：
.（2.7）
可以得出影长之比M，记附件中的影长为
，可以得出实际影长之比M’：
.（2.8）
n表示可能的拍摄地点，则该染色体的适应度记为：
.（2.9）
即适应度为影长与原始影长相似度，优化的目标就是选择适应度函数尽可能大的染色体，适应度函数值越大的染色体越优质，反之越劣质。
（4）交叉变异操作
本文的交叉和变异概率选用动态参数，参数的变化与进化代数相关。
在种群早期，高交叉率，低变异率；在种群后期，低交叉率，高变异率。
用这种方法减少种群早熟的可能性。
经过多次试验发现参数为：
（2.10）
（2.11）
其中，Pc表示交叉概率，Pm表示变异概率，
表示当世代数，
表示总代数，Pc在0.875-0.958之间变化，Pm在0.0208-0.0625之间变化。
时，算法可以很好的收敛。
（5）选择
选择操作即从旧群体中以一定概率选择个体到新群体中，个体被选中的概率跟适应度值有关，个体被选中的概率跟适应度值有关，个体适应度值越大，被选中的概率越大，从而实现优胜劣汰的过程。
（6）解码
对二进制染色体解码，即选取适应度最高的染色体，其编码为1所对应的经纬度即为可能的拍摄地点。
表2附件2大致经度及地点
经度
12:
41
12:
44
12:
47
12:
50
12:
53
12:
56
12:
59
地点（海南省）
1.247
1.223
1.199
1.175
1.152
1.130
1.108
4.3问题三的求解
4.3.1天文背景知识
虽然地球绕太阳运行，但相对来说，在地球上观察到的却是太阳在天空中移动。
为了确切地描述太阳在天空中的移动与位置，假定地球不动，以地球为中心，以任意长为半径作一假想球面，天空中包括太阳在内的一切星体，均在这个球面上绕地轴转动，这个假想的球体称为天球。
延长地轴线与天球相交的两点称为天极，N为北天极，S为南天极，NS即为天轴。
扩展地球赤道面与天球相交所成的圆PP′，称为天球赤道，如图6所示。
由于黄道面（地球运行的轨道面，也就是太阳视运行的轨道面）与天轴的夹角为66°33′，则黄道面与天球赤道面夹角为23°27′，太阳沿着天球黄道周而复始地绕地球运行[2]。
图6黄道在天球上的位置图
太阳在天球上的位置每日、每时都有变化。
为了确定其位置，常用赤道坐标系和地平坐标系从不同角度来表示。
赤道坐标系是把地球上的经、纬度坐标系扩展至天球，在地球上与赤道面平行的纬度圈，在天球上则叫赤纬圈；在地球上通过南北极的经度圈，在天球上叫做时圈。
以赤纬δ和时角Ψ表示太阳的位置，如图8所示。
所谓时角，是指太阳所在的时圈与通过南点的时圈构成的夹角，单位为度。
自天球北极看，顺时针方向为正，逆时针方向为负。
时角表示太阳的方位，因为天球在一天24h内旋转360，所以每小时旋转150。
图7赤道坐标系—赤纬和时角
4.3.2模型三的建立
基于某固定直杆在水平地面上的太阳影子顶点坐标数据，利用投影角度和太阳角速度之间的关系，参数估计算法，建立模型直杆所处的观测日期和经纬度，从而求出直杆所处的地点的参数方程。
（1）太阳高度角（h）计算公式近似为
.（3.1）
其中h为太阳高度角，t为时角，δ为太阳的赤纬角，φ为当地纬度。
（2）当时角为正值时（下午）太阳的方位角（A）为，
（3.2）
当时角为负值时（上午）太阳的方位角（A）为，
（3.3）
其中A为太阳的方位角，h为太阳高度角，t为时角，δ为当时的太阳赤纬角，φ为当地纬度。
4.3.3方位角求解
根据“天球理论”，可导出太阳方位角与观测地点经度、纬度以及观测时刻的关系如下：
（3.4）
式中A为太阳方位角，h为太阳高度角，δ为太阳的赤纬，t为太阳时角。
利用三角函数关系式，有
（3.5）
地球的自传以24h为周期。
在一个周期内，太阳运行方位角距离为2π。
因此可以只考察k=0的情况。
于是式（3.4）和式（3.5）简化为：
（3.6）
（3.7）
因此，式
（1）存在3个解，即：
（3.8）
（3.9）
（3.10）
可见，解决问题的关键在于确定选择A1，A2，A3作为A的取值区间。
可以想到，对于任意时刻的太阳位角，可以取3个值中的任意一个作为结果，所以不同的是所确定的0度点存在不同。
与太阳的实际运行趋势相一致，取值必须满足两个原则，即对应于以24h为周期的时间标尺，太阳运行的方位角必须保持向正向增加或负向减小的变化趋势，为叙述简便，将此原则称为条件
（一）。
同时，太阳方位角连续变化，不存在跳跃，称为条件
（二）。
尽管式
（1）中方位角的推导以正南方向为0度点，但是，考虑到太阳东升西落的视觉感受，在跟踪系统的实际应用中，以正东方向作为0度点更加符合日常的经验。
但是，由于作者求解太阳方位角的主要目的是满足太阳跟踪设备研制的需要，考虑到日出时刻太阳位置存在着处于东北方向的可能，为了避免系统在初始位置左右跨越，带来控制上的不便，在实际设计中选取正北方向为方位角的0度点。
由此，可以确定，在真太阳时0时，A1是A的解。
通过综上分析，可以得到下面的太阳方位角求解公式，公式所解出的0度方位角为正北方向[3]。
（3.11）
上式中：
（3.12）
式中，t为真太阳时。
4.3.4太阳高度角及其计算方法
图8太阳直线轨迹变化图
太阳高度角简称太阳高度,系指太阳相对于地平面的角距离。
赤纬，一般用6表示，一年中在±23°27′的范围内变动，在从春分至秋分期间为正,在从秋分至春分期间为负[4]。
（1）正午太阳高度角的变化
对于地球上的某个地点，太阳高度是指太阳光的入射方向和地平面之间的夹角。
太阳高度是决定地球表面获得太阳热能数量的最重要的因素。
我们用h来表示这个角度，它在数值上等于太阳在天球地平坐标系中的地平高度。
太阳高度角随着地方时和太阳的赤纬的变化而变化。
太阳赤纬以δ表示，观测地地理纬度用φ表示，地方时（时角）以t表示，有太阳高度角的计算公式：
（3.13）
日升日落，同一地点一天内太阳高度角是不断变化的。
日出日落时角度都为零度，正午时太阳高度角最大。
正午时时角为0，以上公式可以简化为：
（3.14）
其中，h表示正午太阳高度角。
由两角和与差的三角函数公式,可得，
（3.15）
因此，
对于北半球而言，
；（3.16）
对于南半球而方，
；（3.17）
一地正午太阳高度随距太阳直射点的远近而变化。
（φ为某地的纬度，δ为太阳直射点的纬度，当地夏半年取+，冬半年取—）
（2）太阳高度的分布规律
①太阳直射点上的太阳高度为90°，晨昏线上的太阳高度等于0，昼半球的太阳高度大于0，夜半球的太阳高度小于0。
②正午太阳高度的分布规律。
一是纬度分布规律：
同一时刻，正午太阳高度由太阳直射点向南北两侧递减,二是季节分布规律：
夏至北回归线及其以北各纬度，正午太阳高度达到一年中的最大值；冬至日北半球各纬度，正午太阳高度达到一年中的最小值。
正午太阳高度角的计算。
首先要明确两个纬度：
太阳直射点的纬度和所求地的纬度，根据以下公式计算：
某地正午太阳高度=90°—（太阳直射点和该地间的纬度差，可用线段法分析计算）
4.3.5直杆所处的地点和日期求解
由于太阳不断运动着，太阳的运动轨迹可用太阳高度角（h）和太阳方位角（A）来描述，它们的计算公式如下：
（3.18）
（3.19）
式中，h为太阳高度角，t为时角，
为太阳的赤纬角，
为当地纬度，A为太阳的方位角。
当一地确定后，纬度便是定值，太阳倾角与日期与日期有关，t则决定于每日个时刻。
再将以上模型分别应用于附件2和附件3的影子顶点数据坐标，就能得出可能的地点与日期。
其次用MATLAB软件对附件2和附件3的影子顶点数据坐标的x轴和y轴数据进行处理，绘制出影子顶点轨迹图如图9.1和图9.2所示：
图9.1附件2影子顶点轨迹图
图9.2附件2影子顶点轨迹图
4.4问题四的求解
已知直杆的高度为2米，利用Matlabr2014b工具箱中的视频查看器播放视频获得截图，使用图片查看器中的测距工具测量视频中的直杆长度为680，从而得到视频中尺度和实际尺度的比例尺为340：
1。
从北京时间8:
54:
20每间隔3051帧（2分钟）截取图片，测量出对应时刻的视频影长l，利用Excel记录数据并计算出实际影长L列表如表3所示：
表3视频中所获的数据
北京时间
l（视频影长）
L（实际影长/米）
8:
54:
20
801.28
2.36
8:
56:
20
796.16
2.34
8:
58:
20
781.16
2.30
9:
00:
20
771.17
2.27
9:
02:
20
762.13
2.24
9:
04:
20
753.08
2.21
9:
06:
20
745.11
2.19
9:
08:
20
732.08
2.15
9:
10:
20
724.07
2.13
9:
12:
20
717.06
2.11
9:
14:
20
702.06
2.06
9:
16:
20
698.04
2.05
9:
18:
20
687.04
2.02
9:
20:
20
升级会员 