A201525003029_杨俊龙_苏婷_李存课件

1、太阳影子定位摘要本文对太阳影子定位问题进行了深入探讨，建立了太阳影子和时间以及变位参数之间关系的数学模型，应用了matlab进行求解。针对问题一：已知时间与地点，求解影长与时间的关系。通过影长与太阳高度角、太阳方位角的关系，建立参数模型，运用matlab画出影子长度变化曲线，影子长度与太阳高度角、太阳方位角、时间之间的关系如表1。表1 2015年10月22日9:00-15:00天安门广场直杆（3米）的太阳影子长度当地时间9:0010:0011:0012:0013:0014:0015:00太阳高度角23.319531.00136.155837.99136.155831.00123.32太阳方位角

2、A40.021139.61939.375939.29439.375939.619240.021影子高度L13.8417573.8413643.8411283.8410493.8411283.8413643.841757针对问题二，已知日期n与时间t以及影子的观测坐标P，经度、纬度、杆长h未知。分析太阳高度角、太阳方位角A与参数h,的关系，求解理论坐标Q。根据运筹学理论，以各PQ2之和最小为目标函数，以h,的取值范围为约束条件，建立非线性规划模型。应用最小二乘参数估计法思想，运用穷举法和放缩法，利用matlab求解，得出,的值。针对问题三，已知t与P，在问题二的基础上增加一个未知参数日期序号n，

3、分析、A与参数h,n的关系，求解理论坐标Q。以各PQ2之和最小为目标函数，以h,n的取值范围为约束条件，建立非线性规划模型。应用最小二乘参数估计法，利用matlab求解，得出,的值。针对问题四，利用视频帧实现经纬度估计的方法，根据未经过校准的影子的位置，得出所拍摄图像的经纬度信息，利用matlab求解，得到相应的经纬度。关键词：影子定位，参数方程，非线性规划，最小二乘参数估计法一、 问题重述太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。问题一：（1）建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律。（2）应用建立的模型画出2015年10

4、月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。问题二：根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。问题三：根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。问题四：（1）附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数

5、学模型，并应用你们的模型给出若干个可能的拍摄地点。（2）如果拍摄日期未知，如何根据视频确定出拍摄地点与日期？二、 问题分析2.1问题一的分析（1）建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律。从太阳影子长度随着时间的变化而改变的角度出发，利用不同时刻的太阳高度角及太阳方位角来计算太阳影子长度。建立太阳影子顶点坐标的参数模型。（2）应用建立的模型，代入数据，运用matlab进行拟合，得到2015年10月22日北京时间9:00-15:00之间天安门广场3米高的直杆的太阳影子长度的变化曲线。2.2问题二的分析已知直杆顶点坐标投影和时间（即已知日期序号和时间t），直杆高度固定，建立数

6、学模型确定直杆所处的地点。要根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。要将太阳高度角和太阳方位角分别进行讨论，考虑到实际观测值P所在的坐标系和正北方向存在夹角的话会使表达式过于复杂，于是假设实际观测值P所在的坐标系和正北方向存在夹角为零。根据太阳影子长度与时间、直杆长度、经纬度之间的关系，建立S关于参数h,的非线性规划模型，运用最小二乘参数估计法求解，利用matlab求解找到最佳的h,值，从而确定附件1直杆所在地点。2.3问题三的分析已知直杆顶点坐标投影和时间（年月日时分秒，即已知日期序号和时间t），直杆高度固定，建立数学模型确定直杆所处的地点。在问题二

7、的基础上根据太阳影子长度与时间、直杆长度、经纬度、日期之间的关系，建立 S 关于参数h,n的最小二乘拟合模型, 运用最小二乘参数估计法求解，利用matlab求解找到最佳的h,n值，从而确定附件2、附件3中直杆所在地点以及日期。2.4问题四的分析（1）根据高度为2米直杆在太阳下的影子变化的视频，确定视频拍摄地点的数学模型。利用视频帧实现经纬度估计的方法，根据未经过校准的影子的位置，得出所拍摄图像的经纬度信息5，建立相应的灰色模型，利用matlab求解，得到相应的经纬度，得到拍摄地点。（2）若拍摄日期未知，根据视频确定出拍摄地点与日期。三、问题假设（1）假设不考虑大气折射率；（2）假设不考虑海拔、

8、大气压强、气候等因素；（3）假设直杆是细杆，不考虑杆的粗细对影长的影响；（4）假设所给数据真实可靠。四、符号定义和说明h杆高L1影长t北京时间ts太阳时太阳高度角A太阳方位角当地纬度当地经度当时的太阳赤纬w时角n日期序号1n365五、模型的建立5.1问题一的模型建立参数模型的建立思路：太阳方位角A太阳高度角时差ET经度赤纬角纬度时角w太阳时ts时间t日期n 5.1.1太阳高度角（1）计算时角w1 w=15（ts-12） （1）其中ts为太阳时2 ts=t-4120-+ET （2）其中ET为时角 ET=229.2（0.000075+0.001868cos-0.032077sin-0.014615

9、cos2-0.04089sin2 （3） =360n-13651n365 （4）（2）计算赤纬角1 =23.45sin2(284+n)365 度 （5）（3）确定太阳高度角3太阳高度角是太阳相对于地平线的高度角，这是以太阳视盘面的几何中心和理想地平线所夹的角度。太阳高度角可以使用下面的算式，经由计算得到很好的近似值： sin=sinsin+coscoscosw （6）5.1.2太阳方位角 4太阳方位角是太阳在方位上的角度，它通常被定义为从北方沿着地平线顺时针量度的角。它可以利用下面的公式，经由计算得到良好的近似值，使用正弦函数： sinA=-sinwcoscos （7）使用余弦函数：当时角为负

10、值时 (上午)，方位角的角度小于180度，时角为正值时 (下午)，方位角应该大于180度，即要取补角的值。 cosA=sincos-coswcossincos （8） cosA=sin-sinsincoscos （9）5.1.3建立参数方程定义图1中y轴方向为正北方向，影子顶点坐标为x，yx，yL1hz y x A图1 “立竿见影”坐标系由影子长度与太阳高度角的关系 L1=htan （10）可得出参数方程： x=sinAL1y=cosAL1 t为参数 （11）根据建立的参数方程，运用matlab拟合出影子长度与时间的变化曲线。5.2问题二的模型建立假设太阳影子顶点坐标实际观测值为PPxi,Py

11、i,太阳影子顶点坐标理论值为Q(QxiQyi),假设P所在的坐标系和正北方向的夹角为0，则有 =hsh, A=hah, Qxi=L1sinAQyi=L1cosA (12) L1=htan (13)由上述三个式子得到： Qxi=htansinA1Qyi=htancosA1 (14) 化解得到： Qxi=-sinwcoshsin (15) Qyi=sincos-coswcossinsinhQyi=sin-sinsinsincosh （16）注：当时角为负值时 (上午)，方位角的角度小于180度，时角为正值时 (下午)，方位角应该大于180度，即要取补角的值。通过建立理论坐标QQxi，Qyi与观测坐

12、标PPxi，Pyi的距离的平方和S最小的非线性规模型进行求解。即： minS=i=121（Pxi-Qxi）2+（Qxi-Qyi）2 (17)s.t.0360-9090 2h3 5.3问题三的模型建立在第二问的基础上增加一个未知参数n ,建立非线性规划模型： min S=i=121Qxi-Pxi2+Qyi-Pyi2 （18）s.t.0360-90902h31n3655.4问题四的模型建立 要根据一根高度为2米的直杆在太阳下的影子变化的视频，建立确定视频拍摄地点。就是要利用视频帧实现经纬度估计的方法，能够根据未经过校准的影子的位置，得出所拍摄图像的经纬度信息5。具体步骤5如下：（1） 获取视频中含

13、有至少两个影子轨迹的视频帧；（2） 对每一帧图像，检测影子轨迹；（3） 找出灭点，计算灭点；（4） 由灭点拟合出地平线；（5） 在图像沿着有垂直关系的物体画出两条直线，根据这两条直线分别与地平线的交点坐标，这两个坐标就是两个相互垂直的灭点vx,vy的坐标。（6） 根据公式f=v0vx3vy2+u0vx1vy3+v0vx2vy3+u0vx3vy1-vx1vy1-vx2vy2vx3vy3-u02-v02计算出焦距f,式中，u0,v0为图像中心点的坐标，根据公式g=10-u001-v0-u0-v0u02+v02+f2计算出g（7） 在地平线上找到另一对满足vxTgvy=0限制的点vp,vq；（8）

14、实现从图像到经过度量纠正的世界坐标的转换；（9） 根据X=Hx还原出经过度量纠正的世界坐标中的点，其中X是经过度量纠正的世界坐标中的点，x是图像坐标中的点；（10） 拟合影子轨迹；（11） 利用日晷原理，计算纬度；（12） 运用极值法求出经度。（13） 根据得到的经纬度给出可能的视频拍摄地点。六、模型的求解6.1问题一的模型的求解曲线拟合根据参数方程模型，代入数据日期序号n=295，当地东经=116.8667，北纬=39.9072，杆高h=3米。以y轴为正北方，求解直杆的顶点坐标（x，y），得到2015年10月22日北京时间9:00-15:00之间天安门广场3米高的直杆的太阳影子的运动轨迹如图

15、2，影子长度随时间和方位角变化的曲线如图3，影子长度随方位角的变化曲线如图4，影子长度随时间的曲线如图5，不同时刻对应的影长如表1。图2 直杆影子的运动轨迹（附录1）图3影子长度随时间和方位角变化的曲线（附录3）图4 影子长度随方位角的变化曲线（附录3）图5 直杆影子随时间的变化曲线图（附录2）表1 2015年10月22日9:00-15:00天安门广场直杆（3米）的太阳影子长度（附录4）当地时间9:0010:0011:0012:0013:0014:0015:00太阳高度角23.319531.00136.155837.99136.155831.00123.32太阳方位角A40.021139.61

16、939.375939.29439.375939.619240.021影子高度L13.8417573.8413643.8411283.8410493.8411283.8413643.841757观察图形，分析数据：太阳高度越小，影子越长。中午12:00太阳高度最大，影子最短。早晚太阳高度最小，中午最大，早上9：00到中午12：00太阳高度慢慢变大，影子变短，中午12：00到下午15：00太阳高度慢慢变小，影子变长。此数据与日常生活中实际现象相符合，可见此模型是合理的。6.2问题二的模型求解最小二乘法用最小二乘参数估计法确定参数。最小二乘参数估计法基本思想：根据Sh,t的关系表达式求出太阳影子的理

17、论坐标，通过计算理论坐标与观测坐标的距离，建立距离平方和最小模型。通过对,h进行等间距的穷举最终求得理论坐标QQxi，Qyi与观测坐标PPxi，Pyi的距离的平方和S，当S取得最小值时，此时，h的值即为所求的最佳值。即求解如下最小二乘拟合模型。算法描述：输入：时间t、日期序号n输出：经度、纬度、杆高h的值Step1: 建立以y轴为正北方向的空间直角坐标系，以地面为水平面，直杆垂直于地面，建立理论坐标QQxi，Qyi的参数方程，根据太阳高度角、太阳方位角与影长的关系，求出理论值QQxi，Qyi关于，h的关系式。Step2:时间t1，t2，t3t21，日期序号n=108，等间距穷举（附录5）经度0

18、,360，纬度-90,90，杆高h2，3，根据理论值QQxi，Qyi关于，h的关系式，求得理论坐标QQxi，Qyi。Step3:根据求出的理论坐标QQxi，Qyi与附录1中的观测坐标 PiPxi，Pyi （i0,21）,求出距离平方和S=i=121（Pxi-Qxi）2+（Qxi-Qyi）2令S0初始值为，SS0,令S0=S，循环找出S0的最小值，最小值S=S0。当S取得最小值时，求得，h的值。算法结束运用穷举法与放缩法，用matlab求解得出最优解，如表2。取值范围0：30：360-90：30：902：0.1：3Step_1所得值120032.9145取值范围90：5：120-30:5:302

19、：0.1：:3Step_2所得值120-1032.6843取值范围115:0.1:125-15:0.1：-52：0.1:3Step_3所得值118.8-9.130.1086取值范围118.7:0.01:118.9-9.2:0.01：-92：0.1：:3Step_4所得值118.83-92.80.0007表2穷举法的缩放过程（附录6）得到，h可能值为=118.83，=9，h=2.8m。6.3问题三的模型求解在问题二的基础上，增加一个未知参数n，求解如下非线性规划模型：minSh,n=i=121（Pxi-Qxi）2+（Qxi-Qyi）2在matlab中根据n的取值范围1，365作为限制条件，引入一

20、个for循环，同理运用穷举法和放缩法进行求解。七、模型的评价7.1模型的优点（1）通过建立非线性规划模型,在利用未知参数取值范围的条件下，运用matlab编程得出了影子顶点坐标理论值与实际值的距离平方之和，最小二乘法求解最小值，该方法结果合理可信。（2）运用matlab绘制出图形，用数形结合的方法来进行分析，使模型思路更加清晰，更有说服力。（3）建立参数方程，可直接求解与逆向求解，模型具有普遍性。（4）通过对参数范围的缩小，减小了大量的运算，同时提高了精度，使目标函数的最小值达到一个很小的水平，可以认为实际值和理论值几乎相等S00.0007.7.2模型的缺点（1）由于未考虑实际观测值P所在的坐

21、标系和正北方向存在的夹角，导致问题二的求解出的值存在一定的误差。（2）由于运用穷举法，未知参数过多，等间距穷举参数有限，运算量过大，存在误差。（3）模型假设过于理想化，和实际情况存在一定的差异。八、模型推广本文所建的模型可以作为太阳影子定位方面的参考，具有借鉴意义；可以通过该地经纬度、时间、日期、杆高拟合影长随时间的变化曲线，可通过该地的影长顶点坐标与时间求解所处位置等。可以应用在野外求生时地理位置的预测、某地的经纬度预测等九、参考文献1 方荣生，太阳能应用技术.北京：中国农业机械出版社，1985.2 邱国全，太阳时的计算.科普资讯，1998. 3 Wikipedia，https:/en.m.

22、/wiki/Solar_zenith_angle，2015,9,144 Wikipedia，/wiki/Solar_azimuth_angle，2015,9,14.5操晓春，曲彦龄，孙济洲，武琳，郭晓杰，张炜.基于视频中太阳影子轨迹的经纬度估计方法.天津：中华人民共和国国家知识产权局，2009.6地球在线，7韩中庚等.数学建模使用教程.北京：高等教育出版社，2012.十、附录附录1%problem_1程序如下：clcclear allclose allbjtime = 9:1:15;jingdu = 116+23/60+2

23、9/3600n = 295;tao = 360*(n-1)/365;ET = 229.2*(0.000075+0.001868*cosd(tao)-0.032077*sind(tao)-0.014615*cosd(2*tao)-0.04089*sind(2*tao);ts = bjtime - 4*(120-jingdu)+ ET; fai = 39+54/60+26/3600sigma = 23.45 .* sin(2.*pi*(284+n)/365)ts1 = 9:1/6:12;w = 15.*(ts1-12)alpha = asind(sind(fai).*sind(sigma)+cos

24、d(fai).*cosd(sigma).*cosd(w)A1 =acosd(sind(sigma)-sind(alpha).*sind(fai)./(cosd(alpha).*cosd(fai) -180ts2 = 12:1/6:15;A2 =180- acosd(sind(sigma)-sind(alpha).*sind(fai)./(cosd(alpha).*cosd(fai)hold onL1 = 3./tand(alpha)x1 = -sind(A2).*3./tand(alpha);x2 = sind(A2).*3./tand(alpha);y1 = cosd(A1).*3./tan

25、d(alpha);y2 = cosd(A2).*3./tand(alpha);axis(-6,6,-1,6)plot(x1,y1,*)plot(x2,y2,*)plot(0,0,ro)xlabel(x)ylabel(y) 附录2程序如下：clcclear allclose allbjtime = 9:1/6:15;jingdu = 116+23/60+29/3600n = 295;tao = 360*(n-1)/365;ET = 229.2*(0.000075+0.001868*cosd(tao)-0.032077*sind(tao)-0.014615*cosd(2*tao)-0.04089*

26、sind(2*tao);ts = bjtime - 4*(120-jingdu)+ ET; fai = 39+54/60+26/3600sigma = 23.45 .* sin(2.*pi*(284+n)/365)ts1 = 9:1/6:15; w = 15.*(ts1-12)alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w)L=3./tand(alpha)plot(bjtime,L)xlabel(北京时间)ylabel(影子长度)grid on附录3时间、太阳方位角与影子长度的关系图：将三维图旋转一定角度后

27、可得到时间、太阳方位角与影子长度的关系图：太阳方位角与影子长度的关系图：clcclear allclose allbjtime = 9:1:15;jingdu = 116+23/60+29/3600n = 295;tao = 360*(n-1)/365;ET = 229.2*(0.000075+0.001868*cosd(tao)-0.032077*sind(tao)-0.014615*cosd(2*tao)-0.04089*sind(2*tao);ts = bjtime - 4*(120-jingdu)+ ET; fai = 39+54/60+26/3600sigma = 23.45 .*

28、sin(2.*pi*(284+n)/365)ts = 9:1/6:15;w = 15.*(ts-12).*pi./180alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w) A1=acosd(sind(sigma)-sin(alpha).*sin(fai)./(cos(alpha).*cos(fai) hold ongrid onts = 9:1/6:15;w = 15.*(ts-12)alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*co

29、sd(w)L1 = 3./tand(alpha)plot3(ts,A1,L1,r*-)xlabel(时刻)ylabel(方位角/度)zlabel(影子高度/米)附录4时间、太阳高度角、太阳方位角以及相应影子长度的数据如下表：时间9:009:109:209:309:409:5010:0010:1010:2010:30高度角23.319524.739126.10827.422428.678629.872931.00132.05933.04333.948A140.021139.943339.869839.800639.735839.675339.61939.56839.5239.477L16.959

30、396.510756.121615.782065.484465.222884.99264.794.6124.4563时间10:4010:5011:0011:1011:2011:3011:4011:5012:0012:10高度角34.771635.508636.155836.709737.167337.526137.78437.93937.99137.939A139.439239.405339.375939.351139.330739.314939.30439.29739.29439.297L14.3214.204514.105634.023393.957043.905993.86983.848

31、23.84113.8482时间12:2012:3012:4012:5013:0013:1013:2013:3013:4013:50高度角37.783937.526137.167336.709736.155835.508634.77233.94833.04332.059A139.303639.314939.330739.351139.375939.405339.43939.47739.5239.568L13.869823.905993.957044.023394.105634.204514.3214.45634.6124.79时间14:0014:1014:2014:3014:4014:5015:

32、00高度角31.00129.872928.678627.422426.10824.739123.32A139.619239.675339.735839.800639.869839.943340.021L14.992645.222885.484465.782066.121616.510756.9594附录5%Problem_2clcclear allclose allbjtime = 14.70:0.05:15.70; n = 108;sigma = 23.45.*sind(2.*pi.*(284+n)./365);tao = 360*(n-1)/365;a=0;b=0;c=0;d=0;ET =

33、 229.2*(0.000075+0.001868*cosd(tao)-0.032077*sind(tao)-0.014615*cosd(2*tao)-0.04089*sind(2*tao);X = 1.0365 1.0699 1.1038 1.1383 1.1732 1.2087 1.2448 1.2815 1.3189 1.3568 1.3955 1.4349 1.4751 1.516 1.5577 1.6003 1.6438 1.6882 1.7337 1.7801 1.8277;Y = 0.4973 0.5029 0.5085 0.5142 0.5198 0.5255 0.5311 0

34、.5368 0.5426 0.5483 0.5541 0.5598 0.5657 0.5715 0.5774 0.5833 0.5892 0.5952 0.6013 0.6074 0.6135;a=0;min=inf; %step_1for jingdu = 0:30:360 ts = bjtime - 4.*(120-jingdu)+ET; w = 15./(ts-12); for fai= -90:30:90 alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w); A = acosd(sind(sigma)

35、.*cosd(fai)-cosd(w).*cosd(sigma).*sind(fai)./cosd(alpha); for h = 2:0.1:3 x = sind(A).*h./tand(alpha); y = cosd(A).*h./tand(alpha); c=c+1; for i=1:21 D = (x(i)-X(i)2+(y(i)-Y(i)2; end s = sum(sum(D(1:1,:); if smin min = s ; d = jingdu fai h; end end end endmin,d%step_2for jingdu = 90:5:150 ts = bjtim

36、e - 4.*(120-jingdu)+ET; w = 15./(ts-12); for fai= -30:5:30 alpha = asind(sind(fai).*sind(sigma)+cosd(fai).*cosd(sigma).*cosd(w); A = acosd(sind(sigma).*cosd(fai)-cosd(w).*cosd(sigma).*sind(fai)./cosd(alpha); for h = 2:0.1:3 x = sind(A).*h./tand(alpha); y = cosd(A).*h./tand(alpha); c=c+1; for i=1:21 D = (x(i)-X(i)2+(y(i)-Y(i)2; end s = sum(sum(D(1:1,:); if smin min = s ; d =