太阳影子定位

太阳影子定位 ------------------------------------------------------------------------------------------------ 太阳影子定位 摘要 针对问题一，建立直杆影子长度变化的数学模型，通过分析影子长度关于经度，纬度等参数的变化规律以及与太阳高度角、太阳赤纬和太阳时角等中间因素存在的联系，建立起数学模型。根据所给的日期、时间、地点、杆高等数据，使用模型绘制位于天安门，2015年10月22日的直杆影子长度的变化曲线，得出最小长度出现在正午12时，最短长度约为3.68米。 针对问题二，建立平面坐标，根据给定的平面影子顶点坐标，用最小二乘法思想进行曲线拟合，得出影长随时间的变化规律。通过几何原理计算太阳高度角、赤纬角、地理纬度、时角和太阳方位角之间的关系，计算出目标地的地理纬度。同时选取北京正午十二时作为一个参考值，联系拟合的曲线得出目标地太阳影子最短的时刻，与参考值进行对比得出经度。经纬度的双锁定确定直杆所在地的具体位置大致是广东省茂州市。 针对问题三，建立与第二问类似的目标规划模型，但是加上日期因素，为了提高遍历速度，使用粒子群算法减小时间复杂度。计算得出的结果:附件二地点坐标(80.34?E,33.32?N)，8月13日或4月27日，约在西藏山南地区;附件三地点坐标(113.21?E，23.54?N)，12月26日或12月12日，约在广西河池市。 针对问题四，对视频进行处理以得到必须的数据，由于拍摄时角—————————————————————————————————————— ------------------------------------------------------------------------------------------------ 度的存在，得到的影长并非实际长度，而是投影长度，所以采用基于Hough变换和透视变换的图像矫正法，对斜视图像矫正，得出实际长度，然后将得到的数据代入第二问中的模型，得出的结果是视频拍摄地点在内蒙古包头市境内;在拍摄日期未知的情况下，将图像变换得出的实际影长代入第三问中的基于粒子群算法的规划模型里，可得出日期为6月11日和7月13日. 对于模型的推广，根据某地的阳光信息与地理信息，可以应用到建筑行业日照 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 ，设计住宅布置、建筑间距、房屋朝向等方面，甚至还可以应用到农业大棚养殖安排蔬菜的种植位置等。 关键词:二次函数拟合 最小二乘法 粒子群算法 图形矫正法 1 目录 工贸企业有限空间作业目录特种设备作业人员作业种类与目录特种设备作业人员目录1类医疗器械目录高值医用耗材参考目录 一、问题重述 1.1问题的背景 1.2问题的提出 4.1 问题一分析 4.2问题二分析 4.3问题三分析 4.4 问题四分析 5.1 问题一模型建立与求解 5.1.1 概念引入 5.1.1.1 太阳高度角 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 5.1.1.2 太阳赤纬 5.1.1.3 太阳时角 5.1.2模型建立 5.1.3模型求解 5.2 问题二模型建立与求解 5.2.1 概念引入 5.2.2 模型建立 二、模型假设 三、符号说明 四、问题分析 五、模型建立与求 解 2 5.2.2.1纬度求解模型 5.2.2.2 经度求解模型 5.2.3 模型求解 5.3问题三模型建立与求解 5.3.1 模型的建立与求解 5.3.2粒子群算法目标规划模型 5.3.3粒子群算法的简介 5.3.4粒子群算法模型的求解过程 5.3.5 模型的分析比较 5.4问题四模型建立与求解 5.4.1 概念引入 5.4.2 模型建立 5.4.3 模型求解 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 5.4.4 结果分析 六、误差分析与模型检验 6.1 误差分析 6.2模型的检验 七、模型评价 八、参考文献 附录 3 一、问题重述 1.1问题的背景 对于直立在地面的直杆，其影长与方向会随着日光方向的变化而改变，而且与一天中的时刻有关，因此对太阳光影的变化进行研究以获得时间和地理信息是可行的，特别是使用对图像的处理方法，获得拍摄地点的各种信息，对采光，建筑等 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 提出改进，具有现实应用的意义。 1.2问题的提出 1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场(北纬39度54分26秒,东经116度23分29秒)3米高的直杆的太阳影子长度的变化曲线。 2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶—————————————————————————————————————— ------------------------------------------------------------------------------------------------ 点坐标数据，给出若干个可能的地点。 3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。 4(附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。 如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期, 二、模型假设 假设一:假设地球为一个球体，忽略实际椭圆体对解题的影响; 假设二:忽略海拔、天气等因素对测量和计算造成的影响; 假设三:假设阳光照射到地球上时为平行光; 假设四:假设直杆本身理想笔直，无弯曲情况; 假设五:忽略行星轨道差异; 假设六:太阳影子的长度变化只与太阳照射有关，不考虑大气对太阳的折射; 4 三、符号说明 四、问题分析 问题一要求建立反映影子长度变化的数学模型，引入太阳高度角、太阳赤纬和太阳时角概念，根据三者之间的关系，利用三角函数关系和几何关系，建立不同时刻太阳下影子长度的变化的数学模型，并根—————————————————————————————————————— ------------------------------------------------------------------------------------------------ 据已给的日期、时间、地点、杆高数据，应用所建立的模型画出影子长度变化曲线。 4.2问题二分析 问题二给出了一组数据，要求建立关于直杆所处的地点的数学模型，并应用于附件1中影子顶点坐标数据，得出直杆所在地点。基于直杆影子与时间之间的关系，建立时间与直杆影子顶点坐标的数学模型。同时利用地理上的时区计算方法，来进一步求出该地的经纬度，进而通过经纬度的双重锁定确定所求的地点。 对于纬度求解，我们在问题一中已经有了相应的说明，找出太阳方位角、太阳高度角、时角和太阳赤纬角之间的关系即可求出。根据经线的划分原则，则可知地球自转每分钟为0.25 5 度，然后结合问题一中的方法计算出经度。根据所建立的模型，利用最小二乘法的思维，利用MATLAB软件拟合出直杆影长与时间的关系图像。 4.3问题三分析 对于问题三，可视为已知投影顶点的坐标和投影时刻，要求得出直杆所在位置的经纬度与拍摄的日期。首先，将附件2与3中的影长和时刻数据带入问题一求得的影长变化函数中，得到方程组。其次该问题是一个双目标最优化问题，可以建立以影长差值最小，太阳方位角差值最小为目标的双目标最优化模型对附件2,3 的数据进行求解，得出位置的经纬度与拍摄的日期，其中约束条件考虑经纬度范围。由此，可建立双目标最优化模型，利用粒子群算法优化遍历过程，并对—————————————————————————————————————— ------------------------------------------------------------------------------------------------ 模型的正确性进行检验。 4.4 问题四分析 对问题四所给的视频附件进行一定的预处理，由于给出视频的清晰度不能满足要求，还需要对图片进行灰度处理，以达到减少计算量和提高计算精确度的目的。 利用MATLAB对处理后的数据进行分析，与问题二的模型相互联系，得出直杆影子长度变换规律，继而得出可能的拍摄地点。 五、模型建立与求解 5.1 问题一模型建立与求解 5.1.1 概念引入 5.1.1.1 太阳高度角 对地球上所选取的某个地点来说，太阳高度角是指太阳光照射的方向和地平面之间的夹角，从专业术语上讲太阳高度角是指某一个地方太阳光线与通过该地与地心相连的地表切线的夹角。 太阳高度角随着地方时间和太阳的赤纬的变化而变化。而且当太阳高度角为90度时，此时太阳辐射强度是最大的;当太阳斜射在地面时，太阳辐射强度就变小。 5.1.1.2 太阳赤纬 太阳赤纬又称赤纬角，是地球赤道平面与太阳和地球中心的连线之间的夹角。赤纬角以年为周期，在+23?26′与-23?26′的范围内移动，成为季节的标志。(图1) 6 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 图1 赤纬角的定义 5.1.1.3 太阳时角 日面中心的时角，即从观测点天球子午圈沿天赤道量至太阳所在时圈的角距离。 以地球为例，在地球上，同一时刻，对同一经度，不同纬度的人来说，太阳对应的时角是相同的。单位时间地球自转的角度定义为时角w，规定正午时角为0，上午时角为负值，下午时角为正值。地球自转一周360度，对应的时间为24小时，即每小时相应的时角为15度。 5.1.2模型建立 引入了太阳高度角、太阳赤纬、太阳时角的相应概念，根据三者之间的定义，我们可以得出太阳高度角随着太阳赤纬和太阳时角的变化而变化的结论。根据三者的空间几何联系，建立合适的模型来反映影子长度变化。(下图2) 图2:高度角、赤纬角、时角关系图 已知太阳高度角随太阳赤纬和太阳时角的变化而变化，三者之间存在以下的运算关系: sin??sin?sin??cos?cos?cos?(1) 其中α为太阳高度角，Φ为直射点地理纬度，θ为太阳赤纬角，ω为太阳时角。 7 通过查阅 资料 新概念英语资料下载李居明饿命改运学pdf成本会计期末资料社会工作导论资料工程结算所需资料清单 ，我们发现，因太阳赤纬值日变化很小，所以一年内任意一天的赤纬角可以通过下式进行计算: —————————————————————————————————————— ------------------------------------------------------------------------------------------------ sin??0.39795cos(0.98563(N?173))(2) 其中，N为日数，自每年第一天计算。 ??15(T?12)??/180(3) 其中，T为时刻，自每日凌晨计算。 将查阅到的太阳赤纬角以及太阳时角的的计算方法代入太阳高度角的计算公式，可以得到以下的方程: sin??0.39795sin?cos(0.98563(N?173)) ?cos?cos?cos((T?12)?/12)(4) 根据题目要求画出北京天安门广场的3米高的直杆的太阳影子直射简易图: 图3:太阳影子直射图 阳光与地面的夹角α即为此地的太阳高度角，AC为目标杆，杆长为3m，在三角形ABC中，由正切定理可知: tan??AC/BC (5) 要想求出直杆AC的影长BC随时间的变化规律，既要求出该地太阳直射角α在指定时间内的变化规律，即可得出影长: BC?AC/tan?(6) 由(4)式可知太阳直射角α: ?=arcsin(sin? ) (7) 8 BC?AC/tan? (8) —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 5.1.3模型求解 已知目标对象为2015年10月22日的北京天安门广场(北纬39度54分26秒，东经116度23分29秒)。根据以上信息可得地理纬度Φ为39度54分26秒，日数N为295天，并且根据目前已有信息可求得该地点赤纬度数。 以时刻为单位，并且联立(1)-(8)式，求得各个时刻的影长数据，通过编写程序以及MATLAB软件的应用，画出对应的影长图像，从而得出太阳直射的直杆影子长度变化曲线。(图4) 图4:直杆影长变化图 根据MATLAB所得图像，可知2015年10月22日北京时间9:00-15:00之间天安门广场(北纬39度54分26秒,东经116度23分29秒)3米高的直杆的太阳影子长度的变化曲线为一开口向上的抛物线，直杆影子长度随时间呈现出先减小后增大的规律，且最小长度出现在正午12时，最短长度约为3.68米。 5.2 问题二模型建立与求解 5.2.1 概念引入 太阳方位角:即太阳所在的方位，指太阳光线在地平面上的投影与当地子午线的夹角，可近似地看作是竖立在地面上的直线在阳光下的阴影与正南方的夹角。方位角以目标物正北方向为零，顺时针方向逐渐变大，其取值范围是0—360?。因此太阳方位角一般是以目标物的 北方向为起始方向，以太阳光的入射方向 为终止方向，按顺时针方向所测量的角度。 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 5.2.2 模型建立 已知目标杆各个时刻于水平地面上在太阳下的太阳影子顶点坐标， 可以得出相对应时刻的 9 太阳影长，进而得到在指定时间内的目标杆太阳影子长度变化规律，联系几何原理以及太阳光射角度等变量，得出直杆在该地太阳影子长度的最小值，联系一定的天文知识，计算出该地的经度和纬度，继而查找出该地所在。(图5) 图5;太阳直射下直杆影子投射示意图 以直杆所在点为目标原点，以东西方向为X轴，南北方向为Y轴，建立XY平面坐标系。模拟出太阳直射下直杆影子在平面上的投射，此时影子顶点为S点，其在XY坐标轴上的投影分别为Xs，Ys。(图中α角为太阳高度角，γ角为太阳方位角) 5.2.2.1纬度求解模型 日升日落，同一地点一天内太阳高度角是不断变化的。日出日落时角度都为0，正午时太阳高度角最大，时角为0，以上的公式可以简化为: sin??sin?sin??cos?cos? 由两角和与差的三角函数公式，可得 (9) sin??cos(???)(10) 此时，发现存在两种情况，因此对太阳位于天顶以北的地区和太阳位于天顶以南的地区分别进行讨论: —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 对于太阳位于天顶以北的地区而言，α=90o,(φ,θ); 对于太阳位于天顶以南的地区而言，α=90o,(θ,Φ); 二者合并，因为无论是(Φ,θ)还是(θ,Φ)，都是为了求当地纬度与太阳直射纬度之差，不会是负的，因此都等于它的绝对值，所以得到正午太阳高度角计算公式: 10 ??90????(11) 因此，只需要知道该地的太阳高度角和太阳赤纬即可得出该地的地理纬度，太阳赤纬角可以由式(2)计算得出。引入太阳方位角γ，要得出太阳高度角α的值，经查阅查阅推导得出太阳方位角、太阳高度角、太阳赤纬角以及太阳时角之间满足下列关系: cos??cos?sin?/sin? (12) 经过(11)(12)式的联立，可以得到以下的等式: ??arccos(cos?sin?/sin?)(13) ??90?????arccos(cos?sin?/sin?) 角，ω为太阳时角。 联立(2)(3)(14)即可得出该地的地理纬度。 5.2.2.2 经度求解模型 地球自传一周需要24小时，即将全球被划分成24个时区，可知地球每小时自传15度，即可推出每分钟地球自转了0.25度，根据这个数据既然可以计算出不同经度地方的时间差，反之亦可从时间差推导出经度差。 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 通过问题一的求解，我们已经知道，直杆的影子最短发生在一天中正午时，通过选取已知纬度和正午时刻的地方作为参考地，利用正午时刻时间差即可求出所求地的经度值。 建立目标地点对比点，选取北京正午时间t。和经度ξ。作为参考点(其中t。为正午12:00, ξ。为东经120o)，所求地正午时刻t1与t。对比计算出两地的时间差，再通过数学计算把差值换算成相应的度数，即可以得出所求地的经度值ξ1。两者有如下关系: (14) 在上列处理中，α为太阳高度角，Φ为直射点地理纬度，θ为太阳赤纬度, γ为太阳方位?1??0?t1?t0 (t1<t。) (15) ?1??0?t1?t0(t1>t。) (16) 为了得出附件1所提供的数据反映的可能地点，需要求出该地的经纬度。因此还要对附件1的数据进行模拟得出目标直杆达到影长最短的时刻，即可求出该地的经度数。 5.2.3 模型求解 为了得出问题中直杆所处的地点，我们需要先对附件1直杆的影子顶点坐标数据进行预处理，对给出的数据在应用到XY坐标系的基础上，通过对两坐标轴投影点的数据处理，得出各 11 个时刻的影长，由此得出影长随时间的变化规律，为了较为直观的体现出来，我们采用图标形式进行展现:(图6) 图6:目标地直杆影长随时间变化图标示意图 利用最小二乘法的思想，通过曲线拟合，得出目标直杆影子长度随时间的变化关系满足下式:(图7) —————————————————————————————————————— ------------------------------------------------------------------------------------------------ Y?0.1481t2?3.7438t?24.1371 (17) 图7:直杆影长拟合曲线 并因此式得出拟合曲线图像，通过拟合函数关系找到使得目标直杆影长达到最短的时间，即要求找到使得Y值达到最小的t值，将得到的t值代入式(3)，从而得到此地的经度为东经 12 111度9秒。 在求解目标地纬度的时候，同样应用最小二乘法思维，分别拟合得出直杆影子顶点的X值随时间的函数关系以及直杆影长顶点的Y值随时间的函数关系 :(图8，9) Yx?0.1480t2?3.6541t?24.1371(18) Yy?0.0055t2?0.0466t?0.0280(19) 图8:X坐标拟合曲线 图9:Y坐标拟合曲线 13 使用同样的方法，在MATLAB软件的运行下得出两函数的拟合曲线图，分别找到使函数值达到最小的t值，将其代入(18)(19)两式，得出此时的影子顶点坐标，根据图5算出此时的太阳方位角γ，并因此得出sinγ，联立(14)式求出此地的地理纬度Φ。Φ为北纬22度11分24秒。 综上分析，通过该模型所最终确定的目标地的经纬度为东经111度9秒，北纬22度11分24秒。经查找，该地大概位于广东省茂名市高州市G65包茂高速旁。 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 图10 分析结果所示地点 5.3 问题三模型建立与求解 5.3.1 模型建立 第三问将日期设为未知量，在前两问模型的基础上，要求求解杆子所在位置的经纬度与日期。由地理知识可以得出，日期影响的唯一因素为赤纬角，即太阳直射点的纬度，有以下两个计算某日赤纬角的公式可供选择: (1)θ=23.5?sin(2?(284?n)) 365 式中，n为当日日期序号，1月1日时n=1,以此类推，日期每增加一天，即n加1。 (2) θ=23.5?-(日期-6月22日)×(23.5?×4/365) 此公式可以大致计算一年中某年太阳直射点所在的纬度值。 计算结果若是正值，则为北纬;若为负值，则为南纬;R为某日日期，R-6月22日为该日与6月22日相差的天数，(23?26′*4/365)为太阳直点一日内移动的纬度距离(假设其移动是匀速的)。直射点经度即太阳高度最大(太阳上中天)的经线，地方时12:00的经线; 14 建立目标函数: (1)太阳方位角的差值平方和最小，要求求出的每个时刻方位角与附件2，3中对应时刻方位角差值平方和最小。 (2)影长差值平方和最小,要求求出的每个时刻影长与附2，3中对应时刻影长差值平方和最小。 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 212??1)?min????(?i??b)????(?i?1?21 (20) ?min?l?(l?l)2??????(?2)?ibi?i?1? 该目标函数的约束条件: (1)影长li要求符合问题一所提出的五个基本因素。 (2)影长li非负。 (3)太阳直射点(赤纬角)在南北纬23.5?之间移动。 (4)经度的范围在东西方向各180?范围之内。 ??1)?li?f(?,?,dt,t,h)??( ? l?0??????????????????(?2) ?is.t??3)? -23.5??θ?23.5???( (21)对于经纬度进行取点处理，进行??4)? -180??ω?180?????( 编程求解: 表1:对于杆子所处地点的日期与误差分析数据表 15 5.3.2粒子群算法目标规划模型 在第三问中，加入的新变量日期，增加了模型的可行域，对于遍历算法的实现增加了时间复杂度，因此，我们引入粒子群算法对求解过程进行优化，降低模型求解的时间复杂度。 5.3.3 粒子群算法的简介 假设有一个D维的目标搜索空间，有n个微粒组成了一个粒子群，其中每个微粒都用一个D维的向量描述，将它的空间位置表示为: mi?(mi1,mi2,???,miD), i=1,2,…,n; —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 这可看做目标优化问题中的一个解，带入适应度函数计算出适应度值可以衡量微粒的优劣;第i个微粒的飞行速度也是一个D维度的向量，记为 vi?(vi1,vi2,???,viD);i=1,2,„,n; 第i个微粒所经历过的具有最好适应值的位置称为个体的历史最好位置，记为 pi?(pi1,pi2,???,piD);i=1,2,…,n; 整个微粒群所经历过的最好位置称为全局历史最好位置，记为 pg?(pg1,pg2,???,pgD);i=1,2,„,n; vij?t?1??vij(t)?c1r1(t)(pij(t)?xij(t))?c2r2(t)(pgj(t)?xij(t)) 粒子群的进化方程可描述为: mij(t?1)?mij(t)?vij(t?1)(22) 其中，i,j 分别表示第i个微粒的第j个维度，t代表第t代，c1,c2是两个加速常量，通常的取值范围是(0,2)，r1~U(0,1),r2~U(0,1)是两个相互独立的随机函数。从上式可以看出，c1可以调节微粒去往自身周边的最好位置，c2可以调节微粒去往整个粒子群所能找到的最好位置。 将这种算法引用到本问中。算法中的“位置”就是本文中不同日期的所有离散的坐标值点:mi?(?,?,d)，即“粒子Zj”就是存放并寻找解的动态变量。为了能够寻找较为合适的优化解，并且避免计算量的增大，本模型中设置了20个“粒子”来进行优化解搜索，即令i=20。由于本问中直杆坐标和日期都是未知的，因此将粒子设为储存三个参数: —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 16 ?Z1?经度???Z2?纬度? ?Z?日期d (23) ?3 5.3.4粒子群算法模型的求解过程 (1) 首先要对所有的粒子进行初始化，即让粒子中的不同参数可行范围内的随机值，如 令??rand(?180,?180),??rand(?90,?90),d?rand(1,365)。这个过程相当于在整个“位置”区域内的某些点散布了寻找最优解的“种子”。在何处散布了”种子“与其后来能够寻找到的“最优解”密切相关。当散布位置不理想时，“最优解”很可能会是局部最优解，而非全局最优解。 (2) 其次计算每个粒子的适应值。本问的适应度值可以使用目标函数来进行计算，但由于最优解是令目标函数最小，为最优解的适应度值应该尽量大。所以对适应度值Q做定义为: 2021?V???i??b??li?lbi??i?1i?1Q???R?(1)2?V?(24) (3) 然后记录这些值，对每个粒子，将其适应度与个体最好历史记录和全局最好记录作比较。如果能够比最好记录更好就将此值设置为最好记录。 (4) 接下来，参照式(22)，对每个粒子的去向进行调整。每次调整时，到要对c1?(0,2)，c2?(0,2)，r1?(0,1)，r2?(0,1)进行重新随机取值，保证粒子行为的随机性:既能有向最优解靠近的趋势，也能存在寻找—————————————————————————————————————— ------------------------------------------------------------------------------------------------ 其他更优解的可能。 (5) 最后，查看迭代次数。为了保证找到的最优解，需要机遇粒子足够多的行动时间寻找解。这里令迭代次数为1000。 对经纬度进行离散化处理后，选取步长为0.01，对其进行MATLAB编程求解: 17 5.3.5 模型的分析比较 将两个模型的结果对比，发现结果较为接近，附件1所得出的结果显示杆子处在西藏山南地区，附件2所得出的结果显示杆子所在地为广西河池市内。且得到的日期关于夏至日6/22与冬至日12/22对称，这是由于太阳的运动是关于夏至日以及冬至日对称的。同时在太阳方向角残差比以及影长残差比这两个参量上，粒子群算法的误差并没有显著提升，最大值为3.1%，在合理范围之内;在时间复杂度上，粒子群算法的求解时间大大缩短，所以当未知参量增加时，使用粒子群算法是较合适的算法。 5.4问题四模型建立与求解 5.4.1 概念引入 透视投影:透视投影即使用中心投影法，可以使得延伸到无穷远的物体在照片平面中获得其有效的延伸范围。例如，在现实世界坐标中存在两平行铁轨，则可以投影为照片平面中的相交直线，且其交点存在于两平行直线的消隐点方向。 在水平地面上，消隐点无限延伸后，即为地平线。并据此引出地平线性质:在平面透视图中，从地平线上引出的任意两条直线都是平—————————————————————————————————————— ------------------------------------------------------------------------------------------------ 行的。 5.4.2 模型建立 在题目附件的视频中，我们发现，在假设摄像机自身不随周围气流而发生抖动的基础上，由于摄像机在拍摄时的角度问题，使得从视频中所获取的影长是真实影长向摄像机方向上的投影值，其值必定小于等于真实影长。由于从视频中无法判断影子随时间的真实转动角度，我们需要通过进行一定的透视变换来还原出影子的真实长度及其变换角度。 对视频进行分析，建立相应的示意图。根据前面的相应分析，认定摄像机水平拍摄，且图中左下角的两条绿色直线在现实中处于水平面中且相互平行，将其延伸在远处找到消隐点，且消隐点在水平方线无限延伸后就是地平线。(图11) 18 图11 视频分析模型 确立直杆立点为O，在视角方向，令地平线与直杆的交点为P，P点即为直杆向正前方的消隐点。同时从底部引出一水平辅助线B，则OP连线L与水平线相互垂直。设此时的直杆影子顶点为R，连接PR得辅助线A并延伸。由地平线的性质可知，此时从P点引出的任意两条直线都是平行的，即直线A与直线L相互平行，则A与B必有一交点，记为Q。则?RQO为直角三角形。 利用霍夫变换提取平行直线，当图像中存在很多条直线的时候，霍夫变换直接使用效果并不好，对算法进行改进，舍弃直线宽度在1—————————————————————————————————————— ------------------------------------------------------------------------------------------------ 以上的直线。 同时，由于照片平面中物体因为成像透视变换而产生形变，需要对其进行透视逆变换使其恢复成原来样子。已知图片空间中二维图像经过透视变化变成另一个图像，过程符合下列公式: m=(ax+by+c)/(ox+py+1) n=(dx+ey+f)/(qx+py+1)(25) 式中，斜视图像点的像素坐标为(x,y),正视图像的像素坐标为(m,n)，a,b,c,d,e,f,o,p,q均为透视变换参数。 为了提高算法的精确度，对图像进行预处理，以达到减少计算难度的目的。利用MATLAB对其进行图像提取，并且对视频中提取产生图片进行灰度处理和高斯滤波处理。可以发现，进行高斯滤波处理后的图像噪点明显减少，减少了计算难度。继续通过MATLAB进行矫正处理，提取出实际情况下的影子长度变化。 经过一系列数据处理，我们得到了不同时刻下的直杆影子长度变换情况，代入到问题二的经纬网模型中去，即可得到视频拍摄地点。 5.4.3 模型求解 对经纬度进行离散化处理后，对其进行编程求解: 19 5.4.4 结果分析 通过表格对比发现当把日期作为未知量对模型进行求解时，可得地点坐标与已知日期所求地点坐标经度上基本无变化，在内蒙古包头市内，纬度相差3.4度，变化幅度为1.8%，由于数据的采集和预处理不可避免的造成了误差，可以接受此误差。所求日—————————————————————————————————————— ------------------------------------------------------------------------------------------------ 期为6/11，视频显示日期为7/13，刚好关于夏至日6/22对称，说明模型可靠。 六、 误差分析与模型检验 6.1 误差分析 本文中涉及了大量的角度，由于对角度进行计算时得到结果的小数点位数是较多的，因此还需要对数据进行预先处理，舍弃结果的小数点后四位，这样一来，对于数据的准确性还是有一点影响。不失一般性而且由于同一天中太阳赤纬角的变化很小，因此本文解答过程中进行处理时，将太阳赤纬角当作一天内并没有变化，这样也会对所要求的结果产生一些影响。 6.2模型的检验 本题中所提供的附件本身就是对所建立的模型的一个检验，给模型提供一些检验数据，从理论上来讲，如果所建立的数学模型合理且正确，那么最终求出的地点和时间也必定符合题目的要求。 七、模型评价 本文所建立的数学模型和实际情况相符合，对于问题有合理的猜想、假设、计算以及检验，具有一定的指导作用，其中在计算直杆太阳影子长度的变化曲线时，结合较为大众所接受的几 20 何原理以及现实中“立竿见影”的原理，建立了太阳位置和影子长度的变化的数学模型。并且在此基础上，可以进行进一步研究，例如可以推广到求解在不同的时间条件下或者不同位置时，太阳对直杆影子的长度变化的影响。 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 整体来说，本 论文 政研论文下载论文大学下载论文大学下载关于长拳的论文浙大论文封面下载 中所建立的几个数学模型具有一般推广性，不仅仅可以求解更为复杂的太阳影子为题，同时，利用太阳影子长度可以推算出当地的经纬度、时间、日期和相对位置的思想，可以推广到建筑行业日照设计，设计住宅布置、建筑间距、房屋朝向等现实应用方面，甚至还可以应用到农业大棚养殖方面，利用不同高度植物采光度不同，合理安排蔬菜的种植位置。 在实际纬度求解过程中，对于数据的处理，我们采用了最小二乘法的思想来进行曲线拟合，得出直杆影长的变化与时间的图像，从而找到目标地点正午的时间，后来经过简洁分析，发现用多项式拟合的方法处理曲线会更加精确，但由于时间问题并未实施。 同时，本文在求解过程中，设置了很多假设，例如不考虑大气对阳光的折射，实际上应该重点考虑当地的大气的折射对当地太阳的折射率，因为不同的地方的太阳光线折射率是不同的，根据所考虑的问题，还需要建立气压与太阳折射率的数学模型，但是因为没有时间，没有完成模型的修正。 总之，对于不同的情况，需要要对模型做适当的修正和补充，针对实际情况选择合适的参数，然后利用这些模型中的核心思维，可以比较灵活的去解决所要求解的问题。 八、参考文献 [1]于贺军，气象用太阳赤纬和时差计算方法研究，气象水文海洋仪器，2006年3期:P23-45 [5]霍俊爽 杨柳 吴兆薇 张春转 王柳行，由物体的影子变化定位—————————————————————————————————————— ------------------------------------------------------------------------------------------------ 其位置的数学模型研究，决策与信息，2015年12期:P25-26 21 附录 影长数据处理 附件一: 北京时间 14:42 14:45 14:48 14:51 14:54 14:57 15:00 15:03 15:06 15:09 15:12 15:15 15:18 15:21 15:24 15:27 15:30 15:33 15:36 15:39 15:42 附件二: 北京时间 12:41 12:44 12:47 12:50 12:53 12:56 12:59 13:02 13:05 13:08 13:11 13:14 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 x坐标(米) -1.2352 -1.2081 -1.1813 -1.1546 -1.1281 -1.1018 -1.0756 -1.0496 -1.0237 -0.998 -0.9724 -0.947 y坐标(米) 0.4973 0.5029 0.5085 0.5142 0.5198 0.5255 0.5311 0.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 y坐标(米) 0.173 0.189 0.2048 0.2203 0.2356 0.2505 0.2653 0.2798 0.294 0.308 0.3218 0.3354 22 影长(米) 1.1496 1.1822 1.2153 1.2491 1.2832 1.3180 1.3534 1.3894 1.4262 1.4634 1.5015 1.5402 1.5799 1.6201 1.6613 1.7033 1.7462 1.7901 1.8350 1.8809 1.9279 影长(米) 1.2473 1.2228 1.1989 1.1754 1.1524 1.1299 1.1078 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 1.0863 1.0651 1.0444 1.0243 1.0046 13:17 -0.9217 13:20 -0.8965 13:23 -0.8714 13:26 -0.8464 13:29 -0.8215 13:32 -0.7967 13:35 -0.7719 13:38 -0.7473 13:41 -0.7227 附件三: 北京时间 x坐标(米) 13:09 1.1637 13:12 1.2212 13:15 1.2791 13:18 1.3373 13:21 1.396 13:24 1.4552 13:27 1.5148 13:30 1.575 13:33 1.6357 13:36 1.697 13:39 1.7589 13:42 1.8215 13:45 1.8848 13:48 1.9488 13:51 2.0136 13:54 2.0792 13:57 2.1457 14:00 2.2131 14:03 2.2815 14:06 2.3508 14:09 2.4213 第一问影子长度变化的程序: clc clear close all syms xy_long 0.3488 0.3619 0.3748 0.3876 0.4001 0.4124 0.4246 0.4366 0.4484 y坐标(米) 3.336 3.3299 3.3242 3.3188 3.3137 3.3091 3.3048 3.3007 3.2971 3.2937 3.2907 3.2881 3.2859 3.284 3.2824 3.2813 3.2805 3.2801 3.2801 3.2804 3.2812 23 0.9855 0.9668 0.9486 0.9309 0.9138 0.8971 0.8810 0.8655 0.8505 影长(米) 3.5331 3.5468 3.5618 3.5781 3.5958 3.6149 3.6354 3.6572 3.6805 3.7052 3.7313 3.7589 3.7881 3.8187 3.8508 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 3.8846 3.9199 3.9569 3.9955 4.0358 4.0779 I=imread('sunlight.png'); rotI=imrotate(I,90,'crop'); fig1=imshow(rotI); BW1=rgb2gray(rotI); BW2=im2bw(BW1,0.95); BW3=double(BW2); BW=edge(BW3,'canny'); figure, imshow(BW); [H,theta,rho]=hough(BW); figure,imshow(imadjust(mat2gray(H)),[],'XData',theta,'YData',rho,... 'InitialMagnification','fit'); xlabel('\theta(degrees)'),ylabel('\rho'); axis on, axis normal, hold on; colormap(hot) P=houghpeaks(H,50,'threshold',ceil(0.3*max(H(:)))); x=theta(P(:,2)); y=rho(P(:,1)); plot(x,y,'s','color','black'); —————————————————————————————————————— ------------------------------------------------------------------------------------------------ lines=houghlines(BW,theta,rho,P,'FillGap',18,'MinLength',180); figure, imshow(rotI), hold on max_len=0; for k=1:length(lines) xy=[lines(k).point1;lines(k).point2]; plot(xy(:,1),xy(:,2),'LineWidth',2,'color','green'); plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow'); plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red'); len=norm(lines(k).point1-lines(k).point2); if (len>max_len) max_len=len; xy_long=xy; end end plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','red'); —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 24 第四问: clc clear length=[2.424024024 2.391081081 2.357837838 2.324834835 2.315825826 2.276546547 2.246546547 2.207417417 2.176066066 2.162462462 2.121861862 2.105255255 2.09039039 2.054144144 2.027087087 2.003063063 1.966996997 1.94 1.912942943 1.900900901 1.876876877 1.858858859]; Length=length; sOLen=sum(Length); result=[]; t=8.90166666666667:0.0333333333333332:9.60166666666667; min=2; bc=0.1; lon=105:bc:150; beginl=-30*ones(1,size(lon,2)); endl=70*ones(1,size(lon,2)); n=194; delta=23.45*sin(2*pi*(284+n)/365); for i=1:size(lon,2) ts=t-(120-lon(1,i))/15; —————————————————————————————————————— ------------------------------------------------------------------------------------------------ omega=15*(ts-12); beginl(1,i); endl(1,i); for lat=beginl(1,i):bc:endl(1,i) phi=lat; h=asin(sin(phi*pi/180)*sin(delta*pi/180)+cos(phi*pi/180)*cos(delta*pi/ 180)*cos(omega*pi/180)); l=2./tan(h); if sum(abs(l-length))<min min=abs(l-length); ganH=mean(Length.*tan(h)); result=[lon(l,i) lat ganH abs(l-length) n]; end clear h; end i end plot(result(:,1),result(:,2),'r*') 截取投影的程序: 25 clc —————————————————————————————————————— ------------------------------------------------------------------------------------------------ clear close all syms xy_long I=imread('sunlight.png'); rotI=imrotate(I,90,'crop'); fig1=imshow(rotI); BW1=rgb2gray(rotI); BW2=im2bw(BW1,0.95); BW3=double(BW2); BW=edge(BW3,'canny'); figure, imshow(BW); [H,theta,rho]=hough(BW); figure,imshow(imadjust(mat2gray(H)),[],'XData',theta,'YData',rho,... 'InitialMagnification','fit'); xlabel('\theta(degrees)'),ylabel('\rho'); axis on, axis normal, hold on; colormap(hot) P=houghpeaks(H,50,'threshold',ceil(0.3*max(H(:)))); x=theta(P(:,2)); y=rho(P(:,1)); plot(x,y,'s','color','black'); —————————————————————————————————————— ------------------------------------------------------------------------------------------------ lines=houghlines(BW,theta,rho,P,'FillGap',18,'MinLength',180); figure, imshow(rotI), hold on max_len=0; for k=1:length(lines) xy=[lines(k).point1;lines(k).point2]; plot(xy(:,1),xy(:,2),'LineWidth',2,'color','green'); plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow'); plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red'); len=norm(lines(k).point1-lines(k).point2); if (len>max_len) max_len=len; xy_long=xy; end end 26 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color', 'red'); 第四问视频中影长的值: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 —————————————————————————————————————— ------------------------------------------------------------------------------------------------ 20 21 22 2.424024024 2.391081081 2.357837838 2.324834835 2.315825826 2.276546547 2.246546547 2.207417417 2.176066066 2.162462462 2.121861862 2.105255255 2.09039039 2.045757681 2.017128007 1.988498333 1.95986866 1.931238986 1.902609312 1.873979639 1.845349965 1.816720291 27 ——————————————————————————————————————