基于太阳影子定位的研究数学建模论文

1、基于太阳影子定位的研究摘 要本文研究的是在经纬度、日期、时刻、杆长、影长等变量之间建立起量化关系，构建出模型，能够达到在已知一些量的情况下求解出其余的量。针对问题一：我们先将各个变量之间有数学表达式建立出关系,然后通过球面三角形的运算公式，借助高度角、赤纬角、时角等相关概念求解出实际问题。针对问题二：由于无法确定题目给出的数据是如何建立坐标系的，所有就考虑到可以通过建立函数关系和拟合曲线的方式来解决问题。首先，根据附件坐标，将与其一一对应的影长计算出来，并将这些数据进行拟合后可以发现影子长度与时间所构成的关系式完全满足二次函数，拟合的相关系数为1。所以由二次函数图象的特殊性，图象的最低点就可以

2、直接取到正午影长。再根据时间差与经度差的关系，就可以直接计算出经度。之后又由于构建了杆长和当地纬度的关系式，进行等量取值后再用matlab软件进行循环处理，每循环一次就可得出一个正午影长，将此时循环得到的影长与之前拟合的影长数进行取差比较，绝对值最小影长对应的纬度就为最后所求的结果。整个思路主要是避免了由于坐标系不确定而带来的误差。针对问题三：本题求解经度的方式还是同题二大致一样，都是通过拟合影子长度，获得函数曲线，再由时间差求解经度差。对于纬度和日期求解采用的是间隔取值的方法，运用matlab软件进行循环，解出不同时间的影子长度，在于附件二、三中计算的影长做比较，计算方差，由最小方差来确定最

3、合适的解。针对问题四：本问题中只给了一段视频，所以要通过2d图的坐标转换成3d空间坐标。先假设相机光学中心所在的平面坐标系以及摄像机平面与直杆所在平面的法向量，通过视频间隔时间段截图，获取每幅图上影子端点、直杆与地面的交点以及直杆顶点的二维坐标，通过转换可得到影子端点在像平面的坐标与真实的影子端点的坐标之间的关系，将影子端点在像平面的坐标转换成真实坐标，然后再对数据进行处理得出要求的地点 。最后，我们对模型进行了优缺点分析，以及下一步要进行的工作.关键词 数据拟合；循环计算；matlab；数图分析；一、问题重述a题 太阳影子定位如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子

4、定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。1.找出问题中的变量，分析影子长度关于各个参数的变化规律，从而建立影子长度变化的数学模型，并应用建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。2.根据附件1中给的顶点坐标数据，通过建立数学模型确定直杆所处的地点，给出若干个可能。3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日

5、期。4附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。建立确定视频拍摄地点的数学模型，并应用模型给出若干个可能的拍摄地点。最后，如果拍摄日期未知的情况下，尝试根据视频确定出拍摄地点与日期。二、问题分析针对问题一：首先要先找到所有的变量，再确定了各个变量之间的关系之后，建立出数学表达式。最后把已知的量带入，计算求解出其余的量。针对问题二：本题的需要通过将图像与数据相结合的分析方式，才能在坐标系方向不确定的情况下求出经纬度。主要是通过坐标来求出影长，进行曲线拟合后，通过图像的特殊规律来作为切入点。针对问题三：将附件二和附件三的数据进行处理，运用matlab数据处

6、理软件将处理后的数据拟合出该天的影长变化曲线，以最低点作为切入点，求出所求地点的经度，利用纬度和杆长的关系将杆长用纬度替换，减少变量，以日期为第一变量，纬度为第二变量，分别在北京时间12:41到13:41时及北京时间13:09到14:09时进行循环操作，每次运行都可输出一组影长，再把得出的所有影长与原始数据中影长进行最小二乘法回归拟合，在这两段时间内若对应的拟合结果在一定范围内符合即为可能结果，可能结果对应的日期和纬度就是要求的最优解。针对问题四：若要对问题更进一步的求解，则需对视频进行处理，进而得出影子端点的真实坐标，再根据得出的真实坐标对目标进行求解。三、模型假设假设1：地球在一日当中只考

7、虑自转（即不考虑公转的影响）；假设2：相机光学中心到杆的直线垂直距离保持不变；假设3：忽略由于风速、太阳光等因素对像素的影响。四、符号说明 ：表示不同测量位置的经度；：表示不同测量位置的纬度； ：表示不同的测量时刻； ：表示不同的测量日期；：表示不同的直杆长度；：表示直杆的影子长度； : 表示太阳高度角； ：表示赤纬角； : 表示太阳时角； ：表示真太阳时； ： 表示测量的天数；：表示正午影长；五、模型的建立与求解5.1、问题一的分析与求解：5.1.1 问题一的概述生活中常常会遇到这样的情况，傍晚夕阳西下时，人身后的影子会很长，但是在中午太阳最炙热的时候，影子却明显缩短了。即使连续很多天的正午

8、时分站在同一个位置，也会发现身后的影子也会有不同的长度，更别说站在不同的位置了。由此可以总结规律。在测量影子长度时，经纬度、测量日期、时刻和位置都会直接影响影子长度。容易被忽略的是，不管是人的身高还是杆的长度，同样也会对影子的长度造成影响。通过影响光线与地平面的之间的夹角从而间接影响影子的长度，这个角称为高度角。确定了各变量之间的关系之后，就想设法将各量之间的关系量化，由表达式构建出模型。5.1.2 问题一的模型的建立与求解5.1.2.1 模型一：影子长度变化的模型题目要求是解出3米高的直杆的影子长度变化曲线，其关键在于要先确定出9:00-15:00这个时间段中高度角的变化范围。通过查阅文献1

9、得知在求解高度角的问题上，通常运用天文学中常见的球面三角形来求解。根据球面三角形的余弦公式，将其中的相关的各角度找关系求出后，最后联立方程，就可推出太阳高度角h的计算公式。再将杆长除以高度角的正切值，就得出影子的长度。针对模型的具体分析如下：太阳光照射直杆的图像可以简化为图（1）图（1）太阳照射直杆后投射出影子，其中太阳光线和地平面的夹角h就是高度角。由图可以看出，求解影子长度的关键是在于把上图构成的直角三角形中的高度角求出。因为杆长是已知条件，所以最后可直接通过杆长除以高度角的正切值，求出影长。通过查阅文献1可以得出，高度角的计算式需要通过运用天文学中常用的球面三角形中边的余弦公式推导求出：

10、球面三角形的任意一边的余弦等于其他两边预选的乘积加上这两边的正弦及其夹角的余弦的连乘积。（球面三角形如图（2）图（2）公式表达式为：（以a边为例） 【1】 另外两边直接通过等效代换得出。假设地球近似为一个球体，所以相对于围绕它自东向西的运行的太阳来说，地球即使存在自转也可视为为静止不动的参照物，太阳则就是在绕地球作近似圆周的运动。图（3）各角度参数图2对于所求的问题，可将直杆竖立的原点和影长移动的端点分别看作静止的地球和运动的太阳。如图（3）所示，定义与地球赤道位于同一平面的球面为赤道面；与地球地平线圈位于同一球面上的球面圈为地平圈；与地球上经度圈位于同一平面上的球面圈为时圈（通过南北极的球面

11、圈称为经度圈）3图中的太阳高度角是指某地太阳光线与通过该地与地心相连的地表切线的夹角，也就是图中太阳位置点与地球位置点连线与其在地平圈上的投影的夹角。当太阳高度角为90°时，此时太阳辐射强度最大，影长较短；当太阳斜射地面时，太阳辐射强度就小影长较长。因此这个角度是个变化的量。方位角是从过地球位置点与太阳投影点的连线从指北方向线起，依顺时针方向到赤道圈的水平夹角。赤纬角是地球赤道平面与太阳和地球中心的连线之间的夹角。通常，太阳所处位置s的表示是通过用坐标来说明。一般情况下可建立以地平坐标为横轴，赤道坐标为纵轴的坐标系，再由相关的太阳高度角,方位角，赤纬角和时角来表示。太阳高度角： si

12、n= sinsin+cos 由球面三角形公式结合图（2）可推导出太阳高度角的求解公式，其中涉及到的量分别是赤纬角、纬度和时角。下面就需要继续求解这些量。赤纬角： （为日数，自每年1月1日开始计算 ）太阳时角： 真太阳时： 将式带入式中，即可得出sin的值，进一步即可得出的值，综合上述分析，建立可建立出具体模型，如下：影子长度： 其中，代表的是直杆长度，表示的是太阳高度角。根据式，可以得出 式中代表的是纬度；代表的是赤纬角，可由式求出； 代表的是时角，可由式求出。变量之间的关系都由已经通过数学表达式量化，可以通过已知其中某几个量求解其余的量。5.1.3 问题一的模型求解5.1.3.1 模型一的求

13、解由题目可得：已知日期为2015年10月22日，由此可根据式得到天数为265天，纬角的值为-0.6054；再根据经度为东经116.39139°，可根据式和得到太阳时角与时间的关系；又由于纬度为北纬39.90722°，将所有的已知量带入模型，进行量化计算后，得出影长在9:0015:00的变化曲线图（如图（4）所示）。图（4）影子长度时间变化曲线通过matlab软件进行计算后，得出了影子长度时间的变化曲线，从图（4）中可以看出，2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子在9:00是此

14、时间段中最长的，达到7.75m；随着时间的推移，影子长度在逐渐缩短，到12:10左右是此时间段中最短的，有3.84m。之后，影子的长度又继续增长。总的来说，整个过程呈现先减后增的趋势。5.2问题二的分析与求解5.2.1 问题二的概述由于附件一的坐标系无法确定，若与建立出来的模型坐标系不一致，数据就不可用。所以无法直接利用附件中的数据进行处理，但可以间接计算出附件中所列时刻的影长，进行数据拟合后，从图像得到信息，再根据相关计算分别求出经纬度，最后确定地点。5.2.2 模型二的建立即使附件一已经直接给出了一定的坐标，但是由于考虑到坐标系建立的不确定性，就采取通过影长对经纬度分别求解的方法进行操作。

15、图（5）附件一拟合结果图图（6）附件一拟合曲线图图（7）影子长度时间二次函数关系图首先用附件一中已知的坐标计算出相应的影子长度，再将整个时间段和影子长度拟合，得到一个图像和对应的函数表达式，从图像来看，是近似于二次的一个函数。但根据拟合后得出的相关系数为1，证明这个函数表达式是完全符合于二次函数的。又可根据图（7）中二次函数图像的特性，在横轴为时间，纵轴为影长的情况下，模拟出一天的影长变化曲线。以附件中所给数据拟合后，曲线最低点的纵坐标代表的是当地正午影长为0.4918m，横坐标为12.5856。在求解经度的过程中，不打算采用穷举法和最小二乘法回归拟合，而是根据经度共有360°，1天

16、是24小时，360/24=15°，所以经度每15°，时差就是1小时，直接通过这样的关系，利用时间差求出精度差。将式变形就可得计算当地经度的公式 （北京时间为12时）此时杆顶和影子的顶端依然构成直角三角形，就可利用高度角求解，拟合曲线的最低点的纵坐标表示的是正午时的影长。 (表示的是正午太阳高度角、表示太阳直射点的纬度和当地纬度的差) 可将式带入式，整理后就可建立出杆长与当地纬度的关系式 每一个影子的端点都可以求出一个影长，所以求解影长的关键就在于确定了和的值。由和可得到影长的计算公式 又根据方位角公式 4 可得出方位角和、的关系。 ； 所以每一个纬度对应一个点，即每一个纬度

17、对应一个影长。现取纬度的范围为，纬度的步长设为0.001（rad），计算出每个纬度对应的各自影长，将计算得到的影长分别与拟合得出的正午影长取差的绝对值，绝对值最小的结果对应的纬度即为所求纬度。由于事先已把经度求出，现只需要在该经度上找符合要求的纬度点就可以。5.2.3 模型二的求解分别对数据中的进行计算处理，可得下表。表一 附件一影子长度-北京时间关系表把表一中时间和影长利用matlab软件拟合出该天影长随时间的变化曲线如下图所示：图（8）附件一中影长随时间的变化曲线通过matlab软件进行影长和纬度关系的处理，程序代码详见【附录2】：s,t=min(chazhi); % t为最优解的返回位置

18、根据t的结果可计算出所求地点的纬度=（t×0.001-1.2217305）×180°÷=（2001×0.001-1.2217305）×180°÷=44.64885345°根据式可得出所求地点的经度= =128.784°即所求地方在北纬44.64885345°，东经128.784°处。5.3问题三的模型建立与求解：5.3.1 问题三的概述现已知附件二中北京时间12:41到13:41时的21个影子顶点坐标（，）以及附件三中北京时间13:09到14:09时的21个影子顶点坐标（，）

19、，求可能的地点与日期。5.3.2 模型三的建立与求解首先利用式对已给出的数据算对应的影长，再用matlab软件拟合影长和北京时间的关系曲线，曲线上的最低点对应的横坐标表示当地的正午，纵坐标表示最短影长，利用式可求得当地的经度。利用可得出纬度和杆长的关系，即可将杆长用纬度代替。再根据式可求得影长。设日期和纬度为两个变量，可从1取到365；由于测量地点在陆地，所以排除海域的纬度，即可取20°到74°。以为第一变量，为第二变量，在北京时间为12:41到13:41时代入式中，在此范围内进行双重循环，每次循环可得21个与附件二相对应的影长，同样在北京时间13:09到14:09时代入式

20、中，在此范围内进行双重循环，每次循环可得21个与附件二相对应的影长对每次的循环结果与相对应附件中的21个影长进行最小二乘法拟合，在一定范围内筛选出合适的一组或多组影长，从而即可得到最有可能的纬度和日期。表二 附件二影子长度-北京时间关系表把上表中时间和影长在matlab上进行拟合，可得到一天的影长变化曲线如下图所示：图（9）附件二中影长随时间的变化曲线表三 附件三影子长度-北京时间关系表把上表中时间和影长在matlab上进行拟合，可得到一天的影长变化曲线如下图所示：图（10）附件三中影长随时间的变化曲线通过matlab软件对数据进行处理后得出的纬度，代码源程序详见【附录3】：=21.12

21、76; =92.03° 日期为2015.5.12=40.79° =108.94° 日期为2015.3.2或2015.10.95.4问题四的模型建立与求解：5.4.1 问题四的概述本题求地点的过程可转化为建立像平面坐标系和真实坐标系之间存在的关系，将像平面坐标转化为真实坐标系中的坐标，进而再进行问题二相似的求解过程。所以先在附件四的视频中截取20幅图片，找到每幅图片上影子端点的像平面坐标，并记录下来。5.4.2 模型的建立与求解 设20幅图片的像平面坐标分别为，.,相机像平面与直杆所成的夹角为，相机光学中心到地平面的法向量为n=，法向量为n=与真实坐标平面的交点坐标

22、为m（，），因为点m（，）在球面上，所以满足球面方程。由上面的公式可求得像平面上的坐标投影到地面坐标后应为,.。利用像平面上的坐标和地面上的坐标分别做商比较，得出像平面和地面坐标系之间的关系，即可将像平面坐标系中的坐标转化为地面坐标系中的坐标，再利用即可求出拍摄地点的经纬度。六、优缺点分析6.1 优点我们通过对问题的完整分析考虑，能够再数据条件不完整的前提下，采用建立函数模型、拟合图像的方式，来分析问题。能够有效的避免在坐标系不确定的情况下带来的误差，可以让我们以性能快速解决问题。6.2 缺点未考虑地球一天中的自转影响，导致计算的结果会有误差。七、参考文献1 赖月喜，中学政史地，福建宁化：河南

23、大学出版社，2007年4月2庄肃，太阳能，知网，19933西安冶金建筑学院、清华大学,建筑物理m, 北京：中国建筑工业出版社，1979:101-1104吉林省建筑设计院，建筑日照设计m，北京：中国建筑工业出版社，1979：17280八、附录本次使用软件为matlab r2014a附录1：clear;d=286;jingdu=116.3913;w=0.6965;c=23.45*sin(2*pi*(d+284)/365)/180*pi;t1=9:0.2:15;t=(15*(t1-(120-jingdu)/15-12)/180*pi;sinh=sin(w)*sin(c)+cos(w)*cos(c)*

24、cos(t);tanh=tan(asin(sinh);gc_juzhen=3*ones(1,31);yc=gc_juzhen./tanh;plot(t1,yc);grid on附录2：clear;kesi=23.45*sin(2*pi*(108+284)/365)/180*pi;a=85.35;b=-89.51;c=23.96;zz=-b/(2*a);x=0.375:0.001:0.625;y=a*x.2+b*x+c; plot(x,y); grid on; %画图 jingdu=120+360*(zz-0.5); zwyc=a*zz2+b*zz+c; sink=sin(kesi); cosk

25、=cos(kesi); beitime=12; n=(120-jingdu)/15+beitime; cosou=cos(15*(n-12)/180*pi);sinou=sin(15*(n-12)/180*pi);m=;n=;for w=-1.2217305:0.001:1.2217305 h=zwyc/(tan(pi/2-kesi-w); h=asin(sin(w)*sink+cos(w)*cosk*cosou); a=asin(cosk*sinou/cos(h); x0=(h*cot(h)/(1+(tan(a)2)(1/2); y0=tan(a)*x0; m=m,x0; n=n,y0; e

26、ndl=(m.2+n.2).(1/2);zwyc_juzhen=0.4918*ones(1,2444);chazhi=l-zwyc_juzhen;for i=1:2444 if(chazhi(1,i)<0) chazhi(1,i)=-chazhi(1,i); endends,t=min(chazhi);附录3：clear；for d2=1:365    d2    for w=0.3490659 :0.01:1.308997     

27、60;   ll=;        for beitime=12+41/60:0.05:13+41/60;                       n2=(120-jingdu2)/15+beitime;   

28、60;        cosou2=cos(15*(n2-12)/180*pi);            sinou2=sin(15*(n2-12)/180*pi);            cwj2=(23.45*sin(2*pi*(284+d2)/365)/180*pi;

29、60;            h=zwyc2/(tan(pi/2-cwj2-w);             h=asin(sin(w)*sin(cwj2)+cos(w)*cos(cwj2)*cosou2);             a=real(sin(cos(cwj2)*sinou2/cos(h);             x0=(h*cot(h)/(1+(tan(a)2)(1/2);&#