太阳影子定位数学建模论文

1、精选文档太阳影子定位摘 要 太阳影子定位对视频拍摄分析至关重要，本文通过建立几何模型、太阳高度角模型和遍历模型，绘制出了影长变化曲线，确定了视频拍摄的地点和日期，解决了太阳影子定位问题。针对问题一，基于地球是球形的假设，建立几何模型和坐标系，用与影长相关的五个参数表示出了杆顶点坐标和光线向量。由点、线、面间的关系求得直线（光线向量所在直线）和平面（过原点且与杆垂直的平面）的方程，联立方程组求得了影子顶点坐标，并用向量的模表示出了影子的长度，由Matlab软件绘制出了影长变化曲线。最终，基于地球的真实外形，用太阳高度角模型对模型进行检验，验证了它的合理性。针对问题二，由于它恰好有两个条件的确定性

2、与问题一相反，所以我们接受问题一精度更高的太阳高度角模型倒推求解。由勾股定理求出21组影长，以影长为纵坐标，北京时间为横坐标作图，得到的最低点同时对应北京时间和当地时间12:00，依据问题一中的时差关系式，反推出当地的经度，再用遍历法，求出了合适的纬度，由经纬度确定了地点为蒙古自治区鄂尔多斯市。针对问题三，由于它与问题二的区分仅是日期未知，所以只需求日期，地点用问题二模型来求。我们接受逆向思维，反推出求日期需要知道五个角，然后正向求解。先通过方位角与影轴角、时角和太阳高度角的关系建立两个等式求出方位角，再利用赤纬角与方位角的关系求出赤纬角，最终利用赤纬角与日期的关系，建立遍历模型求出了日期。最

3、终接受问题二的模型求解地点。针对问题四，我们通过处理图像求得了影长。在日期已知时，我们通过绘制影长-时间图求解出经度为，通过求解高度角与其三个参数的值，建立太阳高度角模型，通过总关系式求出了纬度为。 在日期未知时，我们给出了遍历模型求解纬度，由于经度与日期无关，所以照旧接受日期已知时的方法来求解。本模型考虑了地球不是规章球体的因素，引入了修正值，使结果更加牢靠，且后续问题参照前面的模型来求解，使问题大大简化。关键字：太阳影子定位；几何模型；太阳高度角模型；遍历模型1. 问题重述视频数据分析需要确定视频的拍摄地点和拍摄日期两大方面，太阳影子定位技术就是其中一种确定方法，它通过分析视频中物体太阳影

4、子的变化状况来反推拍摄视频的地点和日期。为了解决太阳影子的定位问题，我们要解决以下四个问题。1.通过设定中间参数，建立影子长度变化的数学模型，并分析出影子长度关于各个参数的变化规律。应用建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。2.由直杆在水平面上的影子顶点坐标数据，建立数学模型，确定出直杆所处的地点。将模型应用于附件1的影子顶点坐标数据，给出若干个符合条件的地点。3.由直杆在水平面上的影子顶点坐标数据，建立数学模型，确定直杆所处的地点和日期。将模型分别应用于附件2和附

5、件3的影子顶点坐标数据，给出若干个符合条件的地点与日期。4.建立确定视频拍摄地点的数学模型，并将其应用于附件4，求出若干个符合条件的拍摄地点。附件4为一根直杆在太阳下的影子变化的视频截图，并且已通过某种方式估量出直杆的高度为2米。在拍摄日期未知的状况下，推断是否可以依据视频确定出拍摄地点与日期，若可以则给出模型和方法。2模型假设1.假设地球是一个规章的球形；2.假设附件所给数据真实可信；3.假设本文引用的定理公式精确度足够高。3通用符号说明序号符号符号说明1直杆的长度2影长3经度4纬度5太阳赤纬：地球赤道平面与太阳和地球的中心的连线之间的夹角6日期：从1月1日开头计时，每天累计+17太阳高度角

6、8时角4问题一模型的建立与求解4.1问题分析问题一要我们建立影子长度变化的数学模型，并分析影长关于各个参数的变化规律。由于在基本假设中已经假设了地球是一个规章的球形，我们打算建立几何模型，接受建立坐标系的方法来求解直杆影子的向量，拟将此坐标系的原点设在坐标原点，赤道设为xoz平面，地球自转轴设为y轴。由于影子的长度与直杆长度、经度、纬度和太阳光线直射角度有关，而太阳直射角度又与日期有关，我们打算设定五个参数分别表示它们。我们拟用这五个参数表示出直杆顶点坐标和太阳光线向量，可能由点、线、面间的关系表示出太阳光线向量所在直线的方程、过原点且与直杆垂直的平面的方程，联立两个方程可能得到影子顶点的坐标

7、。此坐标到原点的距离就是直杆影子的长度，我们打算用向量的模来表示影长1。依据影长的公式可能分析出各参数对影长的影响。我们拟以影长为纵坐标，时间为横坐标作图来绘制曲线。4.2模型的建立我们以地心为坐标原点O，设直杆所在经线与赤道的交点为B点，以OA所在直线在X轴，地球自转轴为y轴，z轴过圆心且与xoy平面垂直。直杆的影子变化与直杆长度、经度、纬度和太阳光线直射角度有关，而太阳直射角度又与日期有关，所以我们设定了杆长、直杆所在经度、直杆所在纬度、太阳赤纬、时角五个参数。为了便于求解，我们将直杆向量平移到圆心，此时直杆底端位于圆心，设直杆的顶端A坐标为。在坐标图中由简洁的几何关系可以得到A的坐标与四

8、个参数间的关系： （1）其中 （2）T是观测日期，它从2015年1月1日开头计时，每天+1。 （3）由四个参数的定义可以得到太阳光线向量为 （4）由式（4）得向量所在直线的斜率为，过点A且斜率为的直线即为从太阳放射经A照射到平面S的光线所在的直线，它的点斜式方程为 （5）设过原点且与OA垂直的平面为平面S，平面S的点法式方程表示为： （6）其中为平面S内已知点的坐标，即原点坐标；为平面的法向量即。带入上述坐标值，可以将式（6）写为 （7）直线与平面S的交点即为直杆影子的顶点,联立方程（5）和（7）即可得到顶点的坐标。求投影点坐标的Matlab程序和运行结果见附录。由运行结果我们得出投影点的坐标

9、为：直杆影子的长度即为向量的模， （8）4.3模型的求解（1）影长与参数间关系分析（ki为常数）我们以其它参数为常数，杆长为变量，对影长与杆长间的关系进行分析，得到可以得出：在其他条件一样的状况下，影长与杆长呈线性关系，影长随杆长的增加而增加。我们以其它参数为常数，纬度为变量，对影长和纬度间关系进行分析，近似得到可以看出影长与纬度近似呈二次函数关系，且存在一个最低点，随纬度的增加影长先减小后增加，最低点对应纬度为0.我们以其它参数为常数，日期为变量，对影长和日期间关系进行分析，近似得到可以看出影长与日期呈二次函数关系，有一个最低点，对应冬（夏）至日，随日期的增大，影长先减小后增大。 同理分析，

10、影长随时刻的增加呈先减小后增加的趋势。（2）坐标图的绘制本题要我们求2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线，我们由题目可以读出直杆长度、经度、纬度和日期，由模型一可以求出太阳赤纬。表1 问题一已知条件表符号数值直杆的长度2米直杆所在经度东经116度23分29秒直杆所在纬度北纬39度54分26秒日期2015年10月22日，即295天将上述已知条件带入公式（2）和（3），用Matlab软件求出和 并将坐标和参数的具体数值带入向量的模的公式（8），然后绘图，程序见附录1。我们得到下图：

11、图1 问题一影长-北京时间图由图得，在北京时间9:00-15:00间，影子长度呈先减小后增加的趋势，二者符合二次函数关系，图像存在一个最低点，此点对应北京时间12:00，由于直杆位于东八区，所以12:00也是最低点对应的当地时间。4.4模型的检验由于模型一是基于地球是一个规章的球形的假设来计算的，存在误差，为了消退误差，我们考虑到地球的实际外形，建立了太阳高度角模型2，来检验模型一。直杆、直杆影子和太阳光线可以构成一个直角三角形，令太阳高度角为，直杆长为，直杆影长为，则它们之间的关系为：太阳高度角与四个因素有关，分别是太阳赤纬、纬度、和时角有关，它们之间的关系为：太阳赤纬的计算公式为：其中，T

12、为观测日期，从每年1月1日开头计算，每天数值累计+1。本题要我们求2015年10月22日对应T=295，将T带入上式得：纬度的计算公式为时角的计算公式为 （9）其中为真太阳时；为北京时间；为北京时间与当地时间的时差；为天数，从1月2日开头计时；为修正参数，由于地球不是一个规章的球形，所以需要修正。是当地纬度，由于北京时间指东京的地方时间，有纬度差，而每经度变化，时间相差1小时。用上述模型求出9:00时的影长为5.67米，而模型一求出此时的影长为5.13米，偏差为-9.53%，偏差不大，所以模型一合理可信。5问题二模型的建立与求解5.1问题分析 考虑到问题二恰好有两个条件的确定性与问题一相反，我

13、们打算接受问题一的模型倒推经纬度，来确定地点。我们拟通过勾股定理对附件1中的影子的坐标进行处理，求21组直杆影子在对应北京时间的长度。以影长为纵坐标，北京时间为横坐标作图，可能得到一个最低点，这个点对应一个北京时间和当地时间（12:00）。我们预备依据问题一模型中的北京时间与当地时间的关系式，反推出当地的经度。对于纬度的确定，我们打算接受遍历的方法，将连续变化的纬度离散化，可能找出直杆的方差符合要求的纬度。经度和纬度求出后，排解出不合适的地点，正确地点就可以确定了。5.2模型的建立（1）影长的求法：用勾股定理求出21组直杆影子在对应北京时间的长度。（2）经度的求法：由问题一知影长与时间呈抛物线

14、，以影长为纵坐标，北京时间为横坐标作图，可以得到一个最低点，这个点对应一个北京时间，此点同时对应当地时间12:00。由于经度每差，时间相差一小时，我们得到北京时间与当地时间的关系为北京时间，为当地时间。由于已知，所以经度可求。（3）太阳高度角的求法：太阳高度角公式为依据时角公式上述各字母意义与问题一相同。将求出的数值带入上式，即可求出太阳高度角与纬度的关系式。（4）纬度的求法：我们用遍历法来确定纬度。以北纬为正，南纬为负，纬度区间为。将连续变化的纬度离散成每一纬度一个值进行遍历，确定出一个可能的纬度区间。每一个可能的纬度与经度结合可以确定出一个当地时间。由问题一可以得到直杆与其影子的关系为由于

15、对一根杆来说，它的高度肯定，所以我们求这21组直杆长度的方差，对应方差最小的纬度是最优解。为组直杆长度的平均值。进行方差检验求出可能区间。然后在可能区间内，细化步长，按上述方法再次求解，直到得出符合要求的结果为止。我们确定目标函数如下：综上所述，问题二的模型为 5.3模型的求解（1）经度的求解求解经度的Matlab程序见附录2，运行结果图如下：图2 问题二影长-时间图由上图得影长与时间呈二次函数关系，由系数可求最低点对应的时间为12:29，经度为110.23度。（2）纬度的求解编写Matlab程序及运行结果见附录2，由结果我们得到方差较小的，满足条件的可能区域为：表2 问题二纬度求解表纬度34

16、3536杆长1.296371848121611.294834288103711.29322806743802纬度373839杆长1.291553595161211.289811295328511.28800160666087纬度404142杆长1.286124982183111.284181888853551.28217280718565纬度434445杆长1.280098230862491.277958666344341.27575463247034由上表我们可以得到可能的纬度区间为北纬。我们再次运用遍历法，将此区间进一步细化，选择直杆方差最小的一组数据对应的纬度，确定此纬度为北纬40度。综

17、上所述，该图像拍摄的地点为东经110.23度，北纬40度，此地点对应内蒙古自治区鄂尔多斯市。6问题三模型的建立与求解6.1问题分析问题三在问题二的基础上增加了日期变量，我们考虑先建立模型求出日期，地点的求解与问题二全都。此问题就转化为了建立模型求日期问题。由于日期与赤纬角存在函数关系，我们打算通过求赤纬角来确定日期。我们拟引入新的变量方位角3，通过几何关系，可能得出方位角与赤纬角的关系。我们打算用太阳高度角和时角来表示方位角。我们预备引入影长与x轴的夹角来表示方位角，此夹角对日期的变化率等于方位角随日期的变化率，由此可建立一个等式。由于方位角又可以由太阳高度角和时角来表示，由方位角的两种表示方

18、法可以建立等式方程。通过这种方法可能求解出日期。6.2模型的建立日期T的求法：我们通过赤纬角求日期，赤纬角通过新变量方位角来确定。太阳方位角是指太阳光线在地平面上的投影与当地经线的夹角，可近似看作是直立在地面上的直线在阳光下的阴影与正南方的夹角s， 对于方位角的确定，我们通过定义影长与x轴的夹角来表示。太阳方位角表示为，i表示附件2第i组数据。由于和关于时间的变化率相等。我们得到方位角还可以由太阳高度角和时角来表示由问题一的公式（9）可以求出时角t。由于日期是连续的，不利于我们逐点分析，所以我们接受遍历法对日期以1天为步进步行遍历。方法类似问题二，这里不再赘述。综上所述，我们建立模型如下： 6

19、.3模型的求解（1）日期的求解求解日期的Matlab程序见附录3，由运行结果我们得到：附件2 数据对应的日期为172和355，分别对应6月21日和12月21日，即冬至和夏至日。附件3 对应的日期为172和355，分别对应6月21日和12月21日。（2）经度的求解图3 附件2影长-时间图由上图，我们得到：对于附件2,我们得到曲线的三个系数分别为0.000245，-0.0253,1.2725，由系数我们求出最低点对应的北京时间为15:16，经度差为49度，所以所求地点的位置位于东经度71度。图4 附件3影子-时间图由上图，我们得到：对于附件3,我们得到曲线的三个系数分别为0.0007411，0.0

20、108, 3.5227，由系数我们求出最低点对应的北京时间为12:47，经度差为11.75度，所以所求地点的位置位于东经度108.25度。（3）纬度的求解纬度的求解同问题二，附件2可能的纬度为表3 附件2纬度求解表纬度日期6月21日12月21日附件3可能的纬度为表4 附件3纬度求解表纬度日期6月21日12月21日（4）地点和日期的求解结果我们将纬度与经度结合，在地图中查找地点，由于存在某些位于海上的点，不符合我们的要求，所以我们将其舍弃。最终得到符合要求的地点。对于附件2，在6月21日，符合的地点为Jalal-Abad，对应东经71度，北纬42度；在12月21日，符合的地点为Karagandy

21、 Region，对应东经71度，北纬48度。对于附件3，在6月21日，符合的地点为Burytia,Sibirsky federating okruh，对应东经108度，北纬53度；在12月21日，符合的地点为Erdene,Tov，对应东经108度，南纬48度。7问题四模型的建立与求解7.1问题分析 此题给出的是直杆在太阳下的影子变化的视频截图，并且已通过某种方式估量出直杆的高度为2米，我们打算用Matlab来处理图像，应当可以依据杆长求出影长。当日期已知时，此题与问题二相像，我们拟接受问题二的作图法，结合时差与经度的关系式来求解经度。对于纬度的求解，由于在问题一的检验模型中我们知道了太阳高度角

22、与包含纬度在内的四个变量之间的总关系式，所以我们打算建立太阳高度角模型，预备求解出高度角与其他三个变量，应当可以通过总关系式得出纬度来。 当日期未知时，经度的求法不涉及日期，所以上一小问的求法照旧可以用。而纬度的求解时，由于上一小问用到了太阳赤纬，它随日期而变，所以我们需要接受其它方法求解赤纬，若求出赤纬，则此题与上一题变得相同，可以接受上一小问的方法求解。结合问题三，我们考虑接受遍历模型，对纬度进行遍历，查找合适的纬度，求出赤纬角。7.2模型的建立7.2.1日期已知时（1）经度的求法：由问题一知影长与北京时间呈抛物线，以影长为纵坐标，北京时间为横坐标作图，可以得到一个最低点，这个点对应一个北

23、京时间，此点同时对应当地时间12:00。由北京时间与当地时间的关系式 我们可以求出经度。（2）纬度的求法：由问题一、二，我们得到如下关系：太阳高度角与太阳赤纬、纬度、时角的关系为： 时角的计算公式为：太阳赤纬的计算公式为：太阳高度角与直杆长、影长的关系：我们通过图像处理，依据上述公式可以得出、和的大小来，由太阳高度角与太阳赤纬、纬度、时角的关系式，我们可以求出纬度的大小。结合经纬度，可求出地点。7.2.2日期未知时经度的求解方法同7.2.1。纬度的求解：参考上一小问的求解过程，我们使用如下关系式： 将纬度以为步进步行遍历，遍历范围为，由每一个纬度，我们可以对应求出一个太阳赤纬角，计算方差值我们

24、保留方差值小的赤纬，由此确定出可能的区域。再对此区域进行遍历，将步长细化十倍，直到方差足够小为止，我们就可以得到一个合适的纬度区间。综合分析纬度与经度，我们可以得到合适的地点。7.3模型的求解由附件4求得15张图片的影长为表5 附件4影长表序号123影长2.389743589743592.329914529914532.29914529914530序号456影长2.264957264957272.222222222222222.18803418803419序号789影长2.140170940170942.111111111111112.07179487179487序号101112影长2.037

25、606837606842.003418803418801.95726495726496序号131415影长1.935042735042741.890598290598291.85982905982906结合上表数据，我们以影长为纵坐标，北京时间为横坐标作图图5 问题四影长-时间图我们由此图的最低点求得经度为东经102.165度。然后由附录4求解纬度的Matlab程序，求得纬度为北纬58.1度。结合地图我们得到此处为Irkutsk Oblast,Sibirsky federalny okruh.8模型的评价8.1模型的优点1.针对几何模型，通过建立坐标系，合理的将影子长度问题转化为了向量模的问题

26、，思路清楚，便于理解，运用了点、线、面之间的关系，很好的解决了影长的求算问题。2.针对太阳高度角模型，它基于地球的真实外形来考虑，引入了修正值，可信度格外高。3.针对遍历模型，它将连续变化的纬度离散化，逐点求解，极大的缩小了纬度区间，再通过细化步长，进一步缩小纬度区间，给出了合理的纬度。8.2模型的缺点由于附件所给的数据太少，使我们的作图过程可用信息太少，图的精确度达不到极高的水平。8.3模型的推广 本文的太阳高度角模型可以推广到建筑领域，对于建筑的采光问题格外有效。同时，我们可以将直杆换成其他任何难以求高度的物体，利用本文模型来求算其高度。参考文献1 高等数学（下）M，北京：高等教育出版社，

27、38-39，2007.62 胡毅华，杨旭龙，刘媛萍，太阳影子定位模型的构建J, 洛阳师范学院学报，Vol.34, No.11:13-18，2015.113 郑鹏飞，基于影子轨迹线反求采光效果的技术争辩J，华东理工高校学报: 自然学科版，Vol.36，No .3 : 458-463，2010附录附录11.求投影点坐标的Matlab程序（由于Matlab程序不能识别和，我们用w和c分别代替）syms x y z w r c l t;eq1=(x-l*cos(w)*cos(t)*tan(c)-(y-l*sin(w);eq2=z-l*cos(w)*sin(t);eq3=l*x*cos(w)*cos(t

28、)+l*y*sin(w)+l*z*cos(w)*sin(t);x y z=solve(eq1,eq2,eq3,x,y,z)2.运行结果x=-(l*sin(w)2+l*cos(w)2*sin(t)2-l*tan(c)*cos(t)*cos(w)*sin(w)/(cos(t)*cos(w)+tan(c)*sin(w)y=-(cos(w)*(l*tan(c)*cos(t)2*cos(w)-l*cos(t)*sin(w)+l*tan(c)*cos(w)*sin(t)2)/(cos(t)*cos(w) + tan(c)*sin(w)z =l*cos(w)*sin(t)3.画图像的Matlab程序w=40

29、/180;l=3;t=-0.8:0.01:0.8;c=0;for i=1:161x(i)=-(l*sin(w)2+l*cos(w)2*sin(t(i)2-l*tan(c)*cos(t(i)*cos(w)*sin(w)/(cos(t(i)*cos(w) + tan(c)*sin(w);y(i)=-(cos(w)*(l*tan(c)*cos(t(i)2*cos(w)-l*cos(t(i)*sin(w)+l*tan(c)*cos(w)*sin(t(i)2)/(cos(t(i)*cos(w) + tan(c)*sin(w);z(i)=l*cos(w)*sin(t(i);H(i)=(x(i)2+y(i)

30、2+z(i)2)0.5;endT=9:0.0375:15;plot(T,H,'r-')ylabel('高度/米');xlabel('时间')附录21.求经度的Matlab程序for i=1:21 y(i)=(S1(i,1)2+S1(i,2)2)0.5;endx=1:21;A=polyfit(x,y,2);y1=polyval(A,x);2.求纬度的Matlab程序B=0.048;for c=-90:90 for i=1:21 t(i)=0.262*(2.7+0.05*i); end for i=1:21 A1=sin(B)*sin(c/180)+cos(B)*cos(c/180)*cos(t(i); af=asin(A1); H(i)=y(i)*tan(af); end k=c+91; Hp(k)=sum(H)/21; for i=1:2