2015年全国数学建模大赛一等奖

赛区评阅编号（由赛区组委会填写）：

2015高教社杯全国大学生数学建模竞赛

承诺书

我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。

我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号（从A/B/C/D中选择一项填写）： A

我们的报名参赛队号（12位数字全国统一编号）：201508011076

参赛学校（完整的学校全称，不含院系名）：哈尔滨理工大学

参赛队员(打印并签名) ：1. 鲁庆豪

2. 孙根

3. 姚朝霞

指导教师或指导教师组负责人(打印并签名)：郑红艳

日期： 2015 年 9 月 13 日

（此承诺书打印签名后作为纸质论文的封面，注意电子版论文中不得出现此页。以上内容请仔细核对，特别是参赛队号，如填写错误，论文可能被取消评奖资格。）

赛区评阅编号（由赛区组委会填写）：

2015高教社杯全国大学生数学建模竞赛

编号专用页

送全国评奖统一编号（由赛区组委会填写）：

全国评阅统一编号（由全国组委会填写）：

此编号专用页仅供赛区和全国评阅使用，参赛队打印后装订到纸质论文的第二页上。注意电子版论文中不得出现此页，即电子版论文的第一页为标题和摘要页。

太阳影子定位

摘要

太阳影子定位技术是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。本文按照题目要求，利用太阳影子变化规律，建立太阳影子长度变化模型，循序渐进地解决了从第一问到第四问，基本解答了太阳影子定位问题。

针对问题一，考虑到影子主要随着太阳高度角的变化而变化的规律，利用地球物理及三角几何知识建立了有关影子长度关于太阳高度角的数学模型。应用此模型，较容易地求得物体在一段时间内的长度，最后通过MATLAB软件作图，画出了2015年10月22日北京时间9:00-15:00之间天安门广场3米高的直杆的太阳影子长度的变化曲线。

针对问题二，考虑到影子增长趋势的变化，首先通过建立曲线拟合模型模拟出影子变化的曲线，找出影长最短的时间点，通过该时间点的当地时间与北京时间的关系，计算得到当地经度。接着对影子顶点的数据进行几何分析，找出杆高与太阳时角的变化关系，依据几何模型求解杆高。在已知当地经度与杆高的基础上，通过建立太阳高度角与观测点纬度等变量的非线性方程模型，求解观测点纬度。最后得到附件1中数据的测量地点为东经112度，北纬23度（位于广东省云浮市）。

针对问题三，观测点的经度，杆高的求解仍使用问题二的拟合模型。但此时由于观测日期未知，使得太阳高度角的计算公式有两个未知量，考虑使用穷举法。可通过对观测日期进行穷举，使用VC++软件编写代码，通过多次代入数据穷举，找出满足方程的值，从而找到满足条件的地点。最后得出：附件二所示地点：东经74.25度，北纬40度，日期：3月22日；附件三所示地点：东经109.5度，北纬32.7度，日期：9月20日。

针对问题四，先使用Matlab从视频中导出图片，并将这些图片导入SketchUp进行配图分析，得到真实的随时间变化的影长。然后分别利用之前求解非线性方程和穷举的方法，找出已知日期的拍摄地点和未知日期的情况下的拍摄地点与日期。最后得到可的拍摄地点为东经115度，北纬40.5度与东经115度北纬41.1度。

关键词：太阳影子定位，SPSS最小二乘法拟合，穷举法，SketchUp照片配图分析

一、问题重述

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

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

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

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

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

二、问题的分析

此题主要是利用太阳影子变化规律确定直杆的观测时间及地点的数学建模问题。进而，利用此太阳影子定位技术对视频中提取来的直杆太阳影子变化数据进行分析，最终确定若干可能的视频拍摄时间及地点。

对于问题一，太阳高度角、被测直杆的长度、观测点的地理位置、日期以及时间等因素直接影响太阳影子长度的变化，为分析影子长度关于各参数的变化规律，可利用地球物理及三角几何知识建立有关影子长度的数学模型。应用此模型，在观测日期、地理位置以及直杆长度确定的情况下，能够较容易地求出一段时间内直杆的影子长度，及画出该时间段内影子长度的变化曲线图。

对于问题二，杆高的计算：虽然题中没有给出杆高的数值，各个量也在不断变化，但可以根据几何关系粗略估测出杆高，以简化之后计算。经度的计算：由问题1中对图2影长变化曲线的观察所得结论可知，该曲线具有二次曲线的普遍特征：具有对称轴，在对称轴处取得极值。因此由附件所给数据可以拟合出影长随时间变化的曲线。又因为影长取得最小值时为当地真太阳时正午，同时已知平太阳时，根据(10)式即可求得经度。纬度的计算：关于当地纬度的计算，由于很难根据影长的变化计算解出纬度，思考根据经度，纬度，杆高与时角，太阳赤纬的数量关系，将时角，赤纬分别用经度纬度的函数关系式表示出来，进而根据第一问的公式，表示出太阳高度角与当地经度，纬度以及其他变量之间的函数关系式，将已经求出的经度，杆高和题目中已知的北京时间，影长等变量代入其中，求解非线性方程，继而解出当地纬度。

对于问题三，问题三的思路与问题二大体相近，分别通过曲线拟合计算经度，利用几何方法计算杆高，最后通过太阳高度角的公式解出纬度以及时间。区别在于最后一步，需要同时解出纬度以及时间。

由于通过解方程组的方法涉及正反三角的运算，即使借助计算机，也难以得到较为准确的结果，所以解方程的方法不可为。这里转换方法，由于本题参数中的日期非连续，且参数均有着确定的区间且区间不大，利用计算机计算快的特点，使用穷举算法进行求解，取得统一较为准确的答案。

对于问题四，只要得出视频中的影子顶点的坐标和长度，即可将问题四的两小问分别转化为问题二和问题三。本问重点是将视频中影子顶点的位置信息转化为真实的坐标和长度。

由于存在透视畸变，依赖于相机和目标物体距离的目标场景像素值会产生很大变化，使得直接测量影子的长度存在很大的误差。这里考虑将视频导出，得到若干张图片后，使用软件进行照片匹配分析，建立三维坐标，将存在畸变的二维图片转化为真实的三维坐标，即可得到真实的影长。

三、模型假设

1、假设太阳折射、海拔等因素对太阳高度角没有影响。

2、假设各测量地点直杆与水平地面垂直。

3、假设视频中各水泥砖边沿平行或垂直且水泥砖左侧边沿与斜坡底边垂直。

4、假设视频中整个水泥砖地面水平。

5、忽略太阳光在大气中的折射率。

四、符号说明

符号表示含义

E表示（东经为正数，西经为负数）

?表示（北纬为正数，南纬为负数）

t表示北京时间（平太阳时）

T表示当地时间（真太阳时）

L表示直杆影子长度

H表示直杆的长度

h表示太阳高度角

s

δ表示太阳赤纬

ω表示太阳时角

五、模型的建立与求解

5.1 问题1——基于太阳影子定位技术的影长变化模型

5.1.1 太阳影子长度模型的建立

影子的长度除了和直杆本身的长度有关之外，还和太阳光的入射角度有关。太阳高度角是指太阳光的入射方向和地平面之间的夹角。而太阳高度角与测量地的纬度、太阳时角及太阳的赤纬度有关。因此设L 为影子长度，H 为直杆长度，s h 表示太阳高度角，

?表示纬度，ω表示太阳时角，δ表示太阳赤纬。由几何关系知

),,(cot δ?ωHf h H L s ==. （1） 太阳时角由测量地的当地时间，即真太阳时决定。若T 表示真太阳时，则太阳时角的计算公式可表示为

)12(15-?=T οω （2）

太阳赤纬是地球赤道平面与太阳和地球中心的连线之间的夹角，即当前太阳直射地的纬度。显然随着日期与季节的不同，同一地点的太阳赤纬不同。通过查找文献[1]知，太阳赤纬角有Cooper 算法，Spencer 算法，Bourges 算法等。可虑到算法复杂度以及硬件实现问题，这里采用较为简单的公式进行粗略计算[2]。由Perrinde Brichambaut 导出： ()()17398563.0cos 39795.0sin -=N δ （3） 其中N 表示日序数，从每年的1月1日开始计算。

太阳高度角的计算则根据文献[3]，观测点对太阳视运动高度角的描述可以表示为向量n 和向量l 夹角的预交，

图1 太阳高度角示意图

令太阳高度角为s h 。由于H 平面完全通过y O -轴，因此易得其平面方程为： ?tan ?-=z x （4） 由上式可得平面H 的单位法向量n 为:

)sin ,0,(cos ??= （5） 根据时角和太阳赤纬可得太阳直射光线的单位向量l 为：

()δδωδωsin ,sin sin ,cos cos =l （6） 由向量的夹角公式得：

>=

<,cos （7）

由于：

><=l n s ,cos sinh （8） 结合式（6）与式（8）得太阳高度的最终表达式：

δ?δ?ωsin sin cos cos cos sin +=s h （9） 至此，建立了完整的太阳定位模型。各参数分别由式(1)，(2)，(3)，(9)给出。

5.1.2 影长变化模型及其关于参数的变化规律分析

由(1)式，结合三角函数变换公式1sin 1

cot 2

-=h

h s ，代入（9）中，就得到了影子长度变化的数学模型为： )1)sin sin cos cos (cos 1

(

2

-+=δ?δ?ωH L (10)

由该模型对影子长度关于各个参数的变化规律如下： 第一，其他变量相同时，影子长度和直杆长度成正比。

第二，太阳时角ω越小，L 越小，反映在初始数据上则为越接近中午（本地时间12点），影子长度越小。

第三，变量太阳赤纬与测量地纬度对影子长度的影响应该分情况而论，其他变量不同，太阳赤纬与测量地纬度对影子长度的影响的正负性也不相同。

5.1.3 影长变化模型的求解

题中情景各项数据已经给出，如下所示：

N=294，H=3，''2454'39ο=?，]15,9[∈t

其中t 为“北京时间”为东八区的地区时，即平太阳时。真太阳时与平太阳时的转换关系如下，其中E 为经度，东经为正，西经为负：

15

120ο

-+=E t T （11）

将(11)代入(3)(8)式，可得出太阳高度角h 关于平太阳时t 的表达式：

1250.0-180)176-15cos 7514.0sin ??

? ??

=πt h s （ (12)

结合(10)(12)式，可得出10月22日北京时间9:00-15:00之间天安门广场3米高的直杆的太阳影子长度的变化函数： 11250.0-)180

)176-15cos(7514.0(1

3

2

-=π

t L （ （13）

通过Matlab 软件作出直杆的太阳影子长度的变化曲线如图2所示：

时间/h

太阳影子长度/m

图2 影长—时间变化曲线图

通过对该曲线的观察可知，影长在真太阳时的正午时刻取得最小值，而且上下午的影长变化关于此时刻对称。

5.2问题2——

5.2.1模型建立

(1)杆高计算的几何模型

图3 杆高计算的几何示图

由图3，在附件2中任取两组数据，得到各自的横纵坐标),(),(2211y x y x ，图示影长21,l l ，影子间的夹角θ?。φ?在误差允许的范围内可以看作太阳的时角变化，即每1个小时对应ο15=?φ.这样可以得到如下等量关系：

H l H l y y x x H ++++=θ （14）

通过代入数据影子顶点横纵坐标数据即可解得杆高H 。 (2)计算经度的曲线拟合模型

由于附件1中所给的数据较多但没有取到最值，由问题分析，通过最小二乘法进行二次函数拟合影长与平太阳时的关系：

c bt at t L ++=2)( （15）

由此得到L 的最小值对应时间a

b

t 2-min =.该时刻对应本地中午时间12：00。由（11）式解得经度计算结果为

ο

ο120)(15min +-=t T E 中午 （16）

(2)关于纬度计算的非线性方程模型

首先确定本文中太阳高度角s h ，影长L ，太阳时角ω，纬度?以及太阳赤纬δ与各变量之间的关系：

1.利用三角函数关系，解得,sinh 2

22

L H H s +=其中H 表示杆高，L 表示影长。 1.设当地经度为E ，结合北京时间经度为E ο120，求出当地时间，进而求出当地时间所对应的太阳时角E t E

t +-=---

?=30015)1215

120(15ω，转化为弧度制为180

)30015(π

ω?

+-=E t ，其中t 为当地所对应的北京时间，即数据中所给出的对应的北

京时间。

2.纬度?转化为弧度制为180

π

??=

y 。其中y 表示纬度（角度制）。 3.太阳赤纬的计算，利用第一问中的公式))173(98563.0cos(39795.0sin -?=N δ，其中N 表示从每年的1月1日起开始计算的日数。

将上文中所求出的太阳高度角，时角，纬度以及太阳赤纬代入公式

δ?δ?ωsin sin cos cos cos sinh +=s 中，得出

δπ

δππsin 180

sin cos 180cos )180)30015cos((2

22??+????+-=+y y E t L H H (17) 式中H E ,由上文模型可求出，s h t ,为已知量，因此上式可转变成关于纬度y 的一元非线性方程，代入数据即可解出纬度。

5.2.2问题二模型的求解

（1）求解关于精度的曲线拟合模型

根据22y x L +=，先求出不同时间所对应的影长。由于测量时刻大都处下午2时附近，因此t 以下午2时为起点，单位为分钟。这里使用SPSS 软件通过最小二乘非线性拟合，可以得到一天内的影长变化。拟合得到的二次函数如下：

786.0007.010136.4)(25++?=-t t t L (18)

拟合图像见图4，通过散点的位置，可见拟合效果比较好。接着，求出影长的最小

值对应，为（分钟）6.84min -=t ，即为此时当地的真太阳时。将结果代入（16）式，得到经

度约为E 112ο。

图4 影长拟合曲线

（2）求解关于杆高的几何模型

在附表1中取多组数据带入（18）式中，计算并求取平均值得m L 1.75≈. （3）求解关于纬度的非线性方程模型

由上文已经求出的当地经度为E ο111，结合北京时间经度为E ο120，求出当地时间，进而求出当地时间所对应的太阳时角18915)1215

111

120(15-=---?=t t ω，为了计算的方便，转化为弧度制为180

)18915(π

ω?

-=t ，并将纬度?也转化为弧度制。且根据公式

（），可知附表中日期4月18日所对应的1740.0sin =δ。

先代入已经确定的数据，得

174

.0180

sin 9847.0180cos )180)18915cos((222??+????-=+π

?π?πt L H H (19) 将L 代入上式，并代入附表1中的时间与对应影长，最终纬度确定为为N.23ο 综上所述，本文确定了直杆所处的地点为N E,23112οο（位于广东省云浮市）。

5.3 问题3

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

（1）经度的求解

对于经度的求解，使用第（2）问的曲线拟合模型，此处不再赘述。

根据22y x L +=，求出附件2与附件3两组数据中不同时间所对应的影长。由于测量时刻大都处下午2时附近，因此t 以下午2时为起点，单位为分钟。这里使用SPSS 软件通过最小二乘非线性拟合，可以得到一天内的影长变化。拟合得到的二次函数如下：

附件2：

25110726.2010.0632.1)(t t t L -?+-= (20) 附件3：25210235.8007.0644.4)(t t t L -?+-= (21) 进而分别求出影长最小值所对应的时间，4.1831=t ，50.422=t ，将结果代入（15）式，可得附件2与附件3中所测地点的经度：E E ο741=，E E ο1092=。

（2）求解关于杆高的几何模型

利用第二问中以求得的几何关系模型，分别取附表二和附表三中的数据代入（16）式中，剔除掉明显有误的计算结果，得到附件二中的杆高结果：1.80m ，附件三的杆高结果：3.34m 。

（3）使用穷举法同时求解纬度与日期

将纬度?转化为弧度制。由上文已经求出的当地经度,109,7421E E E E οο==结合北京时间所对应的经度E ο120，求出当地时间，进而求出当地时间所对应的太阳时角：

75.22515)121525

.74120(151-=---

?=t t ω (22)

5.19015)1215

5

.109120(152-=---?=t t ω (23)

将上文中所求出的各参数代入公式δ?δ?ωsin sin cos cos cos sinh +=s 中，得出

δ

π

?δπ?πsin 180sin cos 180cos )180)5.19015cos((222??+????-=+t L H H (24) 式中太阳赤纬:

))173(98563.0cos(39795.0arcsin(-?=N δ (25)

穷举法的基本思想是根据题目的部分条件确定答案的大致范围，并在此范围内对所有可能的情况逐一验证，直到全部情况验证完毕。若某个情况验证符合题目的全部条件，则为本问题的一个解；若全部情况验证后都不符合题目的全部条件，则本题无解。本题使用穷举法进行计算有两个未知参数的太阳高度角方程。

代入杆长H 后，（23）式中存在两个待定的参数：纬度?、天数N 。输入附表中的时间与影长，便可对这两个量进行穷举求解。为了减少穷举的次数，本题假设测量地在陆地上，通过查看地图，发现E 74ο上的陆地纬度为在N 73-N 14οο之间（显然测量地不在南北极）。可将纬度的区间定位[14,73]，取正数枚举。天数的区间为[1,365]。

本题使用Microsoft Visual C++使用C++语言进行编写。每次输入附件中的一组时间与影长数据，找出大致使（23）式成立的点。这里的“大致”设定为（23）式的左右两边差的绝对值小于设定的值则符合要求并输出（程序界面见图？）。通过代入不同的数据，找出输出频率最大的纬度与天数，作为本题的结果。

图5 穷举算法程序的输入输出界面

综合结果后得到最终结果如下：

附件二：4081==?，N 附件三：33,243==?N 综上所述，可得问题三的最终答案：

附件二所示地点：E N οο25.74,40， 日期：3月22日 附件三所示地点：E N οο5.109,7.32， 日期：9月20日

5.4 问题4——基于SU 图像分析的定位模型

5.4.1 问题分析

只要得出视频中的影子顶点的坐标坐标和长度，即可将问题四的两小问分别转化为问题二和问题三。本问重点是将视频中影子顶点的位置信息转化为真实的坐标和长度。

由于存在透视畸变，依赖于相机和目标物体距离的目标场景像素值会产生很大变化，使得直接测量影子的长度存在很大的误差。这里考虑将视频导出，得到若干张图片后，使用软件进行照片匹配分析，建立三维坐标，将存在畸变的二维图片转化为真实的三维坐标，即可得到真实的影长。

5.4.1 附件处理——基于SU 图像匹配

首先，提取有效信息。对于附件四所给的avi 格式的视频，由于每秒的影长变化可以忽略不计，因此视频中每时每刻的图像中含有大量无效信息。于是将其按照一定时间间隔转化为图片。通过编程将其提取为52张bmp 格式的图片。这些图片就分别记录了每隔一分钟的杆及其影子的瞬时情况，包含了有效的信息。

其次，分析图片信息。本文将这些图片导入SketchUp Pro 进行照片匹配分析。如图

所示：

图6 照片匹配分析图

在图中，以杆子的最底端为坐标中心，建立空间三维坐标系。其中x轴方向由斜坡与地面的交线、右侧洞口上沿共同确定；y轴方向由左侧花坛边缘、以及直杆旁边水泥砖边缘共同确定。

最后，进行长度测量。在确定好地平线之后，使用测量选项即可测出图中各物体的长度。如果要得到绝对长度，需要在图中找到已知长度的参照物。而由题意杆长已知为2米，因此可以求得影长以及太阳高度角，结果如下表所示。

5.4.3 经纬度的确定

同样，按照第(2)问的解法，利用SPSS 软件求解拟合函数。以8：00为原点，单位为分钟，得到拟合函数：

783.2049.010362.9)(25+-?=-t t t L （26）

求出对称轴,得264min =t ，对应的时间为12：34分。由（14）式，可求出经度为E 115ο。 在已知经度和太阳高度角正弦值的情况下，将太阳高度角正弦值、经度、日期代入式（16），由于只有纬度值一个变量，通过MATLAB 软件求解非线性方程即可解出纬度值。然后多次代入太阳高度角正弦值的多组数据，求解多组纬度值，删除明显有误的结果，得到纬度值N ο5.40或N ο1.41

综上所述，得到可能为拍摄地点的位置如下：

E N οο115,5.40与E N οο115,1.41

如果拍摄日期未知，依然可以通过视频确定拍摄地点与日期，首先通过SU 图像匹配首先可以求出在太阳高度角正弦值，并通过影长变化的曲线拟合可以求出经度，杆高已知，接下来的方法与第三问类似，通过穷举法求出纬度与日期。最终得到拍摄地点与日期。

六、模型检验

在根据影长定位或者依据经纬度求得影长的相关研究方面，事实上由于现在的GPS 技术，天文观测技术以及大数据分析技术等的高度发达，已经有了一些智能计算的结果。由于其背后的高科技支撑，计算可谓又快有准。本文因此创新性地提出一种检验方法，与以往的传统模型分析不同，用在线网站[5]进行检验，结果如下。

第一问根据网站得出的结果如下：

图7 问题1影长变化实际曲线

在9点至15点部分与图非常接近，具体数值及误差如下表所示：

本模型计算结果实际情况误差

6.74 6.59 2.2%

5.61 5.50 2.0%

4.85 4.76 1.8%

4.33 4.25 1.8%

3.99 3.92 1.7%

3.79 3.73 1.6%

由杆长为1.75m，得到广东省浮云市地区具体数值及误差如下：

时间数据中影长当地影长误差

14:42 1.15 1.16 0.8%

14:45 1.18 1.19 0.8%

14:48 1.215 1.216 0.08%

14:51 1.25 1.24 0.8%

14:54 1.283 1.278 0.3%

14:57 1.317 1.309 0.6%

时间视频中影长当地影长误差

8:54 2.285 2.235 2.2%

8:55 2.276 2.221 2.4%

8:56 2.265 2.206 2.6%

8:57 2.246 2.192 2.4%

8:58 2.229 2.177 2.3%

8:59 2.213 2.164 2.2%

结果准确度较高。

七、模型的评价

本论文所解决的一系列问题的具体背景是太阳影子定位技术，也就是实现由物体影子的变化与日期及经纬度的相互转化。而由于科学技术的提高，全球定位系统与大数据分析已经取得了较大进展，实际中所解决的问题不可能基于本文所做的方法。从这个意义上讲，本文的模型没有较好的推广性，但这是由于目前数学模型方法本身的局限性决

定的。

本文在第(1)问中基于初等数学的运算推到与实际情况仅仅相差1%，方法简易却得到了很好的结果，相比于其他一些复杂算法而言减少了很大工作量，这种以解决问题为导向而非追求复杂算法的模型值得推广。

本文在处理第(4)问中的视频问题时在没有专业人士的帮助情况下创新性地引入设计软件进行图片处理，最终转化为前两问的处理方法。

参考文献

[1]曾利霞，基于视日运动轨迹的双轴太阳跟踪系统的研究[D]，湖北：湖北工业大学，2012

[2]刘晓燕，余学江，徐巧峰，王海燕，半固定式太阳能集热器角度研究[J]，低温建筑技术，152：116，2011

[3]房淼森，李少华，一种太阳视运动轨迹建模方法及其应用[J]，城市勘测，2015(2):109-112，2015

[4]孙义欣，冯娜，穷举法在程序设计中的应用[J]，计算机时代，2012(8):50-52,2012

程序源码

第一问图

clc

n=[9.00:1/6:15.00];

t=n-12;

oumiga=15*t/180*3.1415926;

theta=2*3.1415926*t/365.2422;

delta=0.3723+23.2567*sin(theta)+0.1149*sin(2*theta)-0.1712*sin(3*theta)-0.7 58*cos(theta)+0.3656*cos(2*theta)+0.0201*cos(3*theta);

delta=delta/180*3.1415926;

fai=(39+54/60+26/3600)/180*3.1415926;

h=asin(sin(fai)*sin(delta)+sin(fai)*cos(delta).*cos(oumiga));

A=asin(cos(delta).*sin(oumiga)./cos(h));

H=3.0;

X=H*cot(h)./sqrt(1+(tan(A)).^2);

Y=tan(A).*X;

L=sqrt(X.^2+Y.^2);

plot(n,L),xlabel('时间/h'),ylabel('太阳影子长度/m')

title('直杆的太阳影子长度的变化曲线');

grid on;

第二问的拟合曲线与散点

x=-250:150;

y=4.136.*10.^(-5)*(x.^2)+0.007.*x+0.786;

plot(x,y)

grid

xlabel('时刻/时');

ylabel('影子长度/米')

hold on

a=xlsread('e:\a');

b=xlsread('e:\b');

plot(a,b,'*')

问题二求解非线性方程代码

syms y;

y=solve('sqrt(1.92*1.92/(1.92*1.92+1.249051*1.249025))=cos((15*14.85-189)*p i/180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.92*1.92/(1.92*1.92+1.149626*1.149626))=cos((15*14.7-189)*pi /180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.92*1.92/(1.92*1.92+1.353364*1.353364))=cos((15*15-189)*pi/1 80)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.92*1.92/(1.92*1.92+1.4634*1.4634))=cos((15*15.15-189)*pi/18 0)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.92*1.92/(1.92*1.92+1.579853*1.579853))=cos((15*15.3-189)*pi /180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.92*1.92/(1.92*1.92+1.703921*1.703921))=cos((15*15.45-189)*p i/180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.92*1.92/(1.92*1.92+1.835014*1.835014))=cos((15*15.6-189)*pi

/180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.17*1.17/(1.17*1.17+1.249051*1.249025))=cos((15*14.85-189)*p i/180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.17*1.17/(1.17*1.17+1.149626*1.149626))=cos((15*14.7-189)*pi /180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.17*1.17/(1.17*1.17+1.353364*1.353364))=cos((15*15-189)*pi/1 80)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.17*1.17/(1.17*1.17+1.4634*1.4634))=cos((15*15.15-189)*pi/18 0)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.17*1.17/(1.17*1.17+1.579853*1.579853))=cos((15*15.3-189)*pi /180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.17*1.17/(1.17*1.17+1.790051*1.790051))=cos((15*15.55-189)*p i/180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

syms y;

y=solve('sqrt(1.17*1.17/(1.17*1.17+1.927918*1.927918))=cos((15*15.7-189)*pi /180)*cos(pi*y/180)*0.9847+sin(y*pi/180)*0.174','y')

问题3穷举算法

#define pi 3.1415926

#include

#include

using namespace std;

int a,b,n,i,p,q,j=0;

float t,l,H;

double c,d,h,f;

int main()

{

cout<<"输入杆长:";

cin>>H;

cout<<"输入纬度的上下限:";

cin>>p>>q;

for(i=0;i<21;i++) //最多重复附件21组数据

{

cout<<"输入时间:";

cin>>t;

cout<<"输入影长:";

cin>>l;

h=sqrt(H*H/(H*H+l*l));

for(a=p;a<=q;a++) //设定纬度范围

{

for(b=1;b<365;b++) //设定日期范围

{

c=sin(0.3975*(cos(0.98563*(b-173))));

d=sqrt(1-c*c);

f=cos((15*t-190.5)*pi/180)*cos(a*pi/180)*d+sin(a*pi/180)*c;

/*cout<

cout<

system("pause");*/ //调试算式结果

if(fabs(f-h)<0.0003) //设定差值

{

cout<<"符合条件的纬度与日期分别为:"<

j++;

}

}

}

if(j<1) {cout<<"不符合"<

}

return 0;

}

问题四求解非线性方程代码：

syms y;

y=solve('0.659=cos((15*8.9-190.62495)*pi/180)*cos(pi*y/180)*0.9847+sin(y*pi /180)*0.174','y')

syms y;

y=solve('0.665=cos((15*8.95-190.62495)*pi/180)*cos(pi*y/180)*0.9847+sin(y*p i/180)*0.174','y')

syms y;

y=solve('0.672=cos((15*9-190.62495)*pi/180)*cos(pi*y/180)*0.9847+sin(y*pi/1 80)*0.174','y')

syms y;

y=solve('0.679=cos((15*9.05-190.62495)*pi/180)*cos(pi*y/180)*0.9847+sin(y*p i/180)*0.174','y')

syms y;

y=solve('0.685=cos((15*9.1-190.62495)*pi/180)*cos(pi*y/180)*0.9847+sin(y*pi /180)*0.174','y')

syms y;

y=solve('0.691=cos((15*9.15-190.62495)*pi/180)*cos(pi*y/180)*0.9847+sin(y*p