太阳影子定位
摘要
本文主要针对不同地点、不同时间影子变化的问题，建立数学模型，进行量化求解。同时对影响影子长度变化的各项因素进行分析。
对于问题一，首先确定太阳赤纬角、时角和太阳高度角等参数，并引入他们的关系式，再依据杆长、影长和太阳高度角间的三角函数关系，初步建立影长变化模型。查阅相关文献，建立太阳高度角模型。因经纬度与物体所处的位置相关，时角与时刻相关，太阳赤纬角与日期相关，可知位置、日期和时刻三个因素同时对太阳高度角产生影响，从而影响影子长度的变化。对于问题一中所给出的天安门的相关信息，可以得到该地区影子长度在北京时间9:00-15:00之间，出现先减小后增加的非线性变化规律，并且在北京时间12:06，影长达到最小，为3.8440米。
对于问题二，已知某固定直杆在水平地面上一段时间内的太阳影子顶点坐标数据，要求确定直杆所在的地理位置。首先建立单目标优化模型，由于影长和杆长为正比关系，因此将直杆不同时刻影长比的真实值与测量值的差平方和作为目标函数，然后采用网格搜索法，以太阳高度角大于零和太阳照射时间段在6:00-18:00之间为约束，确定经度和纬度的大致搜索范围，再进一步进行精确搜索，最后将目标函数值取到最小时所对应的经纬度，作为物体所在地的经纬度。将模型应用于附件1，求得可能地点为：云南省安沛市(104.5°E,22°N) ，海南省东方市(108.6°E,19.3°N) 。
对于问题三，由于影响目标函数的决策变量较多（如经度，纬度，时角，杆长，日期），因此引入影子方向角概念。首先建立关于影长比的优化模型，以不同时刻影长比真实值与测量值的累计误差作为目标函数。然后在此基础上，建立方向角差模型，从而建立多目标优化模型。在求解时以经纬度和日期作为决策变量，代入影长比的优化模型，再将其结果代入方向角差的优化模型，求出最优解。附件2最终确定结果为7月25日新疆维吾尔自治区喀什(79.56°E,39.07°N) ，附件3最终确定结果为2月3日湖北省十堰市(110.19°E,31.99°N) 。
最后，对模型中运用的方法进行了分析，并讨论了模型优缺点。考虑了模型在实际运用中的改进方向，使模型具有较强的通用性和可塑性。
关键词：太阳高度角；最小二乘法；网格搜索法；方位角差模型；多目标优化模型。
一、问题的重述
确定视频的拍摄地点和拍摄日期是视频数据分析中非常重要的方面，太阳影子定位技术，就是通过分析视频中物体在太阳下的影长和影子方位的变化，确定该视频拍摄的地点和日期的一种方法。
1. 建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用所建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒, 东经116度23分29秒）上3米高的直杆的太阳影子长度的变化曲线。
2. 根据某固定直杆在水平地面上不同时刻的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将模型应用于附件1的影子顶点坐标数据（已知拍摄日期），给出若干可能的地点。
3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干可能的地点与日期。
二、问题分析
2.1问题一的分析
对于问题一，要求建立确定日期、地点和杆长的影子随时间变化的数学模型。首先引入太阳高度角、赤纬角、时角、纬度、真太阳时和平太阳时的概念，再根据影长、杆长、太阳高度角三者之间的三角函数关系，建立影长关于太阳高度角和杆长变化的初步模型。然后通过查阅相关文献，得到太阳高度角关于赤纬角、时角和纬度计算公式以及赤纬角关于日期序号的计算公式。由于问题一所给的时间为平太阳时，而真太阳时与平太阳时存在时间差，因此引入真太阳时与平太阳时之间的换算公式，从而建立太阳高度角模型。最后将太阳高度角模型代入所建立的影长变化初步模型，得到最终的影长变化模型。此外，根据影子与其他因素的关系图像，分析其长度关于日期、经纬度、时刻、杆长变化的规律。
2.2问题二的分析
对于问题二，首先建立影长关于经纬度变化的数学模型。由于已知太阳高度角的函数关系而此问题中杆长为未知量，通过问题一所得影长与杆长所成的正比关系，将相邻两个时间点的影子长度相除即可消去未知量杆长。然后采用网格搜索法，通过建立太阳高度角和日照时间的约束条件，确定经度和纬度的大致搜索范围，运用问题一所得的影长关系，求出理论影子的长度。定义理论影长与实际影长的比例，
取21种时刻为一组，求出每组累计误差。然后将误差大小作为指标，得到初步解。再进一步进行精确搜索，得到多组可能位置。将模型运用于处理后的附件1数据，由于得到的结果存在不唯一性，要考虑对结果进行检验和选择。
2.3问题三的分析
对于问题三，与问题二相比，在其基础上增加了日期变量。为了简化方程的求解，可以采用网络搜索方法。首先对问题二所建立的网格搜索模型进行优化。添加日期变量。同时，由于问题三变量的增加，为进一步提高求解的准确性，引入影子方向角的概念，即用两个时间点之间影子的旋转角度作为判别依据，建立方向角差模型。之后结合补充日期变量后的影长变化优化模型，最终建立多目标优化模型。求解时使经纬度逐步以1，0.1度为间隔进行搜索，然后将所得结果带入方向角差优化模型进行进一步精确，将经纬度的间隔分别细化为0.01度，再次进行搜索，得出最优解。
三、模型假设
1. 由于太阳距地球很远，假设太阳为点光源，忽略其大小与体积；
2. 大气云层会导致太阳光折射，假设太阳光不会产生大气折射现象；
3. 海拔对影子存在影响，假设该地所处位置海拔为0；
4. 周围建筑物可能对太阳光存在遮挡，假设太阳光不被遮挡。
四、符号说明符号
H
w
δ
L
l
Φ
ψ
T
Tc
∆T
n
i
ω
ρ意义太阳高度角太阳时角赤纬角直杆的长度影子的长度纬度经度真太阳时平太阳时时间差一年中的第n天第i组数据影子方位角点对误差
五、模型的建立与求解
太阳影子定位技术是确定视频拍摄的地点和日期的方法，一天中，阳光下物体影子的变化遵循一定客观规律，依据此规律建立相关数学模型，求解目标函数，得到结果，如直杆所在的地点或拍摄的日期和时间。
5.1问题一的模型建立与求解
5.1.1模型的准备
明确以下概念：
真太阳时：太阳视圆面中心连续两次上中天的时间间隔。1真太阳日又分为24真太阳时。这个时间系统称为真太阳时。为了和日常生活习惯一致，把真太阳时定义为：真太阳视圆面中心的时角加12小时。
平太阳时：天文学上假定由一个太阳(平太阳) 在天赤道上作等速运行，其速度等于运行在黄赤道上真太阳的平均速度。
平时所用的北京时间是平太阳时，如果考虑地球绕太阳的旋转轨道是椭圆，自转非均匀的，此时得到的时间为真太阳时。真太阳时要求每天的中午12点，太阳处于头顶的最高位置。真太阳时与平太阳时之差称为真平太阳时差。可得公式：
T=Tc+∆T
(1)
其中，T为真太阳时，Tc为平太阳时，∆T为时间差。
下文提到的时间都是指真太阳时，如果考虑平太阳时，只需在相应的真太阳时上减去时间差即可。
时角：一个天体的时角表示该天体是否通过了当地的子午圈。其数值表示该天体与当地子午圈的角距离，并借用时间的单位，以小时来计量，一小时为15度。通常以当地时间12点的时角作为零，前后每隔一小时，增加15度。可得公式：
w=15×(T−12)
其中，w为时角。
5.1.2模型的建立
Step1:影长变化模型的初步建立
为了建立合理的数学模型，需要分析固定直杆在太阳下的成影情况，构建的示意图如下
：(2)
图1:直杆成影示意图
如上图所示，假设直杆高为L，直杆顶点P与地面杆影顶点P′所连成的直线PP′与地面所成的角为太阳高度角H，直杆影长为l。根据三角函数关系，可得影长l的数学表达式为：
l=其中，
tanH=将(4)式代入(3)式, 即得：
√Ll=sinH
(5)sinH(4)LtanH(3)
Step2:太阳高度角模型的建立
太阳高度角是观察太阳时的仰角，也就是太阳光线与地面之间的夹角。由于地球的自转，造成太阳东升西落，太阳高度角在一日内不断发生变化。通过查阅相关文献，假设H为太阳高度角，Φ为纬度，δ为赤纬角，w为太阳时角，得到太阳高度角的计算公式如下：
sinH=sinΦsinδ+cosΦcosδcosw
其中，赤纬角δ与日期序号n之间的关系式如下：
δ=23. 45×sin(2π×
Step3:影长变化模型的完善
将太阳高度角模型代入影长变化模型，即将公式(6)代入公式(5)，得到：
√︀Ll=sinΦsinδ+cosΦcosδcosw
5.1.3模型的求解
Step1:赤纬角δ的求解
问题一中所给日期为2015年10月22日，经换算得该日为一年中的第295天，即n=295, 代入下式
δ=23. 45×sin(2π×
得：δ=-12.1017°
Step2:时角的求解
将北京时间9:00-15:00该时间段内的时间进行转化分解。如10时30分转化为10.5时，13点20分转化为13.3时，以0.1时为单位时长，可分为60个时段，然后利用公式(2)求得各个时间段的时角。
w=15×(T−12)
Step3:太阳高度角的求解
将求解得到的赤纬角δ，时角w，当地纬度Φ代入公式(6)，求得太阳高度角的正弦值：
sinH=sinΦsinδ+cosΦcosδcosw
如北京时间9:00时太阳高度角的正弦值sinH为0.3958。
Step4:影长的求解
根据问题一已知条件，得到杆长L=3m, 将其代入公式(5),即得影长l。284+n) 365(8)284+n) 365(7)(6)
用MATLAB软件进行编程，得到影长l随时间的变化曲线图，如下所示
：
图2:影长变化曲线图
记录上图中一些特殊时刻的影长大小与其对应时间点，如下所示：
表1:特殊点影长
表
可见，北京时间9：00-15:00之间，杆影长度呈现先减小后增大的非线性变化，在12：06左右，影子长度达到最小，符合日常生活实际情况。
Step5:影长变化规律的求解
1. 影长与杆长的关系
图3:影长与杆长的关系图
由图可知，同一地点，同一日期，影长和杆长成正比关系。
2. 影长与日期的关
系
图4:影长与日期的关系图
由图可知，以1月1日为一年中的第一天，全年中影长呈现先减小后增大的非线性变化，且在夏至日左右，影长达到全年最小，冬至日左右影长达到全年最长。
3. 影长与经度的关系
图5:影长与经度的关系图
由图可知，同一纬度，在不考虑建筑物遮挡太阳光的情况下，10月22日正午时分，经度0°-41°E左右的地区太阳高度角小于0，故没有出现影长，经度41°E-51°E左右的地区由于太阳高度角极小，故影长较长，经度为45°E左右时，理论上3m高的直杆，其影长最高可达350m。经度51°E-180°E左右的地区，影长较平均，一般为20m以下。
4. 纬度与经度的关
系
图6:纬度与经度的关系图
由图可知，同一经度，在不考虑建筑物遮挡太阳光的情况下，北半球的10月22日正午时分，纬度0°-64°N左右的地区太阳高度角较大，故影长较短，
经度64°N-82°N左右的地区由于太阳高度角极小，故影长较长，理论上3m高的直杆，其影长最高可达190m。高纬度地区太阳高度角为0°，故影长为0。
5.2问题二的模型建立与求解
5.2.1模型的准备
理论影长和实际影长的定义
在一段时间内连续测得某固定直杆的影子横坐标x和纵坐标y，其中测量时间点(换算成小时为单位) 为ti,
ti=t0+i∆t(i=0, 1, 2,...m)
式中i为测量点数，m+1为测量次数，t0为初始测量时间点，∆t为测量时间间隔。
实际影长为：
′li√︀=(xi) 2+(yi) 2(9)
式中xi, yi分别为ti时刻测得的横纵坐标。
理论影长可由得影长变化模型得：
√︀Ll=sinΦsinδ+cosΦcosδcosw
式中δ为赤纬角，ψ为经度，w为时角。且有
wi=π(ti−12) 12(10)
(11)ψ=120−15(t′0−t0)
式中t0为当地时间，t′0为t0对应的北京时间。
5.2.2模型的建立
Step1:单目标优化模型的建立
为求得纬度Φ和赤纬角δ使理论影长l无限接近于实际影长l′，利用最小二乘法建立非线性优化模型：n∑︁′2min(li−li) (12)
⎧⎨sinH>0
⎩6
太阳始终在地平线以上，太阳高度角大于零，则其正弦值必大于零。且默认太阳在地平线以上的地方时在6点到18点之间。则有以使l和l′差值的平方最小为目标，以太阳高度角大于零为约束条件，建立初步的非线性优化模型。
Step2:影长比模型的建立
在上步的优化模型中，理论影长l与杆长L有关，目标函数中有杆长L、赤纬角δ和纬度Φ三个变量，不易求解，故需对该模型改进。由于理论影长l与杆长L成长比，将相邻两次li相除，可抵消L，故得到以下模型：
n′∑︁li+1li
min(−+1) 2
lili
i=0
(14)
⎧
⎨sinH>0⎩6
目标函数中,
tan Hili+1=litan Hi+1
(15)
化简后可见目标函数只与赤纬角δ和纬度Φ两个变量有关。
Step3:网格搜索
′
①在所需空间范围内，选取一个经纬位置[ψ,Φ]，逐一计算21次测量点li和li差
值的平方，并进行累加，得到点对误差ρij；
②逐一选取空间范围内的所有经纬位置，比较得到点对误差ρij最小的空间位置。
5.2.3模型的求解
附件1中已知2015年4月18日北京时间14：42至15：42每间隔3分钟测得的某固定直杆影子横纵坐标，需确定直杆所处的空间位置。
Step1:数据的处理
①真太阳时(即地方时) 的换算：初始时间点t0为14.7时，时间间隔∆t为0.05小时，其中测量次数m+1为21次。
②实际影长l′的计算：将21组xi，yi代入下式：
′
li
′
得到对应的21个li：
√︀
=(xi) 2+(yi) 2
Step2:定位空间模型的求解
′
将各测量时间点的实际影长li代入太阳影子定位空间模型，纬度δ满足
|δ|≤90
根据附件一数据可得测量时间为当地下午，地方时满足12≤T≤18, 由公式(11)可知，此条件下经度ψ应满足
19. 8≤ψ≤94. 8
在此范围内进行网格搜索。得到满足约束条件的8个区块的目标解(纬度均为北纬，经度均为东经) ：
例如在纬度[0,30]，经度[105,115]范围内以间隔为0.01搜索的结果如下图：
图7:搜索结果图
找出图像在z轴上最低的点为[36,193],即实际经纬度为[108.6°E,19.3°N]。发现仅有两个空间位置的点对误差相对其他位置极小，数量极为10−4。故认为附件1测量的固定直杆可能在2015年4月18日位于以下空间位置(纬度均为北纬，经度均为东经) ：
5.3问题三的模型建立与求解5.3.1模型的准备
以地球中心为原点，赤道平面为xy平面，指向北极的方向为z轴正向，
固定某地S，正午时刻该地所在大圆的平面为xz平面，建立直角坐标系，如下图所示：
图8:地球地标系示意图
通过查阅相关文献，得到上图中正北方向N旋转至太阳影子的方向的角度（即影子方位角）ω满足：
cosω=
cosδcoswsinφ−sinδcosφ
cosHsinω=
sinwcosδcosH
(16)(17)
同时给出cosω和sinω是为了避免出现方位角多解的情况。5.3.2模型的建立
按问题二给出的网格搜索法求解，得到每一天的最优解及相应的优化变量值，然后在得到的365个最优值中选取最优解。因为单纯采用太阳影子的长度信息进行计算，结果并不是最理想，所以考虑在目标函数中利用已知信息进行优化，即引入影子方位角。
Step1:方位角差模型的建立
由于无法确定具体的地理方位，因此如果直接采用方位角作为判别依据，就会增加参数（坐标、旋转角度）。为了优化参数，可以用方位角的差，也就是用两个时间点之间实际影子的旋转角度与理论影子的旋转角度的差方和作为判别依据。由于相邻点之间旋转角误差极小，故取相互间隔为10的点作差，得到公式如下：
min
10∑︁i=0
′′2
[(ωi+10−ωi) −(ωi+10−ωi)]
(18)
Step2:多目标优化模型的建立
在已有问题二关于影长比的单目标优化模型基础上补充日期变量，结合上述方位角差模型，建立多目标优化模型：
365∑︁20′∑︁li+1limin(−+1) 2
lili
n=0i=0365∑︁10∑︁n=0i=0
(19)
min
5.3.3模型的求解
′′2
[(ωi+10−ωi) −(ωi+10−ωi)]
(20)
在问题二中网格搜索法的基础上，增加对于日期范围的设定，确定方位角搜索的范围，搜索流程如下图所示
：
图9:网格搜索法流程图
Step1:影长比优化模型的求解
对于公式(19)，以附件2为例，使经纬度分别以1度为间隔，进行网格搜索，得出结果如下图所示
：
图10:网格搜索法流程图
Step2:方位角优化模型的求解
将所得结果带入公式(20)进行进一步优化，将经纬度的间隔分别细化为0.01度，再次进行搜索，得出最优解，如下图所示：
下面就题目所给的附件2和附件3数据，分别搜索得到以下搜索结果：
表2:附件2
日期7月25日5月26日5月27日
经度79.56°92.01°92.01°
纬度39.07°57.01°57.01°
地点
新疆维吾尔自治区喀什地区
俄罗斯南部俄罗斯南部
表3:附件3
日期2月3日
经度110.19°
纬度31.99°
地点湖北省十堰市
六、模型评价与推广
6.1模型的优点
1. 通过最小二乘法，解决了超定方程无法求出解析解的问题。2. 利用太阳高度角的限定条件和影长变化规律缩小了网格搜索范围。6.2模型的缺点
1. 不考虑太阳光的折射，计算工程量较大。
6.3模型的推广
该模型可以根据一般建筑物的影长变化，较精确地判断该建筑物的所在地以及拍摄日期和时间。
基于该模型可以求解任意建筑物一年中被其他建筑物影响光照的时间，
参考文献
[1]蔡志杰，太阳影子定位，数学建模及其应用，第4期，2015。
[2]孙吉山，洪薇．北京标准时间与各地真太阳时的换算法．上海针灸杂志，1994。[3]司守奎，孙兆亮．数学建模算法与应用．2版．北京：国防工业出版社，2015。
附录一：问题一的MATLAB 程序 d=295; %日期的天数
a=23.45*sind(360*(284+d)/365); %太阳赤纬 w=39.91;%纬度 l=3;%杆长 j=0; x=[]; y=[];
for i=9:0.1:15 %真太阳时 j=j+1;
b=sind(w)*sind(a)+cosd(w)*cosd(a)*cosd(15*(i-12));%求出太阳高度角的正弦值sinH c=b/sqrt(1-b^2);%求出太阳高度角的正切值tanH; ly=l/c;%影长 x(j,1)=ly; end
附录二：问题二的MATLAB 程序 d=108; %日期的天数
a=23.45*sind(360*(284+d)/365); x=[];
y=[1.149625826 1.182198976 1.215296955 1.249051052 1.28319534 1.317993149 1.353364049 ...
1.389387091 1.426152856 1.463399853 1.501481622 1.540231817 1.579853316 1.620144515...
1.661270613 1.703290633 1.74620591 1.790050915 1.835014272 1.880875001 1.927918447]; k=0; Z=[];
for w=0:0.1:30%控制纬度 k=k+1; t=0;
for jing=79.5:0.1:99.5%控制经度 t=t+1; z=0; j=0;
for e=14.7:0.05:15.7 %真太阳时 i=e-(120-jing)/15;%求出当地时间 j=j+1;
b=sind(w)*sind(a)+cosd(w)*cosd(a)*cosd(15*(i-12));%求出太阳高度角的正弦值sinH if b>0
c=b/sqrt(1-b^2);%求出太阳高度角的正切值tanH x(j,1)=1/c; if e>=14.75
z=((y(1,j)/y(1,j-1))-(x(j,1)/x(j-1,1)))^2+z;
Z(k,t)=z;%求累计差值 end else
Z(k,t)=10000; end end end end
[x,y]=find(Z==min(min(Z)))%寻找最小值所在的经纬度 附录三：问题三的MATLAB 程序 （1）确定方位角范围 x=[];
y=[1.247256205 1.22279459 1.198921486 1.175428964 1.152439573 1.1299174... 1.10783548 1.086254206 1.065081072 1.044446265 1.024264126 1.004640314... 0.985490908 0.966790494 0.948584735 0.930927881 0.91375175 0.897109051... 0.880973762 0.865492259 0.850504468]; Z=[]; k=0; o=0; daan=[];
for d=0:1:73%控制日期 o=o+1;
a=23.45*sind(360*(284+d)/365); for w=0:1:90%控制日期 k=k+1; t=0;
for jing=19.8:1:94.8%控制经度 j=0; t=t+1; z=0;
for e=12.68:0.05:13.68 i=e-(120-jing)/15; j=j+1;
b=sind(w)*sind(a)+cosd(w)*cosd(a)*cosd(15*(i-12));%求出太阳高度角的正弦值sinH
if b>0
c=b/sqrt(1-b^2);%求出太阳高度角的正切值tanH x(j,1)=1/c; if e>=12.73
z=((y(1,j)/y(1,j-1))-(x(j,1)/x(j-1,1)))^2+z; Z(k,t)=z;%求累计差值 end else
Z(k,t)=10000;
end end end
[u,v]=min(min(Z));%寻找当天搜索范围内累计差值最小的点 daan(o,1)=u;%输出每日的最小累计差值 end
（2）方位角差方和最小 x=[];
y=[-82.02711405 -81.10849053 -80.16448976 -79.19768875 -78.20352191 -77.19123876... -76.14438496 -75.07333958 -73.9762777 -72.84887252 -71.68887638 -70.4973275 ... -69.27179577 -68.0170072 -66.72685728 -65.39508527 -64.03225482 -62.63232339... -61.18607201 -59.70496474 -58.18242197]; Z=[]; k=0; o=0; daan=[];
for d=120:1:180%控制日期 o=o+1;
a=23.45*sind(360*(284+d)/365); for w=35:0.1:45%控制日期 k=k+1; t=0;
for jing=80:0.1:85%控制经度 j=0; t=t+1; z=0;
for e=12.68:0.05:13.68 i=e-(120-jing)/15; j=j+1;
b=sind(w)*sind(a)+cosd(w)*cosd(a)*cosd(15*(i-12));%求出太阳高度角的正弦值sinH
if b>0
c=sqrt(1-b^2);
C=sind(a)*cosd(15*(i-12))/c;%求出影长变化的正切值 da=asind(C);%求出角度 x(j,1)=da; if e>=13.18
z=((y(1,j)-y(1,j-1))^2-(x(j,1)-x(j-1,1))^2)^2+z; Z(k,t)=z;%求出方差和 end else
Z(k,t)=10000;
end end end
[u,v]=min(min(Z));%寻找当日搜索范围内最小的方差和 daan(o,1)=u;输出每日的最小方差和 end
太阳影子定位
摘要
本文主要针对不同地点、不同时间影子变化的问题，建立数学模型，进行量化求解。同时对影响影子长度变化的各项因素进行分析。
对于问题一，首先确定太阳赤纬角、时角和太阳高度角等参数，并引入他们的关系式，再依据杆长、影长和太阳高度角间的三角函数关系，初步建立影长变化模型。查阅相关文献，建立太阳高度角模型。因经纬度与物体所处的位置相关，时角与时刻相关，太阳赤纬角与日期相关，可知位置、日期和时刻三个因素同时对太阳高度角产生影响，从而影响影子长度的变化。对于问题一中所给出的天安门的相关信息，可以得到该地区影子长度在北京时间9:00-15:00之间，出现先减小后增加的非线性变化规律，并且在北京时间12:06，影长达到最小，为3.8440米。
对于问题二，已知某固定直杆在水平地面上一段时间内的太阳影子顶点坐标数据，要求确定直杆所在的地理位置。首先建立单目标优化模型，由于影长和杆长为正比关系，因此将直杆不同时刻影长比的真实值与测量值的差平方和作为目标函数，然后采用网格搜索法，以太阳高度角大于零和太阳照射时间段在6:00-18:00之间为约束，确定经度和纬度的大致搜索范围，再进一步进行精确搜索，最后将目标函数值取到最小时所对应的经纬度，作为物体所在地的经纬度。将模型应用于附件1，求得可能地点为：云南省安沛市(104.5°E,22°N) ，海南省东方市(108.6°E,19.3°N) 。
对于问题三，由于影响目标函数的决策变量较多（如经度，纬度，时角，杆长，日期），因此引入影子方向角概念。首先建立关于影长比的优化模型，以不同时刻影长比真实值与测量值的累计误差作为目标函数。然后在此基础上，建立方向角差模型，从而建立多目标优化模型。在求解时以经纬度和日期作为决策变量，代入影长比的优化模型，再将其结果代入方向角差的优化模型，求出最优解。附件2最终确定结果为7月25日新疆维吾尔自治区喀什(79.56°E,39.07°N) ，附件3最终确定结果为2月3日湖北省十堰市(110.19°E,31.99°N) 。
最后，对模型中运用的方法进行了分析，并讨论了模型优缺点。考虑了模型在实际运用中的改进方向，使模型具有较强的通用性和可塑性。
关键词：太阳高度角；最小二乘法；网格搜索法；方位角差模型；多目标优化模型。
一、问题的重述
确定视频的拍摄地点和拍摄日期是视频数据分析中非常重要的方面，太阳影子定位技术，就是通过分析视频中物体在太阳下的影长和影子方位的变化，确定该视频拍摄的地点和日期的一种方法。
1. 建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用所建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒, 东经116度23分29秒）上3米高的直杆的太阳影子长度的变化曲线。
2. 根据某固定直杆在水平地面上不同时刻的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将模型应用于附件1的影子顶点坐标数据（已知拍摄日期），给出若干可能的地点。
3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干可能的地点与日期。
二、问题分析
2.1问题一的分析
对于问题一，要求建立确定日期、地点和杆长的影子随时间变化的数学模型。首先引入太阳高度角、赤纬角、时角、纬度、真太阳时和平太阳时的概念，再根据影长、杆长、太阳高度角三者之间的三角函数关系，建立影长关于太阳高度角和杆长变化的初步模型。然后通过查阅相关文献，得到太阳高度角关于赤纬角、时角和纬度计算公式以及赤纬角关于日期序号的计算公式。由于问题一所给的时间为平太阳时，而真太阳时与平太阳时存在时间差，因此引入真太阳时与平太阳时之间的换算公式，从而建立太阳高度角模型。最后将太阳高度角模型代入所建立的影长变化初步模型，得到最终的影长变化模型。此外，根据影子与其他因素的关系图像，分析其长度关于日期、经纬度、时刻、杆长变化的规律。
2.2问题二的分析
对于问题二，首先建立影长关于经纬度变化的数学模型。由于已知太阳高度角的函数关系而此问题中杆长为未知量，通过问题一所得影长与杆长所成的正比关系，将相邻两个时间点的影子长度相除即可消去未知量杆长。然后采用网格搜索法，通过建立太阳高度角和日照时间的约束条件，确定经度和纬度的大致搜索范围，运用问题一所得的影长关系，求出理论影子的长度。定义理论影长与实际影长的比例，
取21种时刻为一组，求出每组累计误差。然后将误差大小作为指标，得到初步解。再进一步进行精确搜索，得到多组可能位置。将模型运用于处理后的附件1数据，由于得到的结果存在不唯一性，要考虑对结果进行检验和选择。
2.3问题三的分析
对于问题三，与问题二相比，在其基础上增加了日期变量。为了简化方程的求解，可以采用网络搜索方法。首先对问题二所建立的网格搜索模型进行优化。添加日期变量。同时，由于问题三变量的增加，为进一步提高求解的准确性，引入影子方向角的概念，即用两个时间点之间影子的旋转角度作为判别依据，建立方向角差模型。之后结合补充日期变量后的影长变化优化模型，最终建立多目标优化模型。求解时使经纬度逐步以1，0.1度为间隔进行搜索，然后将所得结果带入方向角差优化模型进行进一步精确，将经纬度的间隔分别细化为0.01度，再次进行搜索，得出最优解。
三、模型假设
1. 由于太阳距地球很远，假设太阳为点光源，忽略其大小与体积；
2. 大气云层会导致太阳光折射，假设太阳光不会产生大气折射现象；
3. 海拔对影子存在影响，假设该地所处位置海拔为0；
4. 周围建筑物可能对太阳光存在遮挡，假设太阳光不被遮挡。
四、符号说明符号
H
w
δ
L
l
Φ
ψ
T
Tc
∆T
n
i
ω
ρ意义太阳高度角太阳时角赤纬角直杆的长度影子的长度纬度经度真太阳时平太阳时时间差一年中的第n天第i组数据影子方位角点对误差
五、模型的建立与求解
太阳影子定位技术是确定视频拍摄的地点和日期的方法，一天中，阳光下物体影子的变化遵循一定客观规律，依据此规律建立相关数学模型，求解目标函数，得到结果，如直杆所在的地点或拍摄的日期和时间。
5.1问题一的模型建立与求解
5.1.1模型的准备
明确以下概念：
真太阳时：太阳视圆面中心连续两次上中天的时间间隔。1真太阳日又分为24真太阳时。这个时间系统称为真太阳时。为了和日常生活习惯一致，把真太阳时定义为：真太阳视圆面中心的时角加12小时。
平太阳时：天文学上假定由一个太阳(平太阳) 在天赤道上作等速运行，其速度等于运行在黄赤道上真太阳的平均速度。
平时所用的北京时间是平太阳时，如果考虑地球绕太阳的旋转轨道是椭圆，自转非均匀的，此时得到的时间为真太阳时。真太阳时要求每天的中午12点，太阳处于头顶的最高位置。真太阳时与平太阳时之差称为真平太阳时差。可得公式：
T=Tc+∆T
(1)
其中，T为真太阳时，Tc为平太阳时，∆T为时间差。
下文提到的时间都是指真太阳时，如果考虑平太阳时，只需在相应的真太阳时上减去时间差即可。
时角：一个天体的时角表示该天体是否通过了当地的子午圈。其数值表示该天体与当地子午圈的角距离，并借用时间的单位，以小时来计量，一小时为15度。通常以当地时间12点的时角作为零，前后每隔一小时，增加15度。可得公式：
w=15×(T−12)
其中，w为时角。
5.1.2模型的建立
Step1:影长变化模型的初步建立
为了建立合理的数学模型，需要分析固定直杆在太阳下的成影情况，构建的示意图如下
：(2)
图1:直杆成影示意图
如上图所示，假设直杆高为L，直杆顶点P与地面杆影顶点P′所连成的直线PP′与地面所成的角为太阳高度角H，直杆影长为l。根据三角函数关系，可得影长l的数学表达式为：
l=其中，
tanH=将(4)式代入(3)式, 即得：
√Ll=sinH
(5)sinH(4)LtanH(3)
Step2:太阳高度角模型的建立
太阳高度角是观察太阳时的仰角，也就是太阳光线与地面之间的夹角。由于地球的自转，造成太阳东升西落，太阳高度角在一日内不断发生变化。通过查阅相关文献，假设H为太阳高度角，Φ为纬度，δ为赤纬角，w为太阳时角，得到太阳高度角的计算公式如下：
sinH=sinΦsinδ+cosΦcosδcosw
其中，赤纬角δ与日期序号n之间的关系式如下：
δ=23. 45×sin(2π×
Step3:影长变化模型的完善
将太阳高度角模型代入影长变化模型，即将公式(6)代入公式(5)，得到：
√︀Ll=sinΦsinδ+cosΦcosδcosw
5.1.3模型的求解
Step1:赤纬角δ的求解
问题一中所给日期为2015年10月22日，经换算得该日为一年中的第295天，即n=295, 代入下式
δ=23. 45×sin(2π×
得：δ=-12.1017°
Step2:时角的求解
将北京时间9:00-15:00该时间段内的时间进行转化分解。如10时30分转化为10.5时，13点20分转化为13.3时，以0.1时为单位时长，可分为60个时段，然后利用公式(2)求得各个时间段的时角。
w=15×(T−12)
Step3:太阳高度角的求解
将求解得到的赤纬角δ，时角w，当地纬度Φ代入公式(6)，求得太阳高度角的正弦值：
sinH=sinΦsinδ+cosΦcosδcosw
如北京时间9:00时太阳高度角的正弦值sinH为0.3958。
Step4:影长的求解
根据问题一已知条件，得到杆长L=3m, 将其代入公式(5),即得影长l。284+n) 365(8)284+n) 365(7)(6)
用MATLAB软件进行编程，得到影长l随时间的变化曲线图，如下所示
：
图2:影长变化曲线图
记录上图中一些特殊时刻的影长大小与其对应时间点，如下所示：
表1:特殊点影长
表
可见，北京时间9：00-15:00之间，杆影长度呈现先减小后增大的非线性变化，在12：06左右，影子长度达到最小，符合日常生活实际情况。
Step5:影长变化规律的求解
1. 影长与杆长的关系
图3:影长与杆长的关系图
由图可知，同一地点，同一日期，影长和杆长成正比关系。
2. 影长与日期的关
系
图4:影长与日期的关系图
由图可知，以1月1日为一年中的第一天，全年中影长呈现先减小后增大的非线性变化，且在夏至日左右，影长达到全年最小，冬至日左右影长达到全年最长。
3. 影长与经度的关系
图5:影长与经度的关系图
由图可知，同一纬度，在不考虑建筑物遮挡太阳光的情况下，10月22日正午时分，经度0°-41°E左右的地区太阳高度角小于0，故没有出现影长，经度41°E-51°E左右的地区由于太阳高度角极小，故影长较长，经度为45°E左右时，理论上3m高的直杆，其影长最高可达350m。经度51°E-180°E左右的地区，影长较平均，一般为20m以下。
4. 纬度与经度的关
系
图6:纬度与经度的关系图
由图可知，同一经度，在不考虑建筑物遮挡太阳光的情况下，北半球的10月22日正午时分，纬度0°-64°N左右的地区太阳高度角较大，故影长较短，
经度64°N-82°N左右的地区由于太阳高度角极小，故影长较长，理论上3m高的直杆，其影长最高可达190m。高纬度地区太阳高度角为0°，故影长为0。
5.2问题二的模型建立与求解
5.2.1模型的准备
理论影长和实际影长的定义
在一段时间内连续测得某固定直杆的影子横坐标x和纵坐标y，其中测量时间点(换算成小时为单位) 为ti,
ti=t0+i∆t(i=0, 1, 2,...m)
式中i为测量点数，m+1为测量次数，t0为初始测量时间点，∆t为测量时间间隔。
实际影长为：
′li√︀=(xi) 2+(yi) 2(9)
式中xi, yi分别为ti时刻测得的横纵坐标。
理论影长可由得影长变化模型得：
√︀Ll=sinΦsinδ+cosΦcosδcosw
式中δ为赤纬角，ψ为经度，w为时角。且有
wi=π(ti−12) 12(10)
(11)ψ=120−15(t′0−t0)
式中t0为当地时间，t′0为t0对应的北京时间。
5.2.2模型的建立
Step1:单目标优化模型的建立
为求得纬度Φ和赤纬角δ使理论影长l无限接近于实际影长l′，利用最小二乘法建立非线性优化模型：n∑︁′2min(li−li) (12)
⎧⎨sinH>0
⎩6
太阳始终在地平线以上，太阳高度角大于零，则其正弦值必大于零。且默认太阳在地平线以上的地方时在6点到18点之间。则有以使l和l′差值的平方最小为目标，以太阳高度角大于零为约束条件，建立初步的非线性优化模型。
Step2:影长比模型的建立
在上步的优化模型中，理论影长l与杆长L有关，目标函数中有杆长L、赤纬角δ和纬度Φ三个变量，不易求解，故需对该模型改进。由于理论影长l与杆长L成长比，将相邻两次li相除，可抵消L，故得到以下模型：
n′∑︁li+1li
min(−+1) 2
lili
i=0
(14)
⎧
⎨sinH>0⎩6
目标函数中,
tan Hili+1=litan Hi+1
(15)
化简后可见目标函数只与赤纬角δ和纬度Φ两个变量有关。
Step3:网格搜索
′
①在所需空间范围内，选取一个经纬位置[ψ,Φ]，逐一计算21次测量点li和li差
值的平方，并进行累加，得到点对误差ρij；
②逐一选取空间范围内的所有经纬位置，比较得到点对误差ρij最小的空间位置。
5.2.3模型的求解
附件1中已知2015年4月18日北京时间14：42至15：42每间隔3分钟测得的某固定直杆影子横纵坐标，需确定直杆所处的空间位置。
Step1:数据的处理
①真太阳时(即地方时) 的换算：初始时间点t0为14.7时，时间间隔∆t为0.05小时，其中测量次数m+1为21次。
②实际影长l′的计算：将21组xi，yi代入下式：
′
li
′
得到对应的21个li：
√︀
=(xi) 2+(yi) 2
Step2:定位空间模型的求解
′
将各测量时间点的实际影长li代入太阳影子定位空间模型，纬度δ满足
|δ|≤90
根据附件一数据可得测量时间为当地下午，地方时满足12≤T≤18, 由公式(11)可知，此条件下经度ψ应满足
19. 8≤ψ≤94. 8
在此范围内进行网格搜索。得到满足约束条件的8个区块的目标解(纬度均为北纬，经度均为东经) ：
例如在纬度[0,30]，经度[105,115]范围内以间隔为0.01搜索的结果如下图：
图7:搜索结果图
找出图像在z轴上最低的点为[36,193],即实际经纬度为[108.6°E,19.3°N]。发现仅有两个空间位置的点对误差相对其他位置极小，数量极为10−4。故认为附件1测量的固定直杆可能在2015年4月18日位于以下空间位置(纬度均为北纬，经度均为东经) ：
5.3问题三的模型建立与求解5.3.1模型的准备
以地球中心为原点，赤道平面为xy平面，指向北极的方向为z轴正向，
固定某地S，正午时刻该地所在大圆的平面为xz平面，建立直角坐标系，如下图所示：
图8:地球地标系示意图
通过查阅相关文献，得到上图中正北方向N旋转至太阳影子的方向的角度（即影子方位角）ω满足：
cosω=
cosδcoswsinφ−sinδcosφ
cosHsinω=
sinwcosδcosH
(16)(17)
同时给出cosω和sinω是为了避免出现方位角多解的情况。5.3.2模型的建立
按问题二给出的网格搜索法求解，得到每一天的最优解及相应的优化变量值，然后在得到的365个最优值中选取最优解。因为单纯采用太阳影子的长度信息进行计算，结果并不是最理想，所以考虑在目标函数中利用已知信息进行优化，即引入影子方位角。
Step1:方位角差模型的建立
由于无法确定具体的地理方位，因此如果直接采用方位角作为判别依据，就会增加参数（坐标、旋转角度）。为了优化参数，可以用方位角的差，也就是用两个时间点之间实际影子的旋转角度与理论影子的旋转角度的差方和作为判别依据。由于相邻点之间旋转角误差极小，故取相互间隔为10的点作差，得到公式如下：
min
10∑︁i=0
′′2
[(ωi+10−ωi) −(ωi+10−ωi)]
(18)
Step2:多目标优化模型的建立
在已有问题二关于影长比的单目标优化模型基础上补充日期变量，结合上述方位角差模型，建立多目标优化模型：
365∑︁20′∑︁li+1limin(−+1) 2
lili
n=0i=0365∑︁10∑︁n=0i=0
(19)
min
5.3.3模型的求解
′′2
[(ωi+10−ωi) −(ωi+10−ωi)]
(20)
在问题二中网格搜索法的基础上，增加对于日期范围的设定，确定方位角搜索的范围，搜索流程如下图所示
：
图9:网格搜索法流程图
Step1:影长比优化模型的求解
对于公式(19)，以附件2为例，使经纬度分别以1度为间隔，进行网格搜索，得出结果如下图所示
：
图10:网格搜索法流程图
Step2:方位角优化模型的求解
将所得结果带入公式(20)进行进一步优化，将经纬度的间隔分别细化为0.01度，再次进行搜索，得出最优解，如下图所示：
下面就题目所给的附件2和附件3数据，分别搜索得到以下搜索结果：
表2:附件2
日期7月25日5月26日5月27日
经度79.56°92.01°92.01°
纬度39.07°57.01°57.01°
地点
新疆维吾尔自治区喀什地区
俄罗斯南部俄罗斯南部
表3:附件3
日期2月3日
经度110.19°
纬度31.99°
地点湖北省十堰市
六、模型评价与推广
6.1模型的优点
1. 通过最小二乘法，解决了超定方程无法求出解析解的问题。2. 利用太阳高度角的限定条件和影长变化规律缩小了网格搜索范围。6.2模型的缺点
1. 不考虑太阳光的折射，计算工程量较大。
6.3模型的推广
该模型可以根据一般建筑物的影长变化，较精确地判断该建筑物的所在地以及拍摄日期和时间。
基于该模型可以求解任意建筑物一年中被其他建筑物影响光照的时间，
参考文献
[1]蔡志杰，太阳影子定位，数学建模及其应用，第4期，2015。
[2]孙吉山，洪薇．北京标准时间与各地真太阳时的换算法．上海针灸杂志，1994。[3]司守奎，孙兆亮．数学建模算法与应用．2版．北京：国防工业出版社，2015。
附录一：问题一的MATLAB 程序 d=295; %日期的天数
a=23.45*sind(360*(284+d)/365); %太阳赤纬 w=39.91;%纬度 l=3;%杆长 j=0; x=[]; y=[];
for i=9:0.1:15 %真太阳时 j=j+1;
b=sind(w)*sind(a)+cosd(w)*cosd(a)*cosd(15*(i-12));%求出太阳高度角的正弦值sinH c=b/sqrt(1-b^2);%求出太阳高度角的正切值tanH; ly=l/c;%影长 x(j,1)=ly; end
附录二：问题二的MATLAB 程序 d=108; %日期的天数
a=23.45*sind(360*(284+d)/365); x=[];
y=[1.149625826 1.182198976 1.215296955 1.249051052 1.28319534 1.317993149 1.353364049 ...
1.389387091 1.426152856 1.463399853 1.501481622 1.540231817 1.579853316 1.620144515...
1.661270613 1.703290633 1.74620591 1.790050915 1.835014272 1.880875001 1.927918447]; k=0; Z=[];
for w=0:0.1:30%控制纬度 k=k+1; t=0;
for jing=79.5:0.1:99.5%控制经度 t=t+1; z=0; j=0;
for e=14.7:0.05:15.7 %真太阳时 i=e-(120-jing)/15;%求出当地时间 j=j+1;
b=sind(w)*sind(a)+cosd(w)*cosd(a)*cosd(15*(i-12));%求出太阳高度角的正弦值sinH if b>0
c=b/sqrt(1-b^2);%求出太阳高度角的正切值tanH x(j,1)=1/c; if e>=14.75
z=((y(1,j)/y(1,j-1))-(x(j,1)/x(j-1,1)))^2+z;
Z(k,t)=z;%求累计差值 end else
Z(k,t)=10000; end end end end
[x,y]=find(Z==min(min(Z)))%寻找最小值所在的经纬度 附录三：问题三的MATLAB 程序 （1）确定方位角范围 x=[];
y=[1.247256205 1.22279459 1.198921486 1.175428964 1.152439573 1.1299174... 1.10783548 1.086254206 1.065081072 1.044446265 1.024264126 1.004640314... 0.985490908 0.966790494 0.948584735 0.930927881 0.91375175 0.897109051... 0.880973762 0.865492259 0.850504468]; Z=[]; k=0; o=0; daan=[];
for d=0:1:73%控制日期 o=o+1;
a=23.45*sind(360*(284+d)/365); for w=0:1:90%控制日期 k=k+1; t=0;
for jing=19.8:1:94.8%控制经度 j=0; t=t+1; z=0;
for e=12.68:0.05:13.68 i=e-(120-jing)/15; j=j+1;
b=sind(w)*sind(a)+cosd(w)*cosd(a)*cosd(15*(i-12));%求出太阳高度角的正弦值sinH
if b>0
c=b/sqrt(1-b^2);%求出太阳高度角的正切值tanH x(j,1)=1/c; if e>=12.73
z=((y(1,j)/y(1,j-1))-(x(j,1)/x(j-1,1)))^2+z; Z(k,t)=z;%求累计差值 end else
Z(k,t)=10000;
end end end
[u,v]=min(min(Z));%寻找当天搜索范围内累计差值最小的点 daan(o,1)=u;%输出每日的最小累计差值 end
（2）方位角差方和最小 x=[];
y=[-82.02711405 -81.10849053 -80.16448976 -79.19768875 -78.20352191 -77.19123876... -76.14438496 -75.07333958 -73.9762777 -72.84887252 -71.68887638 -70.4973275 ... -69.27179577 -68.0170072 -66.72685728 -65.39508527 -64.03225482 -62.63232339... -61.18607201 -59.70496474 -58.18242197]; Z=[]; k=0; o=0; daan=[];
for d=120:1:180%控制日期 o=o+1;
a=23.45*sind(360*(284+d)/365); for w=35:0.1:45%控制日期 k=k+1; t=0;
for jing=80:0.1:85%控制经度 j=0; t=t+1; z=0;
for e=12.68:0.05:13.68 i=e-(120-jing)/15; j=j+1;
b=sind(w)*sind(a)+cosd(w)*cosd(a)*cosd(15*(i-12));%求出太阳高度角的正弦值sinH
if b>0
c=sqrt(1-b^2);
C=sind(a)*cosd(15*(i-12))/c;%求出影长变化的正切值 da=asind(C);%求出角度 x(j,1)=da; if e>=13.18
z=((y(1,j)-y(1,j-1))^2-(x(j,1)-x(j-1,1))^2)^2+z; Z(k,t)=z;%求出方差和 end else
Z(k,t)=10000;
end end end
[u,v]=min(min(Z));%寻找当日搜索范围内最小的方差和 daan(o,1)=u;输出每日的最小方差和 end