ABSTRACT

In this paper, aiming at the problem of sun shadow localization, by analyzing and calculating the parameters affecting the shadow length of the object, the rod length optimization model based on multi-parameter hierarchical search is established with the aim of inversion of the rod length, and the MATLAB simulation is used to realize the sun shadow positioning technology. Firstly, a straight-line shadow length variation model based on the change of solar height angle is established, and the relationship between the straight shadow length and the declination angle and the time angle is obtained. Secondly, the control variable is used to analyze the length of the straight rod. Using the method of control variable to analyze the variation rule of longitude, latitude, measurement date and measurement time and the curve of the straight shadow length from 9:00 to 15:00 (Beijing time), it is found that the shadow length decreases first and then increases at 11:45. The minimum value obtained is about 2.1. Finally, the model is tested by the actual measurement method, the error of the model result relative to the actual measured value is about 1.732%, and the model has higher precision. For the shadow localization problem, the relationship between the straight shadow length and the shadow vertex coordinate value is obtained through the geometric relationship. For the case where the straight rod location is unknown, this paper applies the long length of the shadow length and the solar elevation angle inversion at the same place. Based on the principle of the parameter, the pole length optimization model based on parameter hierarchical search is established, and the range of the poles in the same group is approximated to zero. The MATLAB is used to perform rough search and fine search on the longitude and latitude. The feasible longitude and latitude values are (108.4566˚E, 19.2771˚N) and (103.1670˚E, 3.2552˚S). The first group of latitude and longitude corresponds to the offshore area of Hainan Province, and its longitude and longitude error of Hainan Province is 1.256%. Therefore, the first location is designated as Hainan Province within the error tolerance. The second group of latitude and longitude corresponds to Indonesia.

Keywords:MATLAB Simulation, Sun Shadow Localization, Hierarchical Search, Single Target Optimization

基于MATLAB仿真的太阳影子定位模型

王者，张吉发，王壮

青岛大学电子信息学院，山东 青岛

收稿日期：2019年4月9日；录用日期：2019年4月19日；发布日期：2019年4月26日

摘 要

本文针对太阳影子定位问题，通过对影响物体影长的各个参数进行分析计算，以反演杆长最优为目标建立基于多参数分层搜索的杆长优化模型，利用MATLAB仿真从而实现太阳影子定位技术。首先建立基于太阳高度角变化的直杆影长变化模型，得到直杆影长与赤纬角、时角等变量的函数关系式；其次，运用控制变量的方法分析影长关于直杆长度、所处经度、纬度、测量日期以及测量时间的变化规律并作出2018年9时至15时(北京时间)的直杆影长变化曲线，发现影长随时间的推移先减小后增大，在11时45分取得最小值约为2.1 m；最后，本文通过实际测量的方法对模型进行检验，得到模型结果相对实际测量值的误差约为1.732%，模型有较高的精度。对于影子定位问题，通过几何关系得出直杆影长与影子顶点坐标数值的关系；针对直杆所在地点未知的情况，本文应用同一地点不同时刻的影长和太阳高度角反演得到的杆长相等的原理建立起基于参数分层搜索的杆长优化模型，以同一地点多组杆长值中极差近似为零为目标；利用MATLAB对经度和纬度先后进行粗搜索和精搜索，得到两组可行经纬值分别为：(108.4566˚E, 19.2771˚N)及(103.1670˚E, 3.2552˚S)；第一组经纬度对应位置为近海南省的海上，其经度与海南省经度误差为1.256%，所以本文在误差允许范围内将第一个地点定为海南省，第二组经纬度对应位置为印度尼西亚。

关键词 :MATLAB仿真，太阳影子定位，分层搜索，单目标优化

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

通过物体影子变化情况得到物体具体所在位置以实现“太阳影子定位”是定位分析非常重要的一个方面。某物体的影子与该物体长度以及所在位置的太阳高度角有关，太阳高度角主要由物体所处经纬度，日期等决定。本文建立物体影长变化模型，进而提出基于多参数分层搜索的杆长优化模型，进而利用MATLAB进行仿真得到较好的结果；本文利用实际测量当地某位置直杆影长的方法对模型进行检验。本文研究问题具有一定的实用价值，利用MATALB进行计算机仿真很好地提高了模型计算效率和精度，并且模型检验方法具有较好的操作性以及实际价值。

2. 太阳影长变化模型

2.1. 影长变化模型

本文以直杆作为研究载体探究直杆影长与各种因素的关系，直杆与其在地面投影的几何关系如图1所示：

图1. 直杆投影几何示意图

图2. 赤纬角和时角几何示意图

就固定在地面某一点的直杆而言，其投影点位置与太阳高度角有关，影子长度l为 $H/\tan \alpha$ ，其中，H为直杆长度， $\alpha$ 为太阳高度角。

地球某点处的太阳高度角 $\alpha$ 取决于该地的纬度 $\phi$ ，某确定时刻太阳赤纬角D以及太阳时角h，其几何关系如图2所示。太阳高度角解析关系式关系如式(1)所示：

$\alpha =\arcsin \left(\sin \phi \sin D+\cos \phi \cos D\cosh \right)$ (1)

首先，本文进行太阳赤纬角的推导 [1] ，赤纬角D有：

$\sin D=\sin {23}^{\circ }26^{\prime }37.8^{″}\sin \sigma$ (2)

其中 [2] ：

$\begin{array}{l}\delta =279.9348+d+1.914827\sin d-0.079525\cos d\\ +0.019938\sin 2d-0.001620\cos 2d\end{array}$ (3)

d为一年中特定日期表示的角分数，有：

$d=\frac{\left(n-1\right)\times 360}{365.242}$ (4)

其中，n为某日期的积日，积日为日期在年内的序列号，如1月1日的积日为1，以此类推。

其次，本文对太阳时角进行推导，时角由公式(5)导出：

$h=15\left(T-M\right)-longitude$ (5)

其中，T为计算时代入的GMT时间，与直杆所在点时间的换算如式(6)所示：

$T=t_B-\frac{longitude+120}{15}-8$ (6)

其中longitude为该点的经度，东经取负，西经取正； $t_B$ 为北京时间，即东八区标准时间，某地时间变化是由地球自转引起的 [1] ，由于地球一天24 h内旋转360˚，则每小时旋转15˚。M为真太阳时，由式(7)得 [3] ：

$M=12+0.123570\sin d-0.004289\cos d+0.153809\sin 2d+0.060783\cos 2d$ (7)

d的定义见式(4)。

至此，本文完成直杆影子长度变化模型的建立，联立公式即可得到影子长度与各个参数的函数关系式。

2.2. 各因素对影长的影响

本文建立直杆影子长度变化模型：直杆影子长度与直杆长度H，积日n，所在经度纬度以及测量时间(北京时间)有关。为分析影子长度关于各个参数的变化规律，本文以2018年10月22日北京时间9时至15时之间天安门广场(北纬39度54分26秒，东经116度23分29秒) 3米高的直杆的影长变化为例，采用控制变量的方法分别研究上述参数对影子长度变化的影响：

(I) 积日、直杆所在点经纬度和测量时间一定，直杆长度变化

当直杆长度变化，其他参数不变时，由公式1可知，太阳高度角 $\alpha$ 保持不变。则直杆影子长度与直杆自身长度成正比。基于此关系，本文作出图3研究直杆影子长度与杆长的关系，其中日期取2018年10月22日(即积日为295)，经度为东经116度23分29秒，纬度为北纬39度54分26秒，测量时间为12:00 (北京时间)，即有：

$t_B=12,\phi ={39.91}^{\circ },longitude=-{116.39}^{\circ },n=295$ 。

利用MATLAB作出直杆影子长度–直杆长度关系图如图3所示：

图3. 影子长度与直杆高度关系图

(II) 积日、直杆所在点经度值和测量时间一定，直杆所在点纬度值变化

当直杆所在点纬度值变化，其他参数不变时，根据一般性的认识，可以得到：在其他参数相同的情况下，纬度值越高，太阳高度角越小，直杆影子长度越大。其中日期取2018年10月22日(即积日为295)，经度为东经116度23分29秒，直杆长度为3 m，测量时间为12:00 (北京时间)，有：

$t_B=12;\phi =0^{\circ },{15}^{\circ },{30}^{\circ },{45}^{\circ },{60}^{\circ };longitude=-{116.39}^{\circ };n=295$

本文利用MATLAB作出直杆影子长度—纬度变化图如图4所示：

图4. 影子长度与所处纬度关系图

从图4可以发现，随着纬度值的增大，直杆影子长度一直增大，且增大的速度越来越快。

(III) 积日、直杆所在点纬度值和测量时间一定，直杆所在点经度值变化

时间变化是由地球的自转引起的，不同时区在同一太阳时所处的时间必然不同，例如当北京时间为9:00时，GMT时间为1:00，此时北京处于白天，格林尼治没有太阳照射，太阳高度角显然为0˚，此时的直杆影子长度趋于无穷大。则在其他参数相同的条件下，经度不同，所在地的白昼情况不同。日期取2018年10月22日(即积日为295)，纬度为北纬39度54分26秒，直杆长度为3 m，测量时间为12:00 (北京时间)，有：

$t_B=12;\phi ={39.91}^{\circ };longitude=-{20}^{\circ },-{30}^{\circ },-{60}^{\circ },-{80}^{\circ },-{100}^{\circ };n=295$

本文利用MATLAB作出直杆影子长度与直杆所在经度值关系图像如图5所示。

从图5可看出，在其他条件相同的前提下，东经0˚到40˚左右范围内，直杆影子长度都在0左右或是很大的异常值，则此地区在北京时间12:00正处于黑夜。为更精确反映二者关系，本文将经度取值范围缩小并作出图6。

从图中可以看出，在有光照的范围内，直杆影子长度随经度的增大而减小，而取得异常值的区域大致为东经30˚到东经40˚，根据世界时区与经度的关系，此区域在北京时间12:00时大概处于6点，在10月份正处于黑夜状态，符合实际。

(IV) 直杆所在点经纬度值和测量时间一定，积日变化

积日这一参数反映的是日期(即月日)对太阳高度角的影响，根据前文定义，本文发现在其他参数一定的条件下，太阳高度角随积日的增加应先增加后减小。且由于本文中的地点天安门广场处于北半球，则

图5. 影子长度与经度关系图

图6. 修正经度后的影长与经度关系图

其他参数一定的情况下，太阳高度角在夏至日(2018年6月21日)取得最大值，则此时直杆影子长度取最小值。利用MATLAB计算得到相应的直杆影子长度，其中纬度为北纬39度54分26秒，经度为东经116度23分29秒，直杆长度为3 m，测量时间为12:00 (北京时间)，有：

$t_B=12;\phi ={39.91}^{\circ };longitude=-{116.39}^{\circ };n=40,100,170,230,290$

本文利用MATLAB作出直杆影子长度–积日关系图如图7所示。

所作图7符合上文提到的影长随积日变化的趋势，本文研究图7得到直杆影子长度取最小时，积日为172。查阅资料得到今年2018年的夏至日为6月21日，对应的积日为173，则本文求得直杆影子长度最小时的积日误差大约为：

$\eta =\frac{n_{xia}-n_{\min }}{n_{xia}}=\frac{173-172}{173}\approx 0.578\%$ (8)

证明模型有相当好的精度。

图7. 影子长度与积日关系示意图

(V) 直杆所在点经纬度值和积日一定，测量时间变化

当其他参数一定时，即对于特定的一天(10月22日)，测量时间不同必然会对直杆影子长度产生影响。根据一般的生活经验，太阳高度角在早上和晚上很小，在中午一般会达到最大值，直杆影子长度会随着一天时间的推移先减小后增大。与上文类似，本文取日期为2018年10月22日，纬度为北纬39度54分26秒，经度为东经116度23分29秒，直杆长度为3 m，积日为295。利用MATLAB作出太阳高度角以及直杆影子长度随时间推移的变化示意图如图8，图9所示：

图8. 太阳高度角与时间关系示意图

从图9可以得到，直杆影子在11时45分左右取得最小值3.6560 m。

3. 太阳影子定位模型

3.1. 定位模型的建立

首先，本文以直杆底端为原点，水平地面为xy平面，作出直杆影子顶点与直杆关系的几何示意图如

图9. 影子长度与时间关系示意图

图10所示：

图10. 影子顶点几何示意图

则，根据图10可得直杆影子长度与x轴，y轴坐标值 $x_i$ 和 $y_i$ 的关系为：

$l_i=\sqrt{x_i^2+y_i^2}$ (9)

本节中本文将针对物体所处地理未知的情况下，根据影子长度变化(数据基于2015年高教社杯全国大学生数学建模竞赛A题)反推出物体所处地理位置。与直杆影子定位模型相比，该模型中直杆长度H，直杆所处地点的经度longitude和纬度 $\phi$ 未知，而任意一组经纬度 $\left[{\phi }_i,longitud{e}_i\right]$ 可定位至不同的直杆所处地点。从而根据上文影长变化模型可得任一北京时间 $t_{Bi}$ 时，经度为 $longitud{e}_i$ ，纬度为 ${\phi }_i$ 处的太阳高度角 ${\alpha }_i$ ，该时刻直杆的影长可以根据直杆影子顶点的坐标数据及公式求得，即可得到该地点处直杆的长度 $H_i$ ，即几个参量有如图11的关系。

基于如图11所示的参量关系，本文建立基于参数分层搜索的杆长优化模型：对于某一未知地点，本文设其经纬度值分别为 $longitud{e}_n$ ， ${\phi }_n$ ，并设测量影长顶点坐标数据时的北京时间数组为 $\left\{t_{Bi}\right\}$ 。

本文首先对经度和纬度进行搜索，在一组确定的经纬度值(即对应某一确定地点)前提下，再将测量时间 $t_{Bi}$ 代入模型，其中 $i=1,2,\cdots ,m$ 。则每一组经纬度值，都会有m个直杆的长度数据，将其存于数组 $A_{ganchang}$

图11. 参量关系示意图

中。由于每一组直杆长度数据都是在同一地点求得，则每一组直杆长度值都应相等，若这组数值差别较大，则说明搜索到的该经度纬度值不符合实际，则再进行下一组经度纬度值的搜索。本文设某一杆长数组中的最小值为 $H_{\min }$ ，最大值为 $H_{\max }$ ，以二者差值最小为优化目标：

$\left|H_{\max }-H_{\min }\right|<0.0005$ (10)

分层遍历搜索经度值，纬度值以及测量时间，对杆长数值进行优化，将每个经纬度下各个测量时间得到的直杆长度生成数组，判断该数组中的最大最小值是否满足公式(10)，若满足公式(10)，则将此时搜索到的经度纬度值存入最优经度数组 $A_{longitude}$ 以及最优纬度数组 $A_{\phi }$ ，将每一组对应的经纬度对应至地图中即可得出直杆所在的地点。最后，将附件1中的数据代入模型，即可得出若干可能的地点。

3.2. 定位模型的求解

首先，本文从数据中可得：

$\{\begin{array}{l}n=108\\ t_{Bi}=14.7:0.05:15.7\end{array}$ (11)

将式(11)代入基于参数分层搜索的杆长优化模型，其中 $m=21$ 。考虑到模型的实际应用以及精确性，经度搜索范围为东经。本文先对经纬度进行粗搜索，经度搜索范围为 $\left[-{180}^{\circ },0\right]$ ，步长为1˚，纬度搜索范围为 $\left[-{90}^{\circ },{90}^{\circ }\right]$ ，取负为南纬，取正为北纬，步长为1˚，时间取值 $t_{Bi}=14.7:0.05:15.7$ 。

针对此算法，本文作出模型求解算法流程图如图12所示。

经过如图12所示的经纬度粗搜索，本文得到符合条件的七组经纬度值如表1所示。

根据表1可以看出，符合条件的地点可能有两个：第一个潜在地点的纬度值范围为2˚S~4˚S，经度值范围为101˚E~106˚E；第二个潜在地点的纬度值范围为17˚N~21˚N，经度值范围为106˚E~111˚E。

为进一步求解两潜在地点的精确位置，本文对图12中的参数搜索算法进行进一步的改进：首先对潜在地点一进行精确定位，本文将经度取值范围修正为 $\left[-{106}^{\circ },-{101}^{\circ }\right]$ ，纬度取值范围修正为 $\left[-4^{\circ },-2^{\circ }\right]$ ，搜索步长均改为0.001˚，以式(10)为目标，利用MATLAB进行模型求解，求出符合条件的经纬度值共64,637组，将这些经纬度分别求平均值，得到潜在地点一的经纬度为(103.1670˚E, 3.2552˚S)，根据Google地图定位系统可以得出此地点在印度尼西亚。

对于潜在地点二，本文利用与潜在地点一相同的修正算法，将其经度范围修正为 $\left[-{111}^{\circ },-{106}^{\circ }\right]$ ，纬度范围修正为 $\left[{17}^{\circ },{21}^{\circ }\right]$ ，搜索步长均为0.001˚，利用MATLAB搜索求得符合条件的经纬度共39,180组，将经纬度分别求平均值得到潜在地点二的经纬度为(108.4566˚E, 19.2771˚N)，根据Google地图定位系统得出该地点位于近海南省的海上，而我国海南省的经度范围 [4] 为 (108.6167˚E, 111.0833˚E)，纬度范围为(18.1667˚N, 20.1667˚N)。潜在地点二的纬度值在海南省的纬度范围内，而经度值与海南经度值的误差为

图12. 影子定位算法流程图

表1. 符合条件的经纬度取值

0.1474%~2.3646%，说明所求地点的经度值与海南经度范围的误差值很小，并且问题二题干中说明此直杆固定在水平地面上，则在误差允许范围内，本文将潜在地点二定位于中国海南省，并且综合两地多组杆长的数据可以得出问题二中的直杆长约为2.1 m。

4. 模型检验

本文模型都是基于直杆影子长度变化模型建立的，所以本文对问题一中的模型进行模型检验，以验证本文所求结果的正确性。

由于直杆影长的测量具有较好的操作性，所以本文利用实际测量的方法对太阳影子长度变化模型进行检验。首先，本文查阅到实测地点的经纬度为(120.3849˚E, 36.0620˚N)，测量日期为2018年9月11日，则对应积日 $n^{\prime }=254$ ，直杆则为长度为1.8 m的小旗杆。由于比赛时间紧张，本文只对9:30，10:47，12:10，14:30，16:00 (上述时间均为北京时间)时的小旗杆影子长度进行测量，测量工具是分度值为1 mm的米尺，所以影长实测值保留4位小数。将上述参数依次代入问题一中的太阳影子定位模型求得影子长度，则五个时间点对应的实际影长测量值和模型求解值见表2：

表2. 实测与模型结果比对

从表2可以看出，影长变化中太阳影子定位模型具有相当好的精度，在其他地点影长的测量也有较好的推广性。

5. 稳定性分析

由于本文建立的模型均是以影长变化模型为基础，所以本文对影长变化模型进行稳定性分析，分别对影长变化中天安门广场纬度值分别作±3%的干扰，分析这纬度因素的波动对影长变化模型求解的直杆影长的影响。

首先，本文对影长变化中的纬度值作±3%的扰动，并取北京时间10:00，12:00以及14:00的直杆影长值如表3所示：

表3. 纬度稳定性分析

从表3中可以看出影长变化的模型对于纬度波动具有较好的稳定性，为更清楚地看出影长变化模型对于纬度波动的稳定性，本文作出纬度稳定性分析图如图13所示。

6. 总结

本文提出的基于多参数分层搜索的杆长优化模型具有较高的使用推广价值，对应一般性的太阳影子

图13. 纬度波动时模型稳定性分析图

定位问题都有较高的求解效率；并且本文对优化模型中优化目标的表示方法对于同一参数有多个求解结果的问题也有着很好的启发性价值；同时，本文利用实际测量当地某位置处直杆影长的方法对模型进行检验，这种检验方法具有较好的操作性以及实际价值，此方法具有一定的参考价值。

文章引用

王 者,张吉发,王 壮. 基于MATLAB仿真的太阳影子定位模型
Sun Shadow Positioning Model Based on MATLAB Simulation[J]. 计算机科学与应用, 2019, 09(04): 790-801. https://doi.org/10.12677/CSA.2019.94090

参考文献