15年全国大学生数学建模竞赛省级一等奖

太阳影子定位

摘要：本文通过利用赤纬角、时角、太阳高度角、时间差等地理概念，建立利用影子的变化情况，确定经纬度、日期的数学模型。

问题一、本问题中日期、时间点、经纬度、杆高都是影响影长变化的因素。利用太阳高度角公式、时角计算公式以及太阳赤纬角计算公式，即可建立日期、时间点、经纬度、杆高与影长之间关系的数学模型。并将所给数据代入上述模型中，画出了天安门广场直杆影长的变化情况（见图1）。

问题二、此问题要求根据影子顶点坐标、日期、北京时间来确定经纬度。通过附件1中时间和影子顶点坐标数据拟合出影子长度随时间变化的二次曲线，求出北京时间与地方时间的时间差，进而利用它们之间关系求出经度值。通过借助Analemmatic 日晷模型，根据日晷产生的椭圆轨迹在垂直和水平方向的分量比值与物体投影在垂直和水平方向的分量之间关系，建立模型求解纬度。将附件1的坐标数据、时间、日期代入上述模型，即可求解出经纬度坐标（见表1）。

问题三、此问题要求根据影子顶点坐标、给定北京时间来求解测量地经纬度和日期。利用问题二中二次曲线拟合的方法可确定经度值，借助太阳高度角、时角、赤纬角、太阳方位角、纬度之间的关系以及它们各自与影子顶点坐标的关系，建立方程组求出纬度值，太阳高度角和赤纬角，再利用赤纬角进而计算出日期。将附件2、3数据代入上述模型，计算出纬度及日期。

问题四、已知日期求解拍摄地点时，通过Matlab 获取附件4中视频每一帧图片的信息并将图片转化为灰度值矩阵。通过比例可算出相应时刻对应的直杆影长，使用问题二中的方法即可求出拍摄地点的经纬度。如果日期未知，可使用问题三中的方法求出地点和日期。

关键词：太阳赤纬角；时角；太阳高度角；Analemmatic 日晷；北京时间；当地时间；太阳方位角；测量地的经纬度；

一、 问题的重述

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。

1. 建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒, 东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。

2. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。

3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。

4．附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。

如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？

二、 模型假设及符号说明

1、 模型假设

(1) 不考虑各附件中坐标的测量误差。

(2) 视频处理时，不受视频中风速、噪音等外界因素的影响。

(3) 不考虑视频拍摄角度对测量结果的影响。

2、 符号说明

N 表示某年月日距离当年1月1日的日数

T 表示具体时刻值

h 表示观测地的太阳高度角

ϕ 表示观测地的地理纬度（北纬为正，南纬为负）

δ 表示赤纬角（北纬为正，南纬为负）

t 表示地方时，即时角

H 表示测量杆的高度

l 表示影子长度

l b 表示北京时间时区的中央子午线的经度值

l d 表示投影物体所在位置的经度值

t c 表示当地时间和北京时间的时间差

A 表示太阳方位角

三、 问题的分析

问题一、本问题中测量日期、测量时间点、测量地的经纬度、测量杆的高度都是影响影子长度变化的因素。其中测量日期会影响到赤纬角的变化；测量时间点会影响到时角的变化；而赤纬角、时角、测量地的经纬度会直接影响太阳高度角的变化，间接影响到影长；杆高也会直接影响到影长的变化。根据此关系，利用一般时间的太阳高度角公式、时角计算公式以及太阳赤纬计算公式，即可建立日期、时间点、经纬度、杆高与影长之间关系的数学模型。

问题二、此问题要求根据影子顶点坐标、日期、北京时间来确定经纬度。通过附件1中北京时间和影子顶点坐标数据，利用Matlab 拟合出以北京时间为横轴，影子长度为纵轴的二次曲线，求出二次曲线最低点即影长最短时对应的时间值，此时间为当地时间正午12时对应的北京时间，求出北京时间和当地时间的偏差值，进而可求出经度值。求解纬度时，可以考虑到利用“Analemmatic 日晷”模型进行求解。“Analemmatic 日晷”是一种灵活并且可以移动的日晷，是人类古代利用日影测得时刻的一种计时仪器，其原理就是利用太阳的投影方向来测定并划分时刻，利用日晷产生的椭圆轨迹在垂直和水平方向的分量的比值以及物体投影到垂直和水平方向的分量之间关系即可求解纬度。

问题三、此问题要求根据影子顶点坐标、北京时间确定经纬度和日期。利用问题二中二次曲线拟合的方法可确定经度值，借助太阳高度角、时角、赤纬角、太阳方位角、纬度之间的关系以及它们各自与影子顶点坐标的关系，建立方程组求出纬度值、太阳高度角和赤纬角，再利用赤纬角进而计算出日期。

问题四、已知日期求解拍摄地点时，通过Matlab 获取附件4中视频每一帧的图片信息并将图片转化为灰度值矩阵。通过比例算出不同时刻对应的直杆影长，使用问题二中的求经度的方法即可求出拍摄地点的经度值，利用问题一中的方法反推可求出纬度值。如果日期未知，可使用问题三中的方法求出地点和日期。

四、 模型的建立及求解

4.1问题一模型的建立及求解

在此问题中，影子长度的变化受日期、时间、观测地的地理纬度和测量杆高度这四个变量的影响，利用各变量与影子长度之间的关系建立模型，求解出影长与各变量之间的函数关系。

4.1.1模型的建立

太阳赤纬角是指地球赤道平面与太阳和地球中心的连线之间的夹角。赤纬角以年为周期，变化范围为+23︒26'至-23︒26'。春分和秋分时期，赤纬角为0︒，此时太阳直射赤道；夏至时期，赤纬角为+23︒26'此时赤纬角达到最大值，太阳直射北回归线；冬至时期，赤纬角为-23︒26'此时赤纬角达到最小值，太阳直射南回归线。赤纬角δ计算公式：

sin δ=0.39795cos[0.98563(N -173)] （1）

其中，N 表示日数，从每年1月1日开始计算； δ表示赤纬角（北纬为正，南纬为负）。

时角是指直杆在地球赤道平面上的投影与当地时间12点时、地中心连线在赤道平面上的投影之间的夹角。其中正午12点的时角为0度，每隔一小时增加15度。时角计算公式：

t =(T -12) ⨯15︒ （2）

其中，T 表示具体时刻值； t 表示地方时，即（时角）。

太阳高度角是指在地球上某地点，太阳光的入射方向和地平面之间的夹角。太阳高度角计算公式：

sinh =sin ϕsin δ+cos ϕcos δcos t （3）

利用太阳高度角确定影子长度L ：

L =H （4） tanh

其中，h 表示观测地的太阳高度角；H 表示测量杆的高度。

将（1）（2）(3)式代入（4）式中，结合观测地所在的地理纬度，即可得到影长与赤纬角、时角、太阳高度角的关系式。

4.1.2模型的求解

从题目所给数据即可获得1月1日至10月22日，N 值为295天；具体时

ππ间T 从9:00至15:00，则时角t 变化范围为-至；观测地（北京）的地理纬44

度ϕ为+39︒54'26"；测量杆高度H 为3m; 将以上数据带入模型一中，即可得到太阳影子长度随时间的变化曲线如图1所示：

图1 太阳影子长度随时间的变化曲线

由图1可以观察到，9点到12点时，随着时间的增加，影长不断减少；在当地时间正午12点时影子长度达到最短；12点到15点时，随着时间的增加，影长不断增加；可得出影长随时间的变化整体呈开口向上的二次曲线。

4.2问题二模型的建立与求解

在本问题中，我们借助“Analemmatic日晷”[1]模型用直杆的影子坐标来计算直杆所处的地点经纬度，计算利用北京时间和当地时间的时间差，进而确定测量地的经度，根据日晷的设计原理来计算，运用几何知识构成比例来求解其纬度值。

“Analemmatic日晷”是一种灵活并且可以移动的日晷，这种日晷可以根据在不同的月份站立在不同位置的行人来作为指针，它的小时指针是赤道日晷在水平面的投影，因此它的轨迹是一个椭圆，并且日晷上面的时刻表示了阴影的方向。

4.2.1 建立Analemmatic 日晷模型

图2 Analemmatic日晷模型

在图2中日晷上的时间点A 指示当地的时间为下午三点，B 点为实际阴影点，O 为日贵的中心，E 为投影物的脚点位置。

4.2.2 建立求解经度模型

由于当地子午线与北京时间地子午线存在着偏差，当地的日晷时间和北京时间之间存在着一个偏差，由于数据中给出的是北京时间而不是当地的时间，因此在计算经度的时候必须考虑它们之间的差值，因此要求当地的经度值，必须求出当地时间和北京时间的偏差值，进而计算出当地经度值，其计算过程如下： 由于已知不同时刻直杆关于x 、y 坐标的投影值，利用两点间距离公式求出影子端点坐标与原点的距离，即某一时刻的影长l ：

l = （5）

由图1画出的曲线可以看出“时间-影长”函数曲线大致为二次曲线，并且曲线最低点所对应的时间为当地时间正午12:00，则我们可以进行拟合“时间-影长”二次函数来求取当地时间，待拟合的二次曲线为：

y =ax 2+bx +c （6）

带入附件1中数据，解出系数a 、b 、c 的值，利用拟合曲线计算出影长最短y min 所对应的x 点的坐标值，即影长最短时对应的北京时间。由于当地时间为正午12点时，影长也会达到最短，可得到时间差t c 。利用北京时间所在时区的中央子午线l b 和时间差t c ，可求出投影物体所在的经度值l d ：

l d =l b -0.25⨯t c （7）

其中，l b 表示北京时间时区的中央子午线；l d 表示投影物体所在位置的经度值；t c 表示当地时间和北京时间的时间差。

4.2.3 建立求解纬度模型

Analemmatic 日晷在当地纬度产生的椭圆轨迹的行成中水平方向可以定义为：

x l =sin t (8)

其中，时间角度t =(T -12) ⨯15︒，T 为（1）中计算出的当地时间。垂直方向可以定义为：

y l =sin ϕcos t (9)

其中，ϕ为当地的纬度值，则比较公式（8）（9）可以看出，水平方向的坐标计算只需要时间角度值，而垂直方向的坐标计算需要当地纬度值，投影物体的脚点相对于Analemmatic 日晷椭圆中心的位置可以定义为：

y o =tan δcos ϕ (10)

其中，δ为太阳直射高度角， δ在一年中的变化范围大致为-23︒26'~23︒26'。当地纬度值可以通过下面的公式进行计算：

y l -y o

x l =y x （11）

其中，y 为影长在y 轴的投影，x 为影长在x 轴的投影。

4.2.4 经度的求解

由附件1中各个时刻的影子坐标可求得影子长度，用Matlab 求出横坐标x 为时间，纵坐标y 为影子长度的二次曲线为：

y =0.1489⨯x 2-3.752⨯x +24.13 (12)

已知数据拟合曲线如图3所示。

图3已知数据拟合曲线图

在图3中，误差平方和=1.649⨯10-5，标准差=0.0009571，其值都较小，可以认为拟合的二次曲线拟合度较高。

由于得出了已知数据拟合的二次曲线（抛物线）图，并且一直其方程则我们可以画出一天时间内，物体影子长度随时间变化图，我们画出早上6:00到下午18：00的影子长度变化如图4所示。

图4 影长关于时刻的拟合曲线图

在图4中，拟合的二次曲线类似于开口向上的抛物线，则它在x 的取值范围内有最小值，其最小值所对应的x 轴的值经计算为约12.599，即对应的时间为12:36,说明北京时间12：36分时对应的当地时间为12：00，则北京时间与当地时间的偏差值t c =36分钟。

已知北京所处位置位于东八区，且东八区的中央子午线为120︒，即l b =120︒。则由（7）式可以计算出l d =111︒

4.2.5 纬度的求解

由（6）式可知附件1中21个影子坐标点对应的当地时间t 1 t 21；当天日期为2015年4月18日，可以算出太阳直射高度角δ=7.4462︒。利用Matlab 根据公式（11）可求得附件1中21个时刻对应的纬度值如下表所示：

表1 各时刻求得纬度值

由表1可以看出，纬度值范围为25.6131︒ 31.09098︒。由4.2.4解得的经纬度为111︒，可大致确定的地点为：广西省，湖南省，湖北省。

4.3问题三模型的建立及求解

此问题中给出多个时间点及各时间点影子的坐标值，要求根据时刻值和坐标值建立数学模型，解出直杆所处的经纬度及测量日期。

4.3.1建立模型求解经度值

在此问中求解投影物体所在经度值的做法与4.2.2中求解经度值的做法相同。由于已知不同时刻直杆关于x 、y 坐标的投影值，利用两点间距离公式求出影子端点坐标与原点的距离，即某一时刻的影长l 。以时间为自变量，影长为因变量建立拟合函数关系，拟合曲线假设为y =ax 2+bx +c 。求出曲线最低点对应x 的值，即为当地时间为正午12点所对应的北京时间，进而算出时间差t c 。

通过北京时间所在时区的中央子午线l b 和时间差t c ，利用如下公式，即可求出投影物体所在的经度值l d ：

l d =l b -0.25⨯t c （13）

4.3.2建立模型求解纬度值及赤纬角

利用公式（3）建立太阳高度角h 与观测地的地理纬度ϕ、赤纬角δ、时角t 之间的关系模型，即：

sinh =sin ϕsin δ+cos ϕcos δcos t （14）

利用问题二中根据日晷求解观测地纬度的方法，建立不同时刻影子坐标点与观测地的地理纬度ϕ、赤纬角δ、时角t 之间的关系模型，即：

y sin ϕcos t -tan δcos ϕ= （15） x sin t

太阳方位角[2]一般是以目标物的北方向起始方向，以太阳光的入射方向为终止方向，按顺时针方向所测量的角度。对于中国区域而言，早上太阳从东方升起，傍晚从西方落下，两者方位角分别为90︒、270︒左右（由于季节变化，方位角也会随之变化）；正中午的太阳方位角在180︒（正南方）。

利用太阳方位角作为中间关系，建立不同时刻影子坐标点、赤纬角δ、时角t 、太阳高度角h 与观测地的太阳方位角A之间的关系模型，即：

y =t a n A [3] （16） x

c o δs s t i n A = s i n （17） c o s h

将公式（14）（15）（16）（17）进行整合，建立一个时角t 、影子坐标点与观测地的太阳高度角h 、观测地的地理纬度ϕ、赤纬角δ、太阳方位角A的关系模型。已知时角t 和影子坐标点，利用所给的方程组即可求出观测地的太阳高度角h 、观测地的纬度值ϕ、赤纬角δ这三个变量。

⎧y sin ϕcos t -tan δcos ϕ⎪x =sin t ⎪sinh =sin ϕsin δ+cos ϕcos δcos t ⎪⎪ ⎨y ⎪x =tan A

⎪⎪sin A =cos δsin t

⎪cosh ⎩ （18）

4.3.3建立模型求解测量日期

根据所求的赤纬角δ，利用公式（1）反推，即可求得日数N ：

sin δ=0.39795cos[0.98563(N -173)] （19）

通过日数N 从而确定测量日期：具体测量日期=1月1日+日数N 。

4.3.4观测地的经度求解

由附件2中各个时刻的影子坐标可求得影子长度，用Matlab 拟合出影子长

度随时间变化的二次曲线：

y =0.09814⨯x 2-2.356⨯x +14.76 （20）

已知数据曲线拟合结果如图5所示。

图5 已知数据曲线拟合图

在图5中，误差平方和= 4.403⨯10-7，标准差=0.0001564，其值都较小，则可以认为拟合的二次曲线拟合度较高。

求解的二次曲线类似于开口向上的抛物线，则它在x 的取值范围内由最小值，其最小值所对应的x 轴的值经计算约为15.2079，即对应的时间为15:12,说明北京时间15:12分时对应的当地时间为12：00。

则可以求出每个北京时间对应的当地时间，也可以知道当地时间和北京时间的偏差值t c =192分钟。

同理可得到附件3拟合的二次曲线方程为：

y =0.2964x 2-7.551x +51.56 （21）

已知数据曲线拟合结果如图6所示。

图6 已知数据曲线拟合图

在图6中，误差平方和= 8.106⨯10-6，标准差=0.0006711，其值都较小，则可以认为拟合的二次曲线拟合度较高。

求解的二次曲线类似于开口向上的抛物线，则它在x 的取值范围内由最小值，其最小值所对应的x 轴的值经计算约为12.7379，即对应的时间为12:44，说明北京时间12:44分时对应的当地时间为12：00。当地时间与北京时间的偏差值t c =44分钟。

已知北京所处位置位于东八区，且东八区的中央子午线为120︒，即l b =120︒。则由（13）式可以计算出附件2中l d =72︒。同理可求得附件3中的l d =109︒

4.4问题四模型的建立及求解

4.4.1已知日期求解拍摄地点

将附件4中给出的阳光下直杆投影随时间变化的视频导入Matlab 中使其视频将每一帧都以图片的形式输出，将输出的图片导入Matlab 中，由于输出的图片为彩色图片，我们首先要将彩色图片转化为黑白图片，得出黑白图片的灰度值矩阵。

已知直杆的长度为两米，我们找出直杆在图片的灰度矩阵中的高度，由于每一个时刻对应一张图片，因此我们找出每张图片中直杆影子在灰度矩阵中的高度，由比例关系可得出每个时刻所对应的影子长度，即：

X 1X 2 （22） =x 1x 2

其中，X 1为直杆真实高度，x 1为直杆在灰度矩阵中高度，X 2为直杆影长真实高度，x 2为直杆影长在灰度矩阵中高度

Step1计算经度

根据公式（22）所计算的影子长度，利用4.2.2中计算经度的方法，用Matlab 拟合出时间影响影子长度变化的二次曲线，求出曲线最低点所对应的时间值即为当地正午12时所对应的北京时间，进一步可求出两地时间差，从而确定测量地的经度值。

Step2计算纬度

根据公式（22）所计算的影子长度l 和直杆高度L 即可确定太阳高度角h

L tanh = （23） l

在附件4所给的视频中可知道视频的拍摄日期，由问题一中的公式（1）可计算出这一天的赤纬角，即太阳直射角度δ。Step1中计算出了当地时间，根据

时角公式计算出每时刻所对应的时角t 。利用太阳高度角公式，已知太阳高度角h 、赤纬角δ、时角t 即可反推出当地纬度值ϕ：

sinh =sin ϕsin δ+cos ϕcos δcos t （24）

4.4.2求解拍摄地点和日期

Step1 计算测量地的经度值

由于经度的计算与日期无关系，则可以根据4.4.1 Step1中计算测量地经度的步骤和方法进行计算。

Step2 计算测量地的纬度值

由于计算纬度的时候没有告诉日期，可以利用4.3.2 中计算测量地纬度的方法进行计算。

⎧y sin ϕcos t -tan δcos ϕ⎪x =sin t ⎪sinh =sin ϕsin δ+cos ϕcos δcos t ⎪⎪⎨y ⎪x =tan A

⎪⎪sin A =cos δsin t

⎪cosh ⎩ (25)

将其整理为

sin ϕcos t -tan δcos ϕ⎧tan A =⎪sin t ⎪⎨sinh =sin ϕsin δ+cos ϕcos δcos t (26)

⎪cos δsin t ⎪sin A =cosh ⎩

Step3 计算测量日期

太阳高度角h 、时角t 已知，代入上式即可求出太阳方位角A 、测量地纬度ϕ、赤纬角δ，利用赤纬角公式即可求出日数N ：

sin δ=0.39795cos[0.98563(N -173)] （27）

通过日数N 从而确定具体测量日期：具体测量日期=1月1日+日数N 。

五、模型的综合评价

此论文中能够根据太阳高度角、赤纬角、太阳方位角和时角等地理概念建立影长、时间、日期、经纬度的数学模型，将地理概念与影响影长的因素密切相连，使模型的建立更加科学、客观。

在问题一中能够充分考虑影响影长的变化因素，根据测量日期、测量时间点、

测量地经纬度以及杆高建立影响影子长度的数学模型；而在问题二根据日晷的设计原理以及杆与影之间的关系，建立模型求解经纬度，使得经纬度的计算更加方便、快捷；问题三和问题四的求解运用方程组简单快捷地解出其包含的未知数即为问题所要求解的值；但是在模型的建立与求解中考虑的因素不足，并且所用到的公式大多数为三角函数，使得方程组的求解难度较大，数值不便统计。 参考文献：

[1] 武琳，基于太阳阴影轨迹的经纬度估计技术研究，天津大学硕士学位论文，第27页至第28页，2010年12月。

[2] 贺晓雷，于贺军，李建英，太阳方位角的公式求解及其应用，太阳能学报，第29卷第1期，第69页，2008年1月。

[3] 郑鹏飞，林大钧，刘小羊，吴志庭，基于影子轨迹线反求采光效果的技术研究，华东理工大学学报（自然科学版），第36卷第3期，第460页，2010年6月。