2015年高教社杯数学建模论文

1基于最小二乘拟合和蒙特卡洛法实现太阳影子定位摘要本文基于太阳高度角的相关原理，建立影子长度与直杆长度、纬度、积日和时角的函数关系，使用最小二乘拟合法、蒙特卡洛法、分类筛选法、误差修正法来实现太阳影子定位。对于问题一，建立多元分析模型来分析影子长度与纬度、积日、地方时间之间的变化规律。固定其中两个变量，分析第三个变量与影子长度的变化规律并画出天安门广场直杆影子长度的变化曲线。对于问题二，建立最小二乘拟合参数模型，使用附件1的数据，得出直杆的地点可能为：北纬19度16分、东经108度39分（中国海南省），南纬3度39分、东经102度41分（南苏门答腊）。对于问题三，建立蒙特卡洛拟合参数模型，得出大量拟合参数，通过拟合优度检验和分类筛选，得出附件2中直杆的可能地点与日期为：北纬34度57分、东经118度53分（江苏连云港）、12月20日或12月22日；北纬37度14分、东经119度31分（渤海）、12月21日；北纬33度14分，东经115度49分（安徽阜阳）、12月20日或12月22日。附件3中可能的地点与日期为：北纬46度25分、东经119度39分（内蒙古自治区）、12月20日或12月22日；北纬36度57分、东经118度26分（山东临沂）、12月21日；北纬33度14分，东经118度4分（江苏宿迁）、12月20日或12月22日。对于问题四，首先使用相似原理计算太阳影子长度，结合问题二模型得出若干可能地点为：北纬39度31分、东经109度04分（内蒙古鄂尔多斯市）。然后建立误差修正模型对结果进行修正。得到修正后的若干可能地点为：北纬49度36分、东经104度22分（蒙古色楞格）；北纬45度36分、东经110度（蒙古赛音山达）；北纬48度21分、东经110度（蒙古鄂嫩）。若日期未知，可应用问题三中的蒙特卡洛拟合模型得出若干可能地点与日期。关键词：太阳高度角最小二乘拟合蒙特卡洛法拟合优度误差修正模型2一、问题重述如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。4．附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？二、问题分析针对问题一，建立多元分析模型来分析影子长度L与纬度、积日V、地方时间T之间的变化规律。太阳高度角与赤纬角、纬度、时角t（t=|18015|T）之间有如下关系式[1]；sin()sin()sin()+cos()cos()cos(),t（22(284V)23.45sin365）。由上可得L与、V、T之间的函数关系式为：201(,,)((,,)sin())(,,)LfTVLfTVfTV，固定其中两个变量可得第三个变量与太阳影子长度L的函数关系式，分析第三个变量与影子长度L的变化规律并画出天安门广场直杆影子长度的变化曲线。3针对问题二，建立最小二乘拟合参数模型。附件1中的时间是北京时间而并非杆子所处的当地时间，设北京时间为0T，0TT-a（a为时差）。在关系式201(,,)(,,)LfTVLfTV中未知参数为0La，，，把影子长度L与北京时间0T导入matlab中，用最小二乘法拟合参数,运用拟合优度newR（newR越接近1，拟合程度越好）检验得出拟合参数0La，，。由时差a可算出测量地点的经度(120)4a，再结合纬度可得到直杆的地点。运用附件1中的数据结合该模型可得出若干可能的地点。针对问题三，建立蒙特卡洛随机拟合参数模型。由于变量增加，拟合出来的参数不唯一，故本文采取蒙特卡洛随机拟合模型产生大量拟合参数，根据参数分布情况并结合拟合优度挑选出符合实际情况的参数。在函数关系式201(,,)(,,)LfTVLfTV中，地方时间0TT-a（0T为北京时间），分析可知在关系式中未知参数为0,a,L,V，把L与0T导入matlab中，运用最小二乘法结合拟合优度检验得出大量拟合参数0,a,L,V，，根据参数分布情况，挑选有代表性的参数0,a,L,V，然后由时差a可算出测量地点的经度，再结合纬度可得到直杆的地点。运用附件2、3中的数据结合该模型可得出若干可能的地点。针对问题四，首先采用相似原理求出影子近似长度，运用问题2中的最小二乘拟合模型求出可能的地点，然后建立误差修正模型修正结果。本问已知杆长L为2m,由相似原理，有00LlLl（l与0l分别为影像中的杆长和影长），使用截屏编辑工具对附件4中的视频进行处理，每隔4分钟，输出一次l和0l，计算得到近似影子长度0L。根据视频中的北京时间0T可以得到一组太阳影子长度0L和0T的数据，使用最小二乘拟合模型，可求出拍摄该视频若干可能的地点。由于实际影子长度与测量出的影子长度之间存在一定误差（存在偏角），此处使用修正模型来减小结果误差。同理，若该视频没有日期，可先求测量出的太阳影子长度0L与4北京时间0T的数据，再使用蒙特卡洛拟合参数模型求出拍摄视频若干可能的地点与日期。三、模型假设1、假设杆子与地面垂直2、假设地球是圆的3、不考虑太阳光线在穿过大气层时的折射、太阳的视角、高山阻挡、海拔高度等因素的影响，照射到地球的太阳光可以看成是平行光线四、符号说明符号含义单位L太阳影子长度米纬度度V积日天a时差分钟太阳高度角度0L杆长米赤纬角度t时角度T地方时时五、模型的建立与求解5.1问题一模型建立与求解-多元分析模型5.1.1模型建立本文建立了多元分析模型来分析影子长度L与纬度、积日V、地方时间T之间的变化规律。查阅文献【1】知太阳影子长度L与太阳高度角有关，它们有如下关系:0tan()LL,其中0L为杆长，即有0tan()LL，具体如下图所示：5太阳高度角与赤纬角、纬度、时角t（t=|18015|T）之间有如下关系式:sin()sin()sin()+cos()cos()cos(),t赤纬角的计算公式如下：2(284V)23.45sin365,V为积日，即日期在年内的顺序号，例如，1月1日其积日为1，平年12月31日的积日为365。将与纬度代入太阳高度角计算公式，推导出sin()(,,)fTV，可得影子长度L与纬度、积日V、时间T之间的函数关系式为：201(,,)(,,)LfTVLfTV，固定其中两个变量可得第三个变量与太阳影子长度L的函数关系式，分析第三个变量与L的变化规律并画出天安门广场3米长直杆影子长度在2015年10月22日北京时间9:00-15:00之间的变化曲线。5.1.2模型求解：求出L与各个参数的函数关系式如下：201(sin()sin()cos()cos()cos(|18015|))2(284V),23.45sinsin()sin()cos()cos()cos(|18015|)365LTLT固定其中两个参数，分析第三个参数与太阳影子长度之间的关系。（1）取定时间T，纬度，由常识可知回归线左右太阳影子的情况略有不同，此处以北半球为例，设杆长03,L取12T，=12度、23度26分、40度来分析太阳影子长度L与积日V之间的变化规律,使用matlab软件(程序见附录)0LL6画出L与V的图像如下：010020030040001234567积日影子长度N23°26'N40°N12°结果分析如下：由图像可知，在北纬40度和北纬23度26分，正午12点，随着积日V的增加，影子长度先变短后变长，在170V左右，北纬40度区域太阳影子长度达到最小约1米，此时北纬40度地区处于夏季，影子长度最短，符合实际情况。170V左右，太阳直射北回归线，在正午12点，太阳影子长度为0。由常识可知在正午12点，因为北纬12度处在南北回归线之间，故它的太阳影子长度随着积日V的增加，在正午12点会经历两次太阳影子长度为0，上图很好的验证了这一点。（2）取定时间T,积日V，使用matlab软件画出L与的图像如下：此处，T=12，即为正午12点，积日V=80、180，横坐标纬度已换算成弧度。杆长0L=3.700.511.505101520253035404550北半球纬度影子长度积日80积日180结果分析如下：由图像（这里纬度的单位已换算成弧度）可知，积日V=80时，在正午12点，随着纬度的增加，太阳影子长度从0开始增加，当达到90°时，太阳影子达到最长。由常识可知V=80左右，太阳直射赤道，随着纬度增加，太阳影子长度逐渐增加，上图符合实际情况。积日V=180时，在正午12点，随着纬度的增加，太阳影子长度先减小，后增加。在V=180左右，太阳直射北回归线，故当纬度从0开始增加时，太阳影子长度先减小，在北回归线处为0，后增加，在90°处达到最大，上图符合实际情况。（3）取定积日V,纬度.此处积日V=295，纬度为北纬39度54分26秒、东经116度23分29秒（天安门广场），经过计算得到太阳影子长度L与时间T之间的函数关系式如下：231(0.0450.7688cos(|15180|))0.0450.7688cos(|15180|)TLT，使用matlab软件（程序见附录）画出天安门广场2015年10月22日北京时间9:00-15:00之间太阳影子长度L与时间T的变化曲线如下：8910111213141522.533.544.5时间影子长度结果分析如下：由上图分析可知，固定积日V和纬度,随着时间的增加，太阳影子长度先减小，在12点左右达到最小2.2米，随后随着时间的增加而增加，与实际情况相符。5.2问题二模型-最小二乘拟合参数模型：5.2.1模型求解由第一问，有如下函数关系式：sin()sin()sin()+cos()cos()cos(),t2(284V)23.45sin,365sin()(,,)fTV，201(,,)(,,)LfTVLfTV附件1中所给的时间是北京时间而并非杆子所处的当地时间，，由测量日期可知积日V，设北京时间为0T，且有0TT-a。分析可知在关系式201(,,)(,,)LfTVLfTV中未知参数为0La，，，将附件1中影子顶点坐标换算成影子长度L，把L与0T导入matlab中，运用最小二乘法求出拟合参数90La，，，由时差a可算出测量地点的经度（1204a），由经纬度,可得到直杆的地点，最后进行拟合优度分析。其原理如下：对非线性方程，计算残差平方和2()cQLL（CL为用拟合参数计算出的影子长度）以及2RL。拟合优度newR=1QR，newR越接近1，说明拟合参数0La，，越可靠，newR越小，说明拟合参数越不可靠。根据newR的大小判断拟合出的参数是否可靠。5.2.2模型求解：使用matlab软件计算得到参数0La，，如下（程序见附录）：参数a0L拟合值1N'1916-45.4min2.0365m拟合值2'S339-69.27min2.2527m拟合优度结果分析：使用matlab（程序见