建模2015A题优秀论文

基于L-M优化算法的太阳影子定位模型

摘要

太阳影子定位技术是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和

日期的一种方法。

本文旨在通过数学建模的方法对太阳影子长度和角度的变化与拍摄地点和日期的关系

进行分析，从而研究该定位技术。

针对问题一，我们利用向量知识和几何关系推导得出影子长度与标杆高度h、拍摄

地纬度、太阳直射点纬度和拍摄地与太阳直射点经度差四个参数之间的函数关

系。而由拍摄日期决定，可由该地地方时表示；同时由于拍摄时刻一般记录的都是北京时间，所以进行地方时和北京时间的转换。分别推导得出随日期变化的公式和由北京时间表示的公式，带入原函数关系，建立影长与杆长、拍摄点纬度、拍摄日期和拍摄北京时刻之间关系的数学模型。

针对问题二，由于所给数据中影长随北京时间呈上升趋势，根据北京时间与地方时

的转换关系可以确定拍摄地大致的经度范围为(E79.50,E154.50)。然后在这个纬度范围内利用Levenberg-Marquart算法[8] ，以影长偏差为目标函数，对参数进行基于最小二乘法的非线性拟合，得出可能地点的地理坐标为：(E108.630,N19.090)。为提高数据利用率，经坐标转换后，以上一步求得的经纬度各增减200作为经纬度搜索范围，将x，y值偏差作为目标函数再次进行拟合，得出拍摄地的地理坐标为：(E109.830,N18.150)。 针对问题三，由于本问与问题二相比差别仅在于拍摄日期未知，且L-M算法能有效

处理冗余参数问题，所以采用与问题二基本相同的处理方法。先根据影长随时间变化趋

势确定附件2地点的经度范围(19.750,94.750)，附件3的经度范围(102.750,1800) (1800,167.250)。然后经过两次参数拟合，得出附件2地点的坐标：(E79.730,N39.670)，拍摄日期为7月23日；附件3地点的坐标：(E111.140,N39.390)，拍摄日期为11月8日。 针对问题四，先对视频每隔2分钟截取一帧，设图像中影长l与拍摄正面的夹角为，

通过像素坐标确定l在拍摄正面的投影lcos与杆高h的比值，进而得出影长。当拍摄日期已知时，问题转化成问题二中以影长偏差为目标函数的非线性拟合模型，区别在于增加了未知量，利用L-M算法求解拍摄地坐标，结果为(E111.140,N39.390)；当拍摄日期未知时，问题转化成问题三中以影长偏差为目标函数的非线性拟合模型，区别也是在于增加了未知量，利用L-M算法求解拍摄地坐标为(E110.570,N39.600)，日期为7月10日。

关键词：太阳影子定位，L-M算法，非线性拟合

一、问题重述与分析

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。在本文中，我们主要对以下问题进行讨论：

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

分析：首先，以立杆点为原点建立平面直角坐标系，利用向量知识和几何关系进行分析发现影长可由标杆高度、拍摄地纬度、太阳直射点纬度和拍摄地与太阳直射点经度差四个量直接表示；而题目所给的参数为拍摄点地点地理坐标、日期和时刻，所以应当寻找题目参数和影长表示参数之间的联系。经进一步的分析，我们发现太阳直射点纬度实际上是由一年中的日期决定，而两地经度差可以由该地时刻表示，所以通过公式的带入可以建立影长与题目所给参数之间的函数关系。

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

分析：在本问题中，已知量为日期和某地从北京时间14:42到15:42不同时刻影子 顶点坐标，要求确定拍摄地的经纬度。值得注意的是，题目中并没有给出杆的长度和坐标系的方向。所以设题目所给坐标系与我们在问题一中建立的坐标系的夹角为，此时待求的未知量为拍摄点经度、纬度、杆长和坐标偏转角度。为求解未知参数，一个方式是将所给数据代入问题一已经得出的影长与参数关系式中进行解方程。在理论上只需要四个数据就可以解含有四个未知量的方程，这就出现了对数据的利用率不高的情况，会降低结果的精确度；因此我们采用非线性拟合的方式求解参数，以理论计算值和实际数据之间的偏差作为目标函数，将问题转化成一个基于最小二乘法的优化模型，利用L-M算法求解。

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

分析：本问与问题二相比差别仅在于拍摄日期未知，相当于多了一个未知参数。由于L-M算法能有效处理冗余参数问题，所以可以采用与问题二基本相同的处理方法。值

得注意的是，因为附件2和附件3所给数据的x坐标正负相反，所以我们认为这两组数据是在不同的拍摄地点。

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

如果拍摄日期未知，要求根据视频确定出拍摄地点与日期。

分析：在本问题中，首先要在视频中提取图像，由于整个视频里没有出现影长的突变，所以我们每隔2分钟均匀提取。然后根据杆和影子的像素坐标确定两者长度的比值，可以得出影长的数据，则当日期已知时，问题蜕化成问题二的模型；日期未知时，问题蜕化成问题三的模型。值得注意的是，影子不一定与拍摄正面平行，所以我们考虑设两者的夹角为，在求解时，未知参数增加了一个，但依然可以用L-M算法求解。

二、模型假设

1、忽略海拔高度的影响。因为即使是海拔较高，但是与太阳和地球之间的距离相比，海拔的影响就非常小了，可以忽略。

2、假设一天中太阳直射纬度不变。一天中地球公转不到1o，所以假设一天内太阳直射纬度不变。

3、假设相对地球而言太阳是一个面光源，照射到地球的太阳光是一组平行光线。

4、假设地球上某地的水平地面是地球球面上过该地的切面。

5、忽略大气折射的影响。

三、符号说明 



x

y拍摄地纬度 拍摄地经度 太阳光线所对应与拍摄地水平地面的法向量夹角 影子顶点横坐标 影子顶点纵坐标

太阳直射点纬度

拍摄地与太阳直射点地的经度差

 

T 北京时间

拍摄点地方时

黄经度

影子长度

杆的高度

拍摄日期距春分日的天数

视频中影子与拍摄正面的夹角

拍摄日期距元旦的天数 t  l h n



m

四、模型的建立与求解

4.1问题一

4.1.1模型的建立

 建立两个坐标系

设地球上某地的水平地面是地球球面上过该地的切面，向量AE是与过A处的经线相切且方向向北的单位向量(如图2)，AK也是A处地平面内方向向北的单位向量(如图

1)；设向量AE是与过A处的纬线相切且方向向东的单位向量(如图2)，AE也是A处地平面内方向向东的单位向量(如图1)；向量AK与AE确定A处地平面．以AE与AK所在直线为坐标轴在A处地平面上建立平面直角坐标系Axy，在坐标系Axy下设日影F(x,y)。

图1 地平面坐标系示意图[1] 图2 拍摄点与太阳直射点经纬度示意图[1]

如图2，设地球O的球面上过A地的经线与赤道交于D点，设某日某时太阳直射地球上B地，过B的经线与赤道交于C点。设地球半径为R，AOD，BOC，DOC,AOB,则为观测地A处的纬度数(一90≤≤90),若A地在北半球，oo

则>0，若A在南半球，则

=(12-t)×15 o(1)

其中，O≤t≤24。

又是太阳光线所对应的向量OB与A地水平地面的法向量(或AH)的夹角，图3是过A，B两地的大圆，因为HF//BO(太阳光线是平行的)，所以AOBAHF, 90o-即A地t时刻太阳高度角。考虑太阳与地球的实际，对A地白昼时刻有

||90o。

如图2，以O为原点，以OD所在直线为x轴，地轴ON所在直线为z轴建立空间直

角坐标系Oxyz，则AE(0,1,0)，AK(sin,0,cos)，A(Rcos,0,Rsin)，

B(Rcoscos,Rcossin,Rsin)，coscosOA,OBcoscoscossinsin。

同时，为统一计算，规定东经为正值，西经度数为负值；北纬度数为正值，南纬度数为负值。

 日影坐标公式推导

图3 拍摄点与太阳直射点纵切面示意图

AHh如图3，在RtAHF中，HF。设HF与AE所成角为，则 coscos

coscosHF,AEcosBO,AEcossin

如图1，对HF在AE上的正射影AJ，有AJHFcoscosh，即F点在平面Axy上 cos

cossinh。设HF与AK所成角为，则 的横坐标xAJcoscoscossinsin

coscosHF,AKcosBO,AKsincoscoscossin，

如图1，对HF在AK上的正射影AG，有AGHFcos

上的纵坐标yAGcosh，即F点在平面Axycossincoscoscossinh，故在平面Axy上， coscoscossinsin

h (2) cossinsincoscoscossinFh,coscoscossinsincoscoscossinsin

 日期与太阳直射点纬度的关系

经查阅资料可知，太阳直射点的纬度可以利用黄经推导[5]，计算公式为：

arcsin(0.397775sin)

其中，为太阳直射点纬度，为黄经度。 (3)

值得注意的是，地球公转的周期是一个回归年(365. 2422日).现行公历的历年是历日的整数倍,它和回归年并不精确相等;另外,由于复杂的历史演变过程,以及一些人为原因,上半年和下半年、冬半年和夏半年以各季、各月之间的天数也并不完全相同,因而日期和黄经并不完全相对应。因此,用日期取代黄经来推算直射点纬度时,需按季节分段时行计算,以确保推算的相对精确[2]。

首先，对于夏半年，从春分日(3月21日前后)到秋分日(9月23日前后)总计186天,在此时段内视太阳在黄道上运转180°.设从春分日开始,视太阳运行了n天,则n天运行了经度:

1800

n 186

代入(3)式,得: 1800

arcsin[0.39775sin(n)] 186

其中，n[0,186]。 (4)

对于冬半年，自秋分日(9月23日前后)至冬至日(12月22日前后)总计90天,在此时段

内视太阳的黄道上运转90°.设从春分日开始,视太阳运行了n天,则n天运行了经度:

900

180(n186) 900

代入(3)式,得:

arcsin[0.39775sin(n186)0]

其中，n[186,276]。 (5)

自冬至日至次年春分日总计89天,在此时段内视太阳在黄道上运转90°.设从春分日开始,视太阳行了n天,则n天运行了经度:

900

270(n276) 890 代入(3)式,得:

900

arcsin{0.39775cos[(n276)]} 89

其中，n[276,365]。

 地方时与北京时的时差换算 (6)

由（1）式可知拍摄地A地与太阳直射点地B的经度差由本地时t决定。而在实际 生活或者是在拍摄视频时，我们所记录下的时间都是北京时间，它与本地实际时间有一定的时差，这个时差在计算中不能忽略。经分析可得，设北京时间为T，A处的经度数为，则

tT(1200)4min 60(7)

 影子长度数学模型

由式（4），（5），（6）可知，在不同的日期范围内，太阳直射点纬度可由日期表示；由式（1）可知拍摄地与太阳直射点经度差可由地方时t表示，而由式（7）可知，地方时t可由拍摄地经度和北京时间表示。

以日期n[0,186]为例，将式（4）、（1）和（7）带入式（2）可得拍摄点坐标：

x

ycossinh coscoscossinsinsincoscoscossinh coscoscossinsin

o=(12-t)×15 （8）

tT(1200)4min 60

1800

arcsin[0.39775sin(n)] 186

此外，当日期n[186,276]时，arcsin[0.39775sin(n186)0]；当日期n[276,365]时，

900

arcsin{0.39775cos[(n276)]}。 89

由上述公式可知，影长实质上是关于杆长、拍摄点经度、纬度、拍摄日期和拍摄时北京时间五个量的一个函数。但是，若拍摄记录的时间即为本地时间，则不需要北京时间与本地时式（7）的转换，影长实质上是关于杆长、拍摄点纬度、拍摄日期和拍摄时本地时刻四个量的一个函数。

4.1.1模型的求解

 影长关于参数的变化规律

1、影长随经度变化的规律

为了研究影长关于经度的变化规律，需要将其它参数取为定值。取日期为9月3日 北纬270处，长为1米的杆，在格林尼治6时、12时和18时绘制影长随经度变化的曲线如下所示：

9月3日北纬27度1米杆随经度变化影长变化

影长/m

经度/度

图4 影长随经度变化曲线图

以单条曲线为例分析：由上图可知，在格林尼治12时，0度经线上影子最短；向南北方向，影长随纬度的增加而减小。在实际情况中，格林尼治12时0度经线上为正午，向南北方向经线分别在上午和下午，所以0度经线影长最短，两侧逐渐边长。这说明图像是符合实际的。

同时从三条曲线观察，发现随着时间的推移，曲线可以认为逐渐向右移动。在0度经线处，从6时至18时，影长经历由长变短再逐渐变长的变化，符合实际。

2、影长随纬度变化的规律

取9月3日，杆长1米，当地时间分别为9时、12时和14时作图；因为所取的时间为本地时间，所以此时影长为关于拍摄点纬度、日期和拍摄时刻三个量的函数，则不需要在给定拍摄地经度值。绘制影长随纬度变化的曲线如下所示：

影长/m纬度/度

图5 影长随纬度变化曲线图

以单条曲线分析，影长曲线大概关于赤道对称。而在实际情况中，9月3日在秋分日左右，秋分日时太阳直射赤道，赤道影长最短，向两侧纬度影长逐渐增加。这说明图像是符合实际的。

同时从三条曲线来看，由图可得本地时12时的影长最短，9时的影长最长，也是符合实际的。

3、影长随日期变化的规律

取北纬27o处，杆长1米，当地时间分别为9时、12时和14时作图；因为所取的时间为本地时间，所以此时影长是关于拍摄点纬度、日期和拍摄时刻三个量的函数，则不需要在给定拍摄地经度值。绘制影长随日期（用距春分日天数表示）变化的曲线如下所

示：

影长/m

距春分日天数/d

图6 影长随日期变化曲线图

以单条曲线为例分析：北纬27在当地时间14时，由图可知影长经历先缩短再变长的变化。在距春分日100天左右（即夏至日附近）达到最短，距春分日280天左右（即冬至日附近）达到最长。在实际情况中，北纬27o处在北回归线附近，所以应当在夏至日左右影子达到最短，冬至日左右影子达到最长；这说明图像是符合实际的。

同时从三条曲线分析，由图可得本地时12时的影长最短，9时的影长最长，符合实际。

 天安门广场影长随时间变化

题目要求绘制天安门2015年10月22日北京时间9:00-15:00之间3米高的直杆的太阳影子长度的变化曲线。需要注意的是，天安门的本地时间并不是北京时间，所以需要利用式（7）进行本地时间和北京时间的转换。根据式（8），将天安门广场的经纬度、拍摄日期带入，做出影长随北京时间变化的曲线图如下所示：

10月22日天安门影长变化

87

影长/m

65

43时间/h

图7 天安门影长随时间变化曲线

4.2问题二 4.2.1模型的建立

 坐标系的转化

在本问题中，数据给出了太阳影子顶点的坐标值，但是题目所给数据只说明了以杆子底部为原点，并未说明x轴和y轴的指向，所以可以引入一个变量来表示题目所给数据的坐标系y轴与我们在问题一中建立的以杆子底部为原点，正北方向为y轴的坐标系Axy中y轴的夹角。如下图所示：

图8 坐标系转换示意图

由图可知，两个坐标系坐标的关系如下：

'xcosysin

y'ycosxsin （9）

 拍摄地经度范围的估计

由附件所给的数据做出影长随时间变化的图像如下：

图9 附件影长随时间变化图

由图可以看出，影长在北京时间14:42至15:42时间段内呈上升趋势。

给定杆长3米，作南北回归线和赤道上一天内影长随地方时变化的曲线图如下：

9月3日1米杆一天的影长

影长/m

时间/h

图10 一天内影长随地方时变化示意图

由图可知，不论在什么纬度，影长随时间推移而上升的情况都出现在地方时12时至18时。

将两个图进行对比，可以知道附件中的北京时间段对应在当地时间的12时至18时时段。因为北京时间和当地时间存在如下转换关系：

tT(1200)

4min

60

当北京时间14:42与本地时12时重合时，经度有最小值；令t12:00，T14:42，得79.50；

当北京时间15:42与本地时18时重合时，经度有最大值；令t18:00，T15:42，得154.50。

所以，拍摄地的经度范围是(E79.50,E154.50)。  基于最小二乘法的非线性参数拟合模型

为避免利用附件数据解方程求未知参数所带来的数据利用率不够，精度不高的影响，我们采用非线性拟合的方法求解未知参数。为衡量拟合效果，运用最小二乘法的思想，将理论计算出的x、y坐标值经坐标转换后与数据提供的的x''、y''值之差的平方和作为目标函数，求其最小值，将问题转化成一个优化模型。

值得注意的是，因为拍摄日期为4月18日，则n=28 ，所以太阳直射点纬度选择如下计算式：

1800

arcsin[0.39775sin(n)]

186

利用问题一模型写出具体表达式如下：

minz(x'ix''i)(y'iy''i)2

2

i0

i0

2020

xi'xicosyisinyi'yicosxisinxi

cossin

h

coscoscossinsinsincoscoscossinyih

coscoscossinsin

(7)

1800

arcsin[0.39775sin(n)]

186

(12t)15otT(1200)

4min60

(79.50,154.50)

其中，,,,h为未知量，T,n,x''i,y''i为已知量。 4.2.2模型的求解

 L-M最优化算法[8]

本题的最优化模型是一个利用最小二乘法原理对多个参数进行拟合的问题。经查阅资料我们发现L-M（Levenberg-Marquardt）算法在解决这类问题时有较好的应用。

Levenberg-Marquardt

算法是介于牛顿法与梯度下降法之间的一种非线性优化方

法，这个算法对于过参数化问题不敏感，能有效处理冗余参数问题，而且寻优速度较快。

在本题中应用该算法的流程如下：

目标：对影子顶点坐标关于四个参数的函数关系(x',y')f(,,h,)，给定观测值坐标(x''i,y''i)，估计参数向量p'(,,h,)。

Step1：取初始向量p0(0,0,h0,0)，终止控制常数，计算0||(x',y')f(p0)||，

k:0,0103,v10（也可以是其他大于1的数）。

Step2：计算Jacobi矩阵Jk，计算NkJTkJkkI，构造增量正规方程

NkkJTkk。

Step3：求解增量正规方程得到k。

（1） 如果||(x',y')f(pkk)||k，则令pk1pkk，若||k||，停止迭代，

输出结果；否则令k1k/v，转到Step2。

（2） 如果||(x',y')f(pkk)||k，则令k1kv，重新解正规方程得到k，

返回步骤（1）。

在L-M算法中，每次迭代式寻找一个合适的阻尼因子，当很小时，迭代步长为牛顿法步长；很大时，蜕化为梯度下降法的最优步长计算式。

由于在计算中，给定未知参数的初始值非常重要，求解的结果对初始值较为敏感。而且随着未知参数数量的增加，算法求解的速度会变慢。 所以我们考虑先减少未知数的个数提高算法速度，得出一些精确度不高的参数解。然后再根据得出的解缩小参数范围，带入原模型（不减少参数个数）求解，得出精度较高的参数值。

 L-M最优化算法的改进 Step1：减少未知参数个数

在本问中，由于数据所给的坐标系方向不确定，所以我们设为数据坐标系和问题 一中设定坐标系的夹角，相当于增加了一个未知参数。这样会降低算法求解和速度。所以我们考虑直接从影长入手，将式（7）模型中的目标函数改为：

minz2 （8）

i020

即以理论计算所得影长和所给数据影长之间差值的平方和作为优化目标。因为无论坐标系如何偏转，

影长l都是定值。这样在方程中就不会用到偏转角，可以减少一个未知量，加快算法计算速度。

值得注意的是，用影长进行优化相当于只运用了题目所给的一组信息，降低了对题

目所给信息的利用率，所以结果的精度不高。

这样得出的参数向量p'(,,h)=（108.63o,19.09o,2.03m），拟合匹配的效果如下：

所给数据与计算位置得到结果影长对比图

影长

时间

图11 影长与数据对比图

此时，偏差平方和z5.305108，拟合效果较好。 Step2：缩小参数范围并求解

根据第一步求出的结果，分别增减20o缩小经度的范围为：（80o，120o）,纬度

的范围为：（0o,40o），利用L-M算法求解式（7）模型，即将理论计算出的x、y坐标值经坐标转换后与数据提供的的x''、y''值之差的平方和作为目标函数，结果为

p(,,h,)=（109.83o,18.15o,1.96m，0.26弧度），偏差平方和z3.07106。x值和y值拟合匹配的效果如下：

所给数据与计算位置得到结果x轴坐标对比图

x轴坐标

时间

图12 x值与数据对比图

所给数据与计算位置得到结果y轴坐标对比图

y轴坐标

时间

图13 y值与数据对比图

由图像可知，求得的结果与数据所给的x值和y值符合的效果都较好，且偏差平方和

z3.07106，到达106数量级，说明结果比较理想。 4.3问题三 4.3.1模型的建立

问题三附件中给出了两组数据。两组的时间相差不到30分钟，但是一组中x坐标全为正值，另一组全为负值，所以推测两组数据并不是在同一个地点。

问题三与问题二类似，不同之处在于增加了日期n这个未知参数，整体的思路是一致的。首先仍是进行坐标系的转化，具体过程与问题二中相同，此处不再赘述。

 拍摄地经度范围的估计 1、附件2

图14 附件2影长随时间变化

由图可以看出，影长在北京时间12:41至13:41时间段内呈下降趋势。

由图10可知，不论在什么纬度，影长随时间推移而下降的情况都出现在地方时6时至12时。

将两个图进行对比，可以知道附件中的北京时间段对应在当地时间的6时至12时时段。因为北京时间和当地时间存在如下转换关系：

tT(1200)

4min

60

当北京时间12:41与本地时6时重合时，经度有最小值；令t06:00，T12:41，得19.750；

当北京时间13:41与本地时12时重合时，经度有最大值；令t12:00，T13:41，得94.750。

所以，在附件2中拍摄地的经度范围是(19.750,94.750)。 2、附件

3

图15 附件3影长随时间变化

由图可以看出，影长在北京时间13:09至14:09时间段内呈上升趋势。

由图10可知，不论在什么纬度，影长随时间推移而上升的情况都出现在地方时12时至18时。

将两个图进行对比，可以知道附件中的北京时间段对应在当地时间的12时至18时时段。因为北京时间和当地时间存在如下转换关系：

tT(1200)

4min

60

当北京时间12:09与本地时12时重合时，经度有最小值；令t12:00，T13:09，

得102.750；

当北京时间13:09与本地时18时重合时，经度有最大值；令t18:00，T13:09，得167.250。

所以，在附件2中拍摄地的经度范围是(102.750,1800)(1800,167.250)。  基于最小二乘法的非线性参数拟合模型

由于问题三和问题二的解决思路相同，也需要建立基于最小二乘法的非线性参数拟合模型。但是由于值得注意的是，因为拍摄日期未知，所以太阳直射点纬度选择计算式时需要分三种情况考虑：

1800

arcsin[0.39775sin(n)]，n[0,186]

186

arcsin[0.39775sin(n186)0]，n[186,276] 900

arcsin{0.39775cos[(n276)]}，n[276,365]

89

利用问题一公式写出本问模型具体表达式如下：

minz(x'ix''i)(y'iy''i)2

2

i0

i0

20

20

xi'xicosyisinyi'yicosxisinxi

cossin

h

coscoscossinsinsincoscoscossinyih

coscoscossinsin

(9)

(12t)15o

4min

tT(1200)

60

1800

arcsin[0.39775sin(n)],n[0,186]

186

arcsin[0.39775sin(n186)0],n[186,276]900

arcsin{0.39775cos[(n276)]},n[276,365]

89

(19.750,94.750),(附件2)

(102.750,1800)(1800,167.250)(附件3)

其中，,,,h,n为未知量，T,x''i,y''i为已知量。

4.3.2模型的求解

与问题二处理类似，由于L-M算法对于过参数化问题不敏感，我们仍然选择L-M算法进行求解，并进行了和问题二中类似的优化。

Step1：减少未知参数个数

直接从影长入手，将式（9）模型中的目标函数改为：

minz2

i0

20

即以理论计算所得影长和所给数据影长之间差值的平方和作为优化目标。方程中不用引入偏转角，减少量一个未知量，加快算法的计算速度。

值得注意的是，用影长进行优化相当于减少了对题目所给信息的利用率，所以结果

的精度不高。

1、对于附件2，得出的参数向量p(,,h,m)=（79.55o,39.83o,1.99m，202.19）

其中，m表示距离元旦的天数，即mn79；拟合匹配的效果如下：

影长

时间

图16 附件2影长与数据对比图

此时，偏差平方和z9.337108，拟合效果较好。

2、对于附件3，得出的参数向量p(,,h,m)=（111.06o,36.1o,2.69m，306.47） 其中，m表示距离元旦的天数；拟合匹配的效果如下：

所给数据与计算位置得到结果x轴坐标对比图

时间

图17 附件3影长与数据对比图

此时，偏差平方和z8.887106，拟合效果较好。

Step2：确定参数范围并求解

1、对于附件2，根据第一步求出的结果，分别增减20o缩小经纬度的范围，同时考虑到(19.750,94.750)的限制，得出经度（60o，95o）；纬度（0o,40o），利用L-M算法求解式（9）模型，即将理论计算出的x、y坐标值经坐标转换后与数据提供的的x''、

y''值之差的平方和作为目标函数，结果为

p1'(,,h,,m)=（79.42o,41.67o,2.05m，-0.33弧度，185.76），偏差平方和

z2.07106；

p2'(,,h,,m)=（79.59o,40.59o,2.02m，-0.34弧度，147.91），偏差平方和

z4.26107；

p3'(,,h,,m)=（79.73o,39.67o,1.99m，-0.35弧度，203.03），偏差平方和

z4.69108；

偏差平方和最小的解为p3'，拟合效果最好。它的x值和y值拟合匹配的效果如下：

所给数据与计算位置得到结果x轴坐标对比图

x轴坐标

时间

图18 附件2 x值与数据对比图

所给数据与计算位置得到结果y轴坐标对比图

y轴坐标

时间

图19 附件2 y值与数据对比图

由图像可知，求得的结果与数据所给的x值和y值符合的效果都较好，结果比较理想。 因此，拍摄点的坐标为（E79.73o,N39.67o），拍摄日期大约为7月23日。

2、对于附件3，根据第一步求出的结果，分别增减20o缩小经纬度的范围，同时考虑到(102.750,1800)(1800,167.250)的限制，得出经度（90o，130o）；纬度（10o,60o），利用L-M算法求解式（9）模型，结果为：

p1'(,,h,,m)=（110.7o,34.05o,2.81m，0.17弧度，309.55），偏差平方和

z1.12105；

p2'(,,h,,m)=（110.9o,35.37o,2.77m，0.16弧度，306.46），偏差平方和

z1.97105；

p3'(,,h,,m)=（110.61o,33.37o,2.83m，0.17弧度，311.20），偏差平方和

z7.91106；

偏差平方和最小的解为p3'，拟合效果最好。它的x值和y值拟合匹配的效果如下：

所给数据与计算位置得到结果x轴坐标对比图

x轴坐标

时间

图20 附件3 x值与数据对比图

所给数据与计算位置得到结果y轴坐标对比图

y轴坐标

时间

图21 附件3 y值与数据对比图

由图像可知，求得的结果与数据所给的x值和y值符合的效果都较好，结果比较理想。因此，拍摄点的坐标为（E110.61o,N33.37o），拍摄日期大约为11月8日。 4.4问题四 4.4.1模型的建立

 图像的提取

问题四给出的是一段时长为40分钟的视频，是连续的。为了提取出其中影长的数

据，就必须先从视频中提取若干帧图像。通过观看视频，我们发现在整个过程中影长并没有发生突然或剧烈的变化，所以选择按时间均匀提取。因为视频时长为40分钟，所以我们每隔2分钟提取视频中图像，共提取20张图像，即20组影长数据。

 杆长与影长比值的获取

在视频过程中，我们发现影子旋转的角度很小，所以我们认为在整个过程里影子只 有长度的变化，没有角度的变化。

由于不能确定影子是否与拍摄的正面相平行，所以设它们之间的夹角为，则影长

l在拍摄正面上的投影长度为lcos。

因为已知杆长为2米，为得到影长，就需要知道杆长与影长比值。我们将图像数字化后，利用杆顶点的像素坐标和影子顶点的像素坐标求得它们在图像中的长度比。值得注意是，利用像素点获得的影长实际上是影子在拍摄正面的投影长度，所以



h

（10）

lcos

这样就可以得到所提取的20张图像中的影长实际长度数据，如下表所示：

表1 lcos随时间变化

T(Min) 8.90 8.93 8.97 9.00 9.03 9,07 9.10 9.13 9.17 9.20

lcos 2.37 2.34 2.32 2.29 2.25 2.22 2.19 2.17 2.14 2.11 T(Min) 9.23 9.27 9.30 9.33 9.37 9.40 9.43 9.47 9.50 9.53

lcos 2.08 2.05 2.03 2.00 1.97 1,94 1.92 1.89 1.86 1.84

 已知拍摄日期的模型

因为可以得出杆长和影长的比值，且杆长已知，则可以计算出拍摄地从8:54到 9:34,40分钟内的20组时刻下的影长数据。又因为拍摄日期已知，这是模型蜕化成问题二中Step1以影长偏差为目标函数的非线性参数拟合模型，只是多加入了一个未知量；

值得注意的是，因为拍摄日期为7月13日，则n= 114，所以太阳直射点纬度选择如下计算式：

1800

arcsin[0.39775sin(n)]

186

则模型为：

minz(li

2

i0

20

i

cossin

h

coscoscossinsinsincoscoscossinyih

coscoscossinsin

1800

arcsin[0.39775sin(n)]

186

(12t)15oT(1200)

4min60

（11）

h

lcos

与问题二类似，这本质上还是一个利用已知数据根据最小二乘法原理进行非线性拟



合求解未知参数的问题。由于L-M算法对于过参数化问题不敏感，能有效处理冗余参数问题，所以我们仍采用L-M算法进行求解。求解得结果为p(,,)=（111.140，39.390，0.14弧度）。拟合的效果如下：

影长

时间

图22 影长与数据对比图

由图像可知，求得的结果与图像所给影长符合的效果都较好。因此，拍摄点的坐标为（E111.140，N39.39o）。

 拍摄日期未知的问题模型

拍摄日期未知的问题的拍摄日期已知模型非常类似，只是增加了日期这个未知量， 其他方面都基本相同。因为拍摄日期不确定，所以太阳直射点纬度有三个不同情况下的计算公式，具体模型如下：

minz(li2

i020

xi

cossin

h

coscoscossinsinsincoscoscossinyih

coscoscossinsin

1800

arcsin[0.39775sin(n)],n[0,186]

186

arcsin[0.39775sin(n186)0],n[186,276]

900

arcsin{0.39775cos[(n276)]},n[276,365]

89

(12t)15o

4min

tT(1200)

60

h

lcos

（12）

由于拍摄日期未知，所以本问题与问题三模型类似，仍采用L-M算法进行求解。求解得结果为p(,,,m)=（110.570，39.600，0.008弧度，190.61）。拟合的效果如下：

影长

时间

图23 影长与数据对比图

由图像可知，求得的结果与图像所给影长符合的效果都较好。拍摄日期距元旦190天，即7月10日，与实际的拍摄日期7月13日只相差了三天，结果较理想。因此，拍摄点的坐标为（E110.570，N39.60o），拍摄日期为7月10日。

五、模型的检验

为了检验模型的正确性，我们选取第一问中的地点北京天安门，已知杆子长度为3米，日期为10月22日，地理位置为北纬39度54分26秒,东经116度23分29秒，根据第一问的模型，可以计算得出北京时间14:00到15:00每个时刻的影长，如图7所示； 5.1对问题二模型的检验

以天安门经纬度为待求的未知量，杆长也未知，拍摄日期为已知量，利用第一问模型求得的影长在不同时刻的数据进行非线性参数拟合，即带入问题二的模型进行求解。

Step1：减少未知参数个数

不引入坐标系偏转角，直接从影长入手，以理论计算所得影长和所给数据影长之间差值的平方和作为优化目标。运用L-M算法得出的参数向量p'(,,h)=（116.03o,38.74o,3.05m），拟合匹配的效果如下：

影长

时间

图24 影长与数据对比图

Step2：确定参数范围并求解

根据第一步求出的结果，引入坐标系偏转角，分别增减20o缩小经纬度的范围，

得出经度（100o，150o）；纬度（20o,60o），利用L-M算法求解式（8）模型，结果为p'(,,h,)=（115.920,38.28o,3.10m，0.01弧度），

它的x值和y值拟合匹配的效果如下：

所给数据与计算位置得到结果x轴坐标对比图

x轴坐标

时间

图25 x值与数据对比图

所给数据与计算位置得到结果y轴坐标对比图

y轴坐标

时间

图26 y值与数据对比图

由图可知，x,y值的拟合效果较好。因为116o23’=116.383o,所以经度的相对误差

116.3830115.920oo

为：||100%0.398%；又3954’=39.9 ,所以纬度的相对误差为：0

116.38339.9038.280||100%4.06%；

39.90

由上述计算可知，利用问题二的模型计算所得的经纬度相对误差都较小，说明模型的精确度较高。 5.2对问题三模型的检验

以天安门经纬度、拍摄日期为待求的未知量，同时杆长未知，利用第一问模型求得

的影长在不同时刻的数据进行非线性参数拟合，即带入问题三的模型进行求解。

Step1：减少未知参数个数

不引入坐标系偏转角，直接从影长入手，以理论计算所得影长和所给数据影长之间差值的平方和作为优化目标。运用L-M算法得出的参数向量p'(,,h,m)=（116.63o,40.66o,2.96m，290.77）（其中，m表示距离元旦的天数），拟合匹配的效果如下：

影长

时间

图27 影长与数据对比图

Step2：确定参数范围并求解

根据第一步求出的结果，引入坐标系偏转角，分别增减20o缩小经纬度的范围，

得出经度（100o，150o）；纬度（20o,60o），利用L-M算法求解式（8）模型，结果为p'(,,h,,m)=（116.420,39.96o,3.00m，0.00弧度，291.75），经换算后日期为10月20日。

它的x值和y值拟合匹配的效果如下：

所给数据与计算位置得到结果x轴坐标对比图

x轴坐标

时间

图28 x值与数据对比图

所给数据与计算位置得到结果y轴坐标对比图

y轴坐标

时间

图29 y值与数据对比图

由图可知，x,y值的拟合效果较好。因为116o23’=116.383o,所以经度的相对误差

116.3830116.420oo

为：||100%0.032%；又3954’=39.9 ,所以纬度的相对误差为：0

116.38339.9039.960||100%0.15%；利用问题三模型计算所得的日期为10月20日，与实

39.90际日期只相差2天。

由上述计算可知，利用问题三的模型计算所得的经纬度和日期的误差都较小，说明模型的精确度较高。

六、模型的优点和缺点

 优点

1、模型采用L-M算法对过参数化问题不敏感，能有效处理冗余参数问题，而且寻优速度较快。

2、通过模型检验可以发现利用模型推导的值与真实值得偏差很小，说明模型较为准确，有一定的精度。

 缺点

1、L-M算法对初始值较为敏感，不同初始值的选取可能会给结果带来偏差。

2、在问题四中确定影长和杆长的比值的方法是人工寻找像素点，效率较低，且存 在误差。

参考文献

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9] 汪和平.探究日影运动轨迹.中学数学学刊. 2010年 屈名，王征兵等.基于交比不变性的太阳定位算法的研究.SILICON VALLEY. 2013年. 林根石. 利用太阳视坐标的计算进行物高测量与定位.南京林业大学学报. 1991年. 张文华，司德亮等，太阳影子倍率的计算方法及其对光伏阵列布局的影响. 蒋洪力.太阳直射点纬度的数学推导和分析.学通报. 2007年. 黄玫生等.利用公式求算太阳直射点和极昼极夜的纬度位置. 潘志荣.太阳视运动路径及物体影子_问题解析.发明与创新. 2007年. 卓金武等.MATLAB在数学建模中的应用. 北京:北京航空航天大学出版社, 2014.9 赵怀森，赵永熹，王宝阔.任意海区的太阳移线定位时机判断研究.中国航海.2004年。

附录

clear all;

num=xlsread('附件1-3.xls','附件1');

%num=xlsread('问题一.xlsx','坐标');

hour=num(:,1)*24;

%y(:,1)=num(:,2);% x坐标

%y(:,2)=num(:,3);% y坐标

y=sqrt(num(:,2).^2+num(:,3).^2);

p0=[110,20,2];

%options = optimset('MaxFunEvals',4000);

[result,resnormal,exitflag]=lsqcurvefit(@qiul,p0,hour,y,[-180,-90,0],[180,90,7]);

test=qiul(result,hour);%用拟合参数结果得到该店的理论阴影

plot(hour,y(:,1),hour,test(:,1),'r');%画出结果图和所给数据对比图x值 title('所给数据与计算位置得到结果x轴坐标对比图');

xlabel('时间');

ylabel('x轴坐标');

legend('所给数据','模型结果','location','best');

clear all;

num=xlsread('附件1-3.xls','附件1');

%num=xlsread('问题一.xlsx','坐标');

hour=num(:,1)*24;

y(:,1)=num(:,2);% x坐标

y(:,2)=num(:,3);% y坐标

p0=[100,25,3,0];

[result,resnormal,exitflag]=lsqcurvefit(@qiuxy,p0,hour,y,[-79.5,-90,0,-pi],

[154.5,90,5,pi]);

test=qiuxy(result,hour);%用拟合参数结果

plot(hour,y(:,1),hour,test(:,1),'r');

title('所给数据与计算位置得到结果x轴坐标对比图');

xlabel('时间');

ylabel('x轴坐标');

legend('所给数据','模型结果','location','best');

figure,plot(hour,y(:,2),hour,test(:,2),'r');%画出结果图和所给数据对比图y值 title('所给数据与计算位置得到结果y轴坐标对比图');

xlabel('时间');

ylabel('y轴坐标');

legend('所给数据','模型结果','location','best');