太阳影子定位问题

1太阳影子定位问题摘要如何确定视频的拍摄地点和拍摄日期是视频数据 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 的重要方面，本文通过研究太阳影子定位技术，建立了基于太阳影子定位的数学模型。第一问，首先，建立太阳高度角和影长的 函 关于工期滞后的函关于工程严重滞后的函关于工程进度滞后的回复函关于征求同志党风廉政意见的函关于征求廉洁自律情况的复函 数关系；接着，利用球面几何关系，推导出太阳高度角与太阳赤纬角、当地纬度和时角的函数关系；然后，推导赤纬角和时角的计算方法，参考有关资料用Bpirges赤纬角算法和Lamm时差算法对赤纬角和时角进行修正，提高模型精度；最后，将题目数据代入影长公式，借助matlab求得太阳影子长度的变化曲线，影长范围为3.67,7.64mm。第二问，首先，根据影子长度和角度建立影子顶点坐标与当地经纬度的函数关系；接着，由于影长及坐标轴方向未知，引入相对杆长和方向角差分别消除杆长和坐标轴方向的影响；然后为了消除误差，使用最小二乘法拟合来进行计算；最后，用matlab编程求解当地经纬度，分别求得利用影长与角度拟合的两组结果，通过相互检验得到角度拟合出的结果更好，其解得直杆所处的地点坐标为19.31,110,18NE。第三问，首先仍然根据影长与角度建立影子顶点坐标与当地经纬度、日期的函数关系，再利用最小二乘法进行曲线拟合，但是所得结果不能满足验证要求，其原因在于仅考虑影长或角度的关系而忽视了日期与经纬度之间的作用关系；因此，选择利用影子的顶点在地面上划过的弧长建立新的函数关系，该函数关系既体现影子长度又体现影子角度；最后，以该函数关系为基础用最小二乘法进行拟合，得到直杆所处的地点和日期，分别用影长和角度的函数关系对其解进行检验，结果可以相互验证，则得到附件二直杆所处的地点坐标为38.99,76.83NE，，日期为2015年7月26日；附件三直杆所处的地点坐标为29.07,112.30NE，，日期为2015年1月21日。第四问，首先进行视频处理，每2分钟取视频一帧获得21张图像；然后，对图像进行裁剪、增加灰度对比度、锐化、边缘检测等处理，得到影子端点的图像坐标；接着，利用成像原理，根据图像坐标求得实际坐标；最后，利用前文方法拟合得到其地点和日期。直杆所处的地点坐标为40.18,122.91NE；当日期未知时，直杆所处的地点坐标为39.86,118.89NE，，日期为2015年7月26日。关键词：太阳高度角太阳方位角最小二乘法相对弧长图像识别2一、问题重述如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。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)所有影长数据均在2015年采集(3)不考虑空气折射对影子长度与方向的影响三、符号说明符号说明L直杆高度l直杆影子长度赤纬角时角F太阳方位角x轴与当地正南方向的夹角3四、问题分析问题一，本问是要求太阳影子长度的变化曲线，首先，考虑到将杆长与影子长度联系起来的是太阳高度角，所以如何求出准确的太阳高度角成为本文的关键；接着，尝试利用几何关系，推导太阳高度角，发现在当地纬度确定的情况下，只要求出太阳赤纬角和时角就可以；然后，我们考虑如何推导出太阳赤纬角与时角的求法，在查阅相关资料后，我们对这两个角度的求法进行了改进；最后，利用改进后的算法，考虑借助matlab求得太阳影子长度的变化曲线。问题二，本问是要在建立以直杆底端为原点的坐标系，水平地面为xy平面的条件下，已知影子端点的坐标来求当地的经纬度。首先考虑根据影子长度建立影子顶点坐标与当地经纬度的函数关系，由于本问中杆长未知，考虑利用比商引入相对杆长消除杆长影响；接着，考虑根据影子方向建立影子顶点坐标与当地经纬度的函数关系，需要利用太阳方位角，由于本问中xy、轴方向与当地地理方向的关系并不已知，考虑引入x轴与当地正南方向的夹角为，列出影子顶点坐标与太阳方位角函数关系，再利用做差引入影角变化消除的影响；然后在求解两个函数关系时，为了更好地利用数据，消除误差，考虑利用最小二乘法拟合来进行计算；最后，考虑用matlab编程求解当地经纬度，并对结果进行分析。问题三，本问与第二问类似，都是已知影子端点的坐标来求当地的经纬度，但区别是本问的拍摄日期也未知。首先，考虑仍然根据影长与影角建立影子顶点坐标与当地经纬度、日期的函数关系，再利用最小二乘法拟合求解；接着，考虑对影长与影角结果进行分析，发现解不能互相验证，考虑其原因在于单纯考虑影长或角度的关系时，日期与经纬度具有耦合作用；然后，考虑找一种既能体现影长又能体现方向的指标，通过这个指标来进行拟合，才能找到正确的位置和日期。考虑到影子的顶点在地面上划过的弧长既包含长度关系又包含角度关系，可以以弧长建立新的函数关系；最后，考虑以该函数关系为基础用最小二乘法进行拟合，来得到直杆所处的地点和日期，并且可以分别用影长和角度对其进行检验。问题四，本问是要通过视频来求当地的经纬度与拍摄日期。首先，考虑进行视频处理，将视频录入matlab中，并将视频每2分钟取一帧，获得21张图像；然后，考虑进行图像处理，对图像进行裁剪、增加灰度对比度、锐化、边缘检测等处理，得到影子端点的图像坐标；接着，考虑利用成像原理，以图像坐标求得实际坐标；最后，考虑利用第三问中的弧长函数关系来得到直杆所处的地点和日期。4五、模型的建立与求解5.1问题一模型的建立与求解5.1.1引入概念1.太阳直射点S太阳直射点是地球和太阳中心连线与地球 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 面的交点。该点的太阳高度角为°90。太阳直射点随地球的自转和公转不断变化，经度坐标以天为单位周期变化，纬度坐标以年为单位周期变化。2.赤纬角赤纬角又称太阳赤纬，它是地球赤道平面与太阳和地球中心的连线之间的夹角，即太阳直射点的纬度。赤纬角是由于地球公转造成的，以年为周期在南回归线与北回归线的范围内移动，即''[2326,2326]，促成了季节的变化。3.时角时角是当地子午圈与过太阳直射点的子午圈之间的夹角，定义当地子午圈在东为正，在西为负，范围°°(180,180]。时角与研究的地点和时间均相关。5.1.2模型的建立hLD图1太阳高度角地面示意图图2太阳高度角天体示意图由图1易知cotlLh，求影长的关键在于求当时当地太阳高度角h，而h可根据图2进行求解。如图2，S为太阳直射点球，N为北极点，D为研究对象，则赤纬角、时角、研究对象纬度分别由、、所示，太阳高度角即为球面弧DS的余角，球面三角形DSN中，边°90SN，°=90DN，=SND，由球面三角形余弦公式可得：°°°°°°cos(90)=cos(90)cos(90)sin(90)sin(90)sin(90)h5化简可得：sin=sinsincoscoscos(1)h从（1）式中可看出，若当地纬度已知，只需求得赤纬角、时角，即可求得当地太阳高度角h，进而求出影子长度l。1.赤纬角的计算赤纬角是由于地球赤道平面与公转平面不重合产生的，它是一个近似以年为周期变化的量，其变化范围为2326。赤纬角全年在南北回归线之间基本呈正弦曲线周期移动。每年夏至6月21日或22日赤纬达到最大值2326，太阳位于地球北回归线正上空，是北半球日照时间最长的一天。随后赤纬角逐渐减少，至秋分9月21日或22日赤纬角等于°0时全球的昼夜时间均相等为。冬至12月21日或22日赤纬角减至最小值2326，此时阳光斜射北半球，北半球昼短夜长。当春分3月21日或22日赤纬角又回到°0，如此周而复始形成四季[1]，如下图：冬至春（秋）分夏至MNPQO(O＇)Smax2r0O’M春分Ot公1图3赤纬度变化图图中箭头1表示上图观察方向，箭头2表示下图观察方向，S为过太阳中心平行于地球赤道平面的参考平面，PQMN、分别是不同位置时地球中心距参考平面的距离，其中PQ表示夏至日时地球到参考平面的距离，其值0maxsinPQr，其中0r为日地距离，6max为回归线纬度。则任意时刻赤纬度0arcsinMNr，max0max'sinsinsinMNOMrt公，其中公为地球的公转角速度，由于地球公转平均用时365.242天，则360/365.2420.9856/dd公，t为地球离开春分点的时间（单位：天）。综上，赤纬度计算公式为：maxarcsin(sinsin)(2)t公2.时角的计算时角是当地子午圈与过太阳直射点的子午圈之间的夹角，设其经度分别为和0，则有0。在一天即24小时内，太阳直射点绕地球转过一周，则其运动角速度等于地球自转角速度，360/2415/hh自，又因为在北京时间12:00时，°120经线过太阳直射点，则在北京时间t时，120-1215t（）经线过太阳直射点，即0=1201512t（）。则时角的计算公式为：120+1512(3)t（）3.对计算公式的修正事实上，由于地球公转轨道不是正圆，而是一个椭圆，日地距离0r及公转角速度公都不是常数，地球日的时间也受公转的影响，赤纬角及时角的精确求解并不如此简单。历史上对赤纬角及时角的求解有多种方法，求解赤纬角的方法包括Cooper法、Spencer法、Stome法、Bpirges法，求解时角的方法包括Lamm法、hillierW法、pencerS法、oolfW法。查阅相关文献，对以上方法进行简单介绍。（1）赤纬角计算方法：①Cooper法284=23.45sin(2)365n式中，n为该日期在一年中的日期序号，例如，当日期为1月1日时，1n，当7日期为3月22日时，81n。②Spencer法0.0069180.399912cos0.070257sin0.006758cos(2)0.000907sin(2)0.002697cos(3)0.00148sin(3)式中，为日角，2(1)365n（n为该日期在一年中的日期序号）。③Stome法12(173)sin0.39795cos365.242n④Bpirges法0.372323.2567sin()0.1149sin(2)0.1712sin(3)0.7580cos()0.3656cos(2)0.0201cos(3)wtwtwtwtwtwt式中，360365.2422w，01tnn，078.8010.2422(1969)0.25(1969)nyearNTyear（2）时角计算方法：当太阳中心连续两次通过地面静止观察者子午圈时称为一个太阳日。为了弥补由于地球绕日公转的速度不均匀以及赤道与黄道不在同一平面，导致一年中太阳日长短不一这一缺点，得到一个相对均匀的时间单位，人们引入一个假想的太阳，称为平太阳[2]。平时生活所用的时间系统是基于平太阳时系统，所以在实时计算太阳时角时，关键是要将平时转换为视时。视时与平时的差称为时差eotequationoftime。所以，求解时角时应用时差修正当地时间在进行计算。与太阳赤纬相同，时差每时每刻都在变化，我们只能用一些不同精度的近似公式进行推算。①Wloof法.258cos7.416sin3.648cos(2)9.228sin(2)eotoxxxx式中，360(1)365.242nx②Spencer法80.0000750.001868cos0.032077sin0.014615cos(2)0.04089(2)eot式中，2(1)365n③Whillier法9.87sin(2)7.53cos1.5sin(10)eotBBB式中，2(81)364nB以上三种时差估算方法均以一年为周期，Lamm法选择时差的周期为四年，更贴近实际。④Lamm法5022()cos()sin()365.25365.25kkkkNkNeotNAB式中，N为从每一个闰年开始为1至四年循环的最后一天，kA和kB的参数值见下表:表1kA和kB参数值k/kAh/kBh02.08704e019.28693e1.22291e25.22582e1.56983e31.30773e5.16023e42.18673e2.98233e51.51004e2.34634e以上是求解赤纬角和时角的典型算法，杜春旭等人在《一种高精度太阳位置算法》一文中对上述算法进行了分类比较，现引述其相关结论[2]。9图4四种赤纬角算法误差比较上图显示近年每天世界时0时四种求解赤纬角算法的计算误差，可以看出，Bpirges法误差最小，最大计算误差为0.025。图5四种时差算法误差比较上图显示近年每天世界时0时四种求解时角算法的计算误差，可以看出，Lamm法误差最小，最大计算误差为0.15。综上所述，Bpirges赤纬角算法和Lamm时差算法精度远优于其他算法，本文将采用上述算法以提高模型的准确性。5.1.3模型的求解编写Bpirges赤纬角算法和Lamm时差算法的程序，通过matlab求解赤纬角和时角，进而求得影长。下图为最终的到的2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线：10图6太阳影子长度的变化曲线从上图中可以看出，在北京时间9:00-15:00之间，位于天安门广场的3米高直杆在9:00时影长有最大值，为7.64m，随着时间推移，影长在9:00-12:18逐渐减小，到12:18影长达到最小值，为3.67m，随后，影长在12:18-15:00逐渐增大。5.2问题二模型的建立与求解5.2.1引入概念1.太阳方位角F太阳直射点和当地的大圆连线与当地南经线所成的球面角叫做太阳方位角，规定其逆时针为正，范围[0360)，，如图中F所示：图7太阳方位角示意图112.方位角的求解同样在球面三角形DSN中运用球面三角形余弦公式可得：°°°°°°cos(90)=cos(90)cos(90)sin(90)sin(90)cos(180)hhF化简可得：sinsinsincos=(4)coscoshFh5.2.2模型的建立1.建立坐标系根据题目要求，建立以直杆底端为原点的坐标系，水平地面为xy平面，但是xy轴的方向与地理方向的关系并不已知，直杆垂直于地面。2.构建影子顶点坐标与当地经纬度的函数关系（1）根据影长建立关系由影子顶点坐标可计算得某时刻直杆影长：22iiilxy根据太阳高度角亦可计算得某时刻直杆影长：cotiilLh则22cot(5)iiiLhxy在(5)式中iixy、已知，ih可根据当地经纬度求得，而杆长L未知。所以引入相对杆长11cot=cotiiilhllh消去杆长影响。最终建立的方程为：2222111cot(6)cotiiixyhhxy（2）根据影子方向建立关系由于坐标系已建立，根据影子顶点坐标，还可以获得影子的方向信息，即影子顶点的坐标可与太阳方向角建立联系。又由于xy、轴方向与当地地理方向的关系并不已知，所以我们假设x轴与当地正南方向的夹角为，逆时针为正，其范围(0,360]。设某时刻影子与x轴夹角为i，则12arctaniiiyx又因为当地影子方向与太阳方向成180，则(arctan180)mod360(7)iiiyFx在（7）式中iixy、已知，iF可根据当地经纬度求得，而未知。所以引入影角变化11iiiFFF消去杆长影响。最终建立的方程为：111arctanarctan(8)iiiyyFFxx可见我们可以根据某一时刻影子顶点的坐标，与当地经纬度建立两个方程，从而解得当地经纬度坐标。但是影子顶点坐标为测量值，各组数据均存在误差，而且只用一组数据求得的结果具有偶然性，没有充分利用所给数据。3.数据拟合为了更高效地利用数据，消除误差影响，提高计算精度，本文采用了最小二乘法对数据进行拟合。下面简要介绍基本最小二乘法的统计学原理:设物理量y与变量12,,,lxxx间有基于1201(,,,;,,,)lnyfxxxaaa的依赖关系。其中01,,,naaa是方程中需要确定的1n个参数。最小二乘法就是通过1mmn个实验点12(,,,,)(1,2,,)iiilixxxyim确定出一组参数值01(,,,)naaa,使由这组参数得出的函数值12101(,,,;,,,)iiinyfxxxaaa与实验值iy间的偏差平方和01(,,,)nsaaa取得极小值。最小二乘法一般采用多点测量来减小随机误差,使待求参数的个数小于方程式个数,即1mn。通过用最小二乘法进行统计处理,使未知数个数和方程个数相等,再进行求解得出01,,,naaa[3]。可见，用最小二乘法拟合需要建立一个关系(,)(,,)xyft，本文之前建立的有关lF、的关系均可以满足其要求。我们选择的做法是先利用其中一个关系进行拟合，并将拟合得到的结果代入另一个关系中检验，通过偏差平方和观察其吻合程度，选择吻合程度较高作为结果。135.2.3模型的求解使用matlab编程，分别用长度和方向关系根据最小二乘法拟合数据，即可求得对当地的经纬度，结果如下：图8拟合结果图可见，用两种关系的拟合结果都不错，具有较高的吻合度。再将拟合结果代入另一个关系中进行检验，结果如下：图9检验结果图可见，虽然利用影长关系所求解的拟合度更高，但其检验结果并不理想，而利用方向关系所求解在拟合和检验中都具有很高的吻合度。所以，该直杆所在地坐标为(19.31,110.18)NE。根据该坐标可以反求得到103.4，即x轴方向为东偏北13.4，14并不是正东方向。下图为直杆所在地坐标在地图中的位置：图10第二问直杆地理位置5.3问题三模型的建立与求解在上一问中，本文先利用根据影长建立的关系进行拟合，并将拟合得到的结果代入根据影子方向建立的关系中检验，通过偏差平方和观察其吻合程度，选择吻合程度较高作为结果。利用影子长度与方向都可以有效地解决影子定位问题，两者相互验证，定位的准确度较高，因而，本问第一步先采用第二问的方法对问题进行求解。由于地球处于公转轨道的同一位置时，不同日期的影长和角度相同，故相同的情况会在各年重复出现，假设题目附件中的数据都是在2015年测的。5.3.1使用第二问模型求解在本问中当地的经纬度及日期未知，所以位置参数为经纬度与日期三个变量根据第二问的方法，依旧使用以下函数关系式：（1）根据影长建立关系2222111cotcotiiixyhxy（2）根据影子方向建立关系15111arctanarctaniiiyyFFxx用matlab编程，以题目附件二中的数据为例，分别用影长和角度关系根据最小二乘法拟合数据，求得当地的经纬度及日期，结果如下：图11利用影长拟合结果图12利用角度拟合结果由上图可以看出，直接使用第二问的方法所得的解不能互相验证，其原因在于单纯考虑影长或角度的关系时，日期与经纬度具有耦合作用，以影长因素为例，影长由太阳高度角决定，太阳高度角与当地纬度和赤纬度的差相关，当日期变化时，日期改变对影长的影响可以通过改变经纬度来抵消，如下图：16图13日期与地点耦合关系图这样就会产生一系列不同日期和地点的组合，并能使影长满足题中所给条件，而在这些日期与纬度的组合中，只有一点可以满足角度条件，但由于在拟合过程中单纯考虑影长条件，故无法找出这一正确组合。相应的单纯考虑角度因素，也会出现相应的问题，所以，解决此题的关键在于找到一种既能体现影长又能体现方向的指标，通过这个指标来进行拟合，才能找到正确的位置和日期。通过以上结果可以看出，在增加一个未知参数后，单独利用影子长度与方向的方法的可靠性均降低了，这样的结果并不能令人信服。5.3.2建立新的函数关系并求解影子顶点与原点的连线反映影子的长度，连线与x轴的夹角反映影子的方向，而影子的顶点在地面上划过的弧长既包含长度关系又包含角度关系，下面根据弧长建立新的函数关系。(,)iixyisil1ili1i11(,)iixyxy图14弧长关系图17如图所示，由影子顶点坐标可以计算得第i段弧长：2211()iiiiisxxyy再由余弦定理可以计算得第i段弧长：221112cos()iiiiiiisllll则222211111()2cos()(9)iiiiiiiiiixxyyllFFll对弧长进行归一化处理得到相对弧长，消去杆长的影响2211121111cotcotcotcot2cos()cotcotcotiiiiiiiishhhhsFFlhhh则，建立的新函数关系式为2222111112111cotcotcotcot()2cos()(10)cotcotcotiiiiiiiiiihhhhxxyyFFhhh上式既能体现影长因素又能体现角度因素，以该式为基础用最小二乘法进行拟合，即可得到直杆所处的地点和日期。结果如下：图15附件二数据拟合结果分别用影长和角度对其进行检验，得到如下结果：18图16附件二数据检验结果从图中可看出，该结果既满足影长条件，又满足角度条件，所以附件2中点的日期和坐标为：201572638.9976.83)NE年月日，(，。同理，使用附件3中的数据可得到其日期和坐标为：2015121年月日，(29.07112.30)NE，，其检验结果见附录。下图为附件二和附件三中直杆在地图中的位置：图17附件二直杆地理位置图18附件三直杆地理位置5.4问题四模型的建立与求解本问不同于前三问，没有直接给出直杆影子顶点的坐标，而是提供了视频，视频中具有大量的信息，包括天气、植被、过往的行人车辆，以及每一帧中都可以提取的时间和影子顶点的信息，结合已做的工作，提取有关影子顶点和时间的信息，并用上文方法求解。195.4.1视频处理这里以matlab为工作平台，将题目所给视频导入到matlab中。下图为导入后得到的结果：图19视频文件信息从图中可知，该视频文件时长40分40秒，帧率25Hz，即每分钟有1500帧图像。由于附件1-3中都给出21组数据，所以以两分钟（3000帧）为间隔对视频进行截取，得到21幅视频截图，以此为基础提取影子顶点和时间信息。图20视频截图5.4.2图像处理首先通过观察的方法在图上点选直杆底部像素，读取其图像坐标为(868,880)，观察法受人为影响较大，具有一定误差，但误差在可接受范围内。然后，提取影子顶点坐标，有以下几个步骤：1.图像裁剪通过对图像的观察，发现影子顶点只在867,1700,780,930xy范围内移动，如图所示20故以此范围对图像进行裁剪以降低计算量和计算复杂度。2.灰度处理用数字图像的直方图代替图像的灰度密度函数。灰度直方图是一个离散函数，表示数字图像每一灰度级与该灰度级出现频率的对应关系。为了使图像清晰，将图像的灰度范围拉开，让灰度频率较小的灰度级变大，让灰度直方图在较大的动态范围内趋于一致。用图像(,)fxy的直方图代替灰度的分布密度函数()fpf，直方图均化后的图像g为[4][5]：0[]()ffgTfpudu直方图均化后的图像(,)gxy在该处的灰度ks为0[]klkklnsTrN上式中，N为像素总数，L为灰度级的个数，kn为频数，kr为(,)fxy在像素(,)xy处的灰度。3.图像锐化处理为凸显影子边缘，将图像转换为黑白图像：4.边缘检测在图像中，边界表明一个区域的结束，另一个区域的开始。边缘检测利用物体差异来实现，包括灰度、颜色或者纹理特征。边缘检测实际上就是检测图像特征发生变化的位置。边缘检测有两个内容：抽出边缘点，剔除间断点，并将边缘连接成完成的线：215.求影子端点图像坐标图像处理至此，影子端点已十分清晰，只需计算白色图像最右端点的坐标即可。5.4.3坐标转化上一步中所得的坐标为图像的像素坐标，想要求解待求量，还需将图像的像素坐标转化成实际坐标，而坐标的转化要涉及到坐标系的转化。这里主要有：1.图像坐标系：在图像上以摄像机光轴与图像平面的交点为原点的坐标系。该坐标系下有两种坐标，一种是以像素为单位的图像坐标，另一种是以毫米为单位的图像坐标，分别以(,)uv和(,)xy来表示，若设每一个像素在x轴与y轴方向上的物理尺寸为,xy，则(,)uv和(,)xy有如下关系：00xuuxyvvy2.摄像机坐标系：以摄像机光心为原点的坐标系，坐标以'''(,,)XYZ来表示，摄像坐标系原点与图像坐标系原点间的距离为焦距f。3.实际坐标系：以实际地物为原点的坐标系，坐标以(,,)XYZ来表示。'''(,,)XYZ和(,,)XYZ的关系是：'''11XXYYZZ其中是44旋转矩阵。这里，本文采用线性成像模型，示意图如下：22Z’P(X’,Y’,Z’)O1pxyX’Y’OZXY图21线性成像示意图从图中可以看出如下的比例关系：''''fXxZfYyZ将三种坐标系下的坐标结合起来，就得到了图像的像素坐标与实际坐标之间的关系：'11XuYZvZ表示一个34的旋转矩阵。在本问中已知杆长，就相当于有了空间上的已知点，利用这些已知点就可以求出旋转矩阵，从而根据影子端点的像素坐标，得到它的实际坐标。用matlab计算，得到影子端点的实际坐标如下表所示：23表2影子端点的实际坐标时间XY8:542.31810.47108:562.29460.41208:582.25920.41209:002.23560.35309:022.23000.35309:042.17800.29509:062.15320.26509:082.11780.23609:102.85400.23609:122.61900.17709:142.26500.17709:162.11800.11809:181.98230.11809:201.95880.59009:221.91750.59009:241.8969-0.29009:261.8740-0.59009:281.8468-0.59009:301.8233-0.59009:321.82700.47109:341.77610.4120根据上表，应用前文的方法解得当地坐标(40.18,122.91)NE，根据该坐标可以反求得到杆长1.72L，与题中数据相近。其拟合结果和地理位置如下：24图22第四问拟合结果图23第四问直杆地理位置六、模型评价6.1模型的优点（1）查阅资料，利用了Bpirges赤纬角算法和Lamm时差算法对求解赤纬角和时角的方法进行了改进，得到了较为精确的结果，提高了模型的准确性。（2）利用影子的顶点在地面上划过的弧长建立新的函数关系，对问题进行求解，能够同时反映长度关系和角度关系，提高了模型的可靠性。（3）在处理图像时，通过图像裁剪、灰度处理等工作使图像分析预处理数据时更加简便、清晰。6.2模型的缺点模型没有考虑大气层对光的折射，导致分析影长与角度时可能会出现误差。七、参考文献[1]互动百科，赤纬角，http://www.baike.com/wiki/%E8%B5%A4%E7%BA%AC%E8%A7%92，2015.09.12。[2]杜春旭，王普，马重芳，吴玉庭，申少青，一种高精度太阳位置算法，新能源工艺，2008（2）：41-48,2008。[3]宋殿瑞，宋文臣，刘鹏振，最小二乘法探讨，青岛化工学院学报，1998（3）：296-301,1998。[4]张毓晋，图像处理和分析，北京：清华大学出版社，1999.3:43-49。[5]阮秋琦，数字图像处理学，北京：电子工业出版社，2001.1:180-208。25八、附录1.附件三拟合及检验结果：图24附件三数据拟合结果图25附件三数据检验结果2.相关程序利用影长拟合：l1=sqrt(X1.^2+Y1.^2);ll1=l1/l1(1);pl0=[philamda];[PL1,var11]=lsqcurvefit(@ratio_l,pl0,th1,ll1);theta1=atan(Y1./X1);dtheta1=(theta1-theta1(1))*180/pi;[PL12,var12]=lsqcurvefit(@d_F,pl0,th1,dtheta1);ll1_check=ratio_l(PL1,th1);plot(th1,ll1);26holdonplot(th1,ll1_check,'*r');title('利用影长拟合结果','fontsize',20);xlabel('时间','fontsize',20);ylabel('相对影长','fontsize',20);lege=legend('原始数据','拟合数据');set(lege,'fontsize',20);text(14.8,1.5,['偏差平方和',10,num2str(var11)],'fontsize',20);text(15.1,1.1,['(',num2str(PL12(1),4),',',num2str(PL12(2),5),')'],...'fontsize',20);利用角度拟合：theta1=atan(Y1./X1);dtheta1=(theta1-theta1(1))*180/pi;pl0=[philamda];[PL12,var12]=lsqcurvefit(@d_F,pl0,th1,dtheta1);dth1_check=d_F(PL12,th1);F1_check=calc_F(T,PL12(1),PL12(2),th1);Delta=F1_check-180-theta1*180/pi;Delta_bar=mean(Delta);plot(th1,dtheta1);holdonplot(th1,dth1_check,'*r');title('利用角度拟合结果','fontsize',20);xlabel('时间','fontsize',20);ylabel('影角变化','fontsize',20);lege=legend('原始数据','拟合数据');set(lege,'fontsize',20);text(15.3,-3,['偏差平方和',10,num2str(var12)],'fontsize',20);text(14.75,-7,['(',num2str(PL1(1),4),',',num2str(PL1(2),5),')'],...'fontsize',20);