基于非线性最小二乘法的太阳影子定位研究

1、一、 摘要本文研究的是利用太阳影子随时间变化的顶点坐标数据或照片进行定位的问题。根据地理知识建立数学模型，研究影长关于多个参数的变化规律，并运用非线性最小二乘法进行拟合。对已知影长地点的经纬度，进行参数最优估计，求解到可能的地点。首先根据地理知识可知，太阳下影子长度只由杆长和该地的太阳高度角决定。而太阳高度角由该时刻的太阳赤纬角、时角和太阳的地理纬度决定。所以根据上述几个参数建立了影子长度变化模型。另一方面，为方便比较建立以旗杆为原点，正北方向为y轴正方向，正东方向为x轴正方向的坐标系。引入太阳方位角概念，太阳方位角与太阳高度角、太阳赤纬角和当地纬度有关。并根据影子长度变化模型的结果与太阳方位

2、角，得到备选城市影子顶点坐标。通过题目所给图片影子顶点位置圆心的像素点位置，获得图片影子顶点坐标。建立了刻画某地一天中影子顶点随时间变化轨迹的数学模型。进而通过比较原始数据与备选城市数据二阶差商比值的方差，判断其离散程度，从而得到最有可能是其拍摄未知的城市。模型检验时，做关于纬度和杆长两个参数的灵敏度分析。考虑将附件中的原始数据影长做±10%的变换。二 、问题的背景与重述及分析 2.1问题的背景早在15世纪时，定位技术就已经随着海洋探索的开始而产生。随着社会和科技的不断发展，我们对定位的需求已不再局限于航海、航空等领域，对于地球上的精确坐标定位已逐渐成为人们关注的热点问题。对于地球表

3、面经纬度的精确定位，可利用变化的太阳影子来进行分析，其作为一种直观简便的定位技术，已受到广泛关注。2.2问题的重述附件是2016年4月18日在某地数码相机拍摄的旗杆顶点下午14点42分开始每隔3分钟的影子的位置，旗杆高约2米，垂直地面。请建立数学模型，确定该拍摄地最可能是昆明、南京、西安、北京或哈尔滨的那个城市？2.3问题的分析 由于问题所给图片中没有出现原本旗杆的位置以及图片比例尺，所以无法判断具体的影长。故不能采用一般直接比较的方法。所以在通过像素点获得了影子顶点坐标后采用了通过拟合其一阶差分曲线，比较二阶差分与所给其它城市的影长顶点坐标的二阶差分离散程度的方法。三、模型假设1.假设一天中

4、的太阳赤纬角保持不变；2.假设大气层对太阳光的折射率保持不变；3.假设影子长度和角度与该点的海拔无关；4.假设地面平整；5.假设杆子长度为严格2米（题目中为约为2米）；6.仅考虑照片拍摄角度为俯拍；7.由于所给题目中未提及所给时间，所以我们采用了分类讨论的方法，在第一个模型中假设所给时间为当地准确时间，在第二个模型中假设为当天的北京时间。四符号说明符号表示含义单位h太阳高度角radH考虑大气折射的太阳高度角radL旗杆长度md影子长度m太阳赤纬角radST真太阳时hMT平太阳时h太阳时角radN积日天地理纬度度地理经度度n大气折射率AS太阳方位角度日角度X1图片影子顶点位置横坐标Y1图片影子顶

5、点位置纵坐标X2各地点影子定点位置横坐标（旗杆为原点）Y2各地点影子定点位置纵坐标（旗杆为原点）五模型的建立与求解5.1影子长度变化模型5.1.1影子长度1图一：影子长度变化模型本图即为太阳位置，旗杆高度与影长关系的模拟图。图中旗杆高为h，影长为L，太阳高度角为H。根据图中所示关系，可以得到数学模型表达式：d=Ltanh (1)5.1.2太阳时角太阳时角，是指日面中心的时角，即从观测点天球子午圈沿天赤道量至太阳所在时圈的角距离。2单位时间地球自转的角度定义为时角，规定正午时角为0，上午时角为负值,下午时角为正值。对于平太阳时来说，地球自转一周360度,对应的时间为24小时,即每小时相应的时角为

6、15度,每4分钟的时角为1度。因为地球自转周期不为准确的24小时，所以为保证时角的准确性，采用修正过的真太阳时。（时差表见附录一）2016年4月18日平太阳时转换为真太阳时的计算公式为（模型一）：ST=MT+16” （2）2016年4月18日平太阳时转换为真太阳时的计算公式为（模型二）： ST=MT-4*（120-）+16” （2）时角的计算公式为： =（ST-12）*15 （3）5.1.3日角日角是决定太阳赤纬的一个因素，其计算公式为： =2(N-79.6764-0.2422×y-1985+INTy-19854)365.2422 5.1.4太阳赤纬太阳赤纬，是地球赤道平面与太阳和地

7、球中心的连线之间的夹角。赤纬角以年为周期，在+23 °26与-23 °26的范围内移动，成为季节的标志。3太阳赤纬的计算公式为：=0.3723+23.2567*sin+0.1149*sin2-0.1712*sin3- 0.758*cos+0.3656*cos2+0.201*cos3 （4）其中，4月18日时，N=109。 经计算，4月18日所对应的赤纬角为11.16332°。5.1.5太阳高度角 太阳高度角是指某地太阳光线与通过该地与地心相连的地表切面的夹角。4太阳高度角计算公式为： h=sin-1(sinsin+coscoscos) （5）综上，可得关于影子长度

8、的数学模型建立如下（一）：&d=Ltanh &ST=MT+16” &=ST-12×15&=0.3723+23.2567*sin+0.1149*sin2-0.1712*sin3-0.758*& cos+0.3656*cos2+0.201*cos3 &h=sin-1(sinsin+coscoscos)可得关于影子长度的数学模型建立如下（二）：&d=Ltanh &ST=MT-4×（120-）+16”&=ST-12×15&=0.3723+23.2567*sin+0.1149*sin2-0.17

9、12*sin3-& 0.758*cos+0.3656*cos2+0.201*cos3 &h=sin-1(sinsin+coscoscos)5.2影子顶端坐标模型为了了解所给数据为所给备选城市中的哪一个，需要首先了解备选城市在同一时段中其影子顶端坐标。在影子长度变化模型中，首先建立了以、坐标系描绘的数学模型：d=Ltan(sin-1(sinsin+coscoscos(MT-4×120-+1'6"-12)×15°(=11.16332°)其中为影长，为太阳方位角。之后以旗杆位置为原点，正东方向为x轴正方向，正北方向为y轴正方向

10、，将建立的影子长度变化模型坐标转化为x,y坐标系坐标5.2.1太阳方位角太阳方位角即太阳所在的方位，指太阳光线在地平面上的投影与当地经线的夹角，可近似地看作是竖立在地面上的直线在阳光下的阴影与正南方的夹角。5AS=cos-1(sinhsin-sincoshcos)5.2.2影子顶端横坐标X=dcosAS5.2.3影子顶端纵坐标Y=dsinAS图二：北京时各城市影子顶点位置坐标图三：地方时各城市影子顶点位置坐标南京北京西安哈尔滨昆明XYXYXYXYXY1.6799320.6380061.8463040.945861.1965070.7863762.3995961.1312650.8612010.

11、4399441.7225890.6321831.8918360.9390031.2318990.7825792.4553031.1209890.8919780.437471.7660270.6261551.9381810.9319111.2677370.7786452.5122051.1103750.923070.4348951.8102820.6199141.9853770.9245751.3040410.7745692.5703261.0993520.9544910.4322161.8553910.6134532.0334590.9169851.3408290.7703482.629721

12、1.0879610.9862560.4294321.9013910.6067632.0824670.9091321.3781230.7659772.6904511.0761651.0183810.4265391.9483250.5998532.1324430.9010061.4159450.7614522.7525771.0639481.0508820.4235351.9962350.5926592.1834310.8925951.4543160.7567682.8161681.051291.0833750.4204172.0451660.5852252.2354760.8838891.493

13、2590.7519192.8812931.0381721.1170790.417182.0951670.5775222.2886270.8748751.5328010.7469012.9480281.0245731.1508110.4138232.146290.5695382.3429350.865541.5729660.741703.0164531.010471.1849880.4103422.1985870.5612622.3984560.8558711.6137810.7363343.0866530.9958391.2196330.4067332.2521170.5526812.4552

14、460.8458531.6552740.7307723.1587190.9806551.2547640.4029922.3069390.5437792.5133670.835471.6974750.7250173.2327480.9648891.2904040.3991152.3631190.5345442.5728840.8247071.7404140.7190623.3088420.9485131.3265730.3950982.4207260.5249592.6338640.8135761.7841240.7128993.3871130.9314951.3632970.3909362.4

15、798320.5150072.6963820.8019671.828640.706523.4676780.9138021.4005980.3866252.5405160.5046712.7605160.7899521.8739970.6999173.5506640.8953981.4385020.382162.602680.4939312.826470.777481.9202330.6930823.6362070.8762441.4770370.3775352.6669540.4827672.8939660.7645261.9673870.6860063.7244530.8562981.516

16、2290.3727442.7328930.4711582.9634650.7510692.0155030.6786793.8155590.8355151.556110.367783表一：北京时各城市影子顶点位置坐标图5.3比较模型5.3.1比较原理由于获取的原始数据为本题中所给图片中影子顶点的像素点坐标，并且旗杆位置未知，图片可能进行了比例尺变换，以及平移旋转变换，故无法采用其数据直接进行欧式几何分析。进而无法确定原始图片在本文所建模型中的具体坐标。基于此，本文对于原始数据的一阶导数变化率进行衡量，通过对其一阶差商和二阶差商进行分析，来确定与题目所给五个城市所得数学模型的近似程度，从而判断图片

17、拍摄地点最有可能所在城市。原始数据与五个目标城市的二阶差商应存在如下关系：d2ydx2=d2yidxi2(i=1,2,3,4,5)最后，通过求解各城市不同时间所对应的的数列的方差，判断其离散程度，方差越小，离散程度越小，即照片拍摄地点越接近所对应城市。5.3.2原数据处理（1）根据图片获得其像素点坐标并绘制坐标图:XYXY83472331632633125614545314280128642993041118328232893103266353751242493785814423240340166215429221881984555210180表二：原始数据影子顶点坐标（像素点）图四：原始数据

18、影子顶点坐标所绘图形及直线对比（2）根据上述模型计算，结果如下：根据（1）中所获得的图片影子顶点坐标数据，求其各点的一阶差商。根据一阶差商的结果我们发现直接求原数据坐标的一阶差商有较大的波动，无法使用。（如下图）图五：北京时原始数据一阶差商描点图为此，我们选择对于这些点进行二阶多项式拟合，然后再对拟合曲线进行对比。（结果如下图）图六：北京时原始数据一阶差商描点图及拟合曲线得到了拟合曲线的公式：y=4×10-7x2-0.0007x+0.8844之后求得此拟合曲线的二阶差商。原数据二阶差商-0.00069-0.0006-0.0005-0.00038-0.00068-0.00058-0.0

19、0048-0.00066-0.00057-0.00046-0.00065-0.00055-0.00044-0.00063-0.00053-0.00042-0.00062-0.00051-0.0004表三：北京时原数据二阶差商5.3.3选项城市数据处理由于选项的计算为理论值，所以其一阶差分曲线不存在波动的情况（见下图），可以直接使用。利用影子顶点坐标可直接求其二阶差分。图七：北京时备选城市影子顶点坐标一阶差商图八：地方时备选城市影子顶点坐标一阶差商南京北京西安哈尔滨 昆明-0.05288-0.05354-0.07067-0.04483-0.07904-0.05153-0.05223-0.0694

20、3-0.04357-0.07789-0.05017-0.05092-0.06818-0.04231-0.07671-0.04882-0.04961-0.06691-0.04106-0.07552-0.04748-0.0483-0.06563-0.03982-0.0743-0.04613-0.047-0.06435-0.03859-0.07306-0.04479-0.04569-0.06305-0.03736-0.07181-0.04346-0.0444-0.06174-0.03614-0.07053-0.04214-0.04311-0.06042-0.03494-0.06925-0.04082

21、-0.04182-0.0591-0.03375-0.06794-0.03951-0.04054-0.05777-0.03256-0.06662-0.03821-0.03927-0.05644-0.0314-0.06529-0.03693-0.03801-0.0551-0.03024-0.06395-0.03565-0.03676-0.05376-0.0291-0.0626-0.03439-0.03552-0.05242-0.02798-0.06124-0.03314-0.03429-0.05108-0.02687-0.05987-0.03191-0.03308-0.04974-0.02578-

22、0.0585-0.03069-0.03188-0.0484-0.0247-0.05712-0.02949-0.03069-0.04707-0.02365-0.05574表四：北京时各城市影子顶点坐标二阶差商5.3.4各城市二阶差商比值方差南京北京西安哈尔滨昆明1.0958852.19033546.198561.122924100.2249表5：北京时各城市二阶差商比值方差南京北京西安哈尔滨昆明0.8268032.1903350.8873896.3559751.79557表6：地方时各城市影子顶点二阶差商比值方差据此可知，最有可能为南京。5.5灵敏度分析图九：太阳高度角对纬度的偏导图像在求最优化

23、解的模型中，我们经常利用灵敏度分析来研究原始数据不准确或发生变化时最优解的稳定性。通过灵敏度分析，我们可以发现哪些参数对模型具有较大的影响。对本题的灵敏度分析主要考虑纬度及杆长的影响。上图是将经度以及杆长固定，对纬度进行±10%的变化探究纬度对于太阳高度角的影响。在32°的理论值区间内，对太阳高度角的偏导近似为一常数，故而两者在此区间内呈线性变化，斜率约为0.8.可见稳定性较差，随纬度的变化，影长的变化性较大。即纬度的变化对各个参数的最优解的稳定性均有较大影响。杆长与影长呈线性关系，故杆长对影长的影响较大。六、模型的评价与改进6.1模型的评价6.1.1模型的优点1.本文考虑

24、了平太阳时与真太阳时对于太阳时角的影响，进行了时差的修正，使结果更为精确。2.本模型考虑了题目所给时间为地方准确时间与北京时间（经度120°）两种情况进行分析，使模型结果更为全面。6.1.2模型的缺点1.由于影长曲线本身的复杂性和波动性，所以采用了二次函数进行拟合，导致本模型存在一定的误差。2.没有考虑到地面不平整以及地点本身海拔高度等因素。6.2模型的改进6.2.1太阳高度角的大气折射修正图九：大气折射率对于太阳高度角的影响如图所示，太阳光在射入大气层时会产生折射。故在计算时增加考虑大气折射可以提高该模型的准确度。由图可知：n=sinRsinr=sin(90°-h)sin

25、(90°-H)根据文献查得公式：n=2.87774×10-9P1+10-10P60.1-0.972T1+0.003661T-2.58×10-9R×e0.059T+16其中P为大气压强，T为温度，R为相对湿度。参考文献：123456刘海宁，李真，张国雄.以相对湿度为参量的大气折射率计算公式J.应用激光,2000,35(5):798-802附表一真平太阳时时差表1月01日-3分9秒1月02日-3分38秒1月03日-4分6秒1月04日-4分33秒1月05日-5分1秒1月06日-5分27秒1月07日-5分54秒1月08日-6分20秒1月09日-6分45秒1月10

26、日-7分10秒1月11日-7分35秒1月12日-7分59秒1月13日-8分22秒1月14日-8分45秒1月15日-9分7秒1月16日-9分28秒1月17日-9分49秒1月18日-10分9秒1月19日-10分28秒1月20日-10分47秒1月21日-11分5秒1月22日-11分22秒1月23日-11分38秒1月24日-11分54秒1月25日-12分8秒1月26日-12分22秒1月27日-12分35秒1月28日-12分59秒1月29日-13分10秒1月30日-13分19秒1月31日-13分37秒2月01日-13分44秒2月02日-13分50秒2月03日-13分56秒2月04日-14分1秒2月05日

27、-14分5秒2月06日-14分9秒2月07日-14分11秒2月08日-14分13秒2月09日-14分14秒2月10日-14分15秒2月11日-14分14秒2月12日-14分13秒2月13日-14分11秒2月14日-14分8秒2月15日-14分5秒2月16日-14分1秒2月17日-13分56秒2月18日-13分51秒2月19日-13分44秒2月20日-13分38秒2月21日-13分30秒2月22日-13分22秒2月23日-13分13秒2月24日-11分4秒2月25日-12分54秒2月26日-12分43秒2月27日-12分32秒2月28日-12分21秒2月29日-12分8秒3月01日-11分56秒

28、3月02日-11分43秒3月03日-11分29秒3月04日-11分15秒3月05日-11分1秒3月06日-10分47秒3月07日-10分32秒3月08日-10分16秒3月09日-10分1秒3月10日-9分45秒3月11日-9分28秒3月12日-9分12秒3月13日-8分55秒3月14日-8分38秒3月15日-8分21秒3月16日-8分4秒3月17日-7分46秒3月18日-7分29秒3月19日-7分11秒3月20日-6分53秒3月21日-6分35秒3月22日-6分17秒3月23日-5分58秒3月24日-5分40秒3月25日-5分22秒3月26日-5分4秒3月27日-4分45秒3月28日-4分27

29、秒3月29日-4分9秒3月30日-3分51秒3月31日-3分33秒 4月01日-3分16秒4月02日-2分58秒4月03日-2分41秒4月04日-2分24秒4月05日-2分7秒4月06日-1分50秒4月07日-1分33秒4月08日-1分17秒4月09日-1分1秒4月10日+0分46秒4月11日+0分30秒4月12日+0分16秒4月13日+0分1秒4月14日+0分13秒4月15日+0分27秒4月16日+0分41秒4月17日+0分54秒4月18日+1分6秒4月19日+1分19秒4月20日+1分31秒4月21日+1分42秒4月22日+1分53秒4月23日+2分4秒4月24日+2分14秒4月25日+2

30、分23秒4月26日+2分33秒4月27日+2分41秒4月28日+2分49秒4月29日+2分57秒4月30日+3分4秒5月01日+1分10秒5月02日+3分16秒5月03日+3分21秒5月04日+3分26秒5月05日+3分30秒5月06日+3分37秒5月07日+3分36秒5月08日+3分39秒5月09日+3分40秒5月10日+3分42秒5月11日+3分42秒5月12日+3分42秒5月13日+3分42秒5月14日+3分41秒5月15日+3分39秒5月16日+3分37秒5月17日+3分34秒5月18日+3分31秒5月19日+3分27秒5月20日+3分23秒5月21日+3分18秒5月22日+3分13秒

31、5月23日+3分7秒5月24日+3分1秒5月25日+2分54秒5月26日+2分47秒5月27日+2分39秒5月28日+2分31秒5月29日+2分22秒5月30日+2分13秒5月31日+2分4秒6月01日+1分54秒6月02日+1分44秒6月03日+1分34秒6月04日+1分23秒6月05日+1分12秒6月06日+1分0秒6月07日+0分48秒6月08日+0分36秒6月09日+0分24秒6月10日+0分12秒6月11日+0分1秒6月12日+0分14秒6月13日+0分39秒6月14日+0分52秒6月15日-1分5秒6月16日-1分18秒6月17日-1分31秒6月18日-1分45秒6月19日-1分5

32、7秒6月20日-2分10秒真平太阳时时差表6月21日-2分23秒6月22日-2分36秒6月23日-2分48秒6月24日-3分1秒6月25日-3分13秒6月26日-3分25秒6月27日-3分37秒6月28日-3分49秒6月29日-4分0秒6月30日-4分11秒 7月01日-4分22秒7月02日-4分33秒7月03日-4分43秒7月04日-4分53秒7月05日-5分2秒7月06日-5分11秒7月07日-5分20秒7月08日-5分28秒7月09日-5分36秒7月10日-5分43秒7月11日-5分50秒7月12日-5分56秒7月13日-6分2秒7月14日-6分8秒7月15日-6分12秒7月16日-6分

33、16秒7月17日-6分20秒7月18日-6分23秒7月19日-6分25秒7月20日-6分27秒7月21日-6分29秒7月22日-6分29秒7月23日-6分29秒7月24日-6分29秒7月25日-6分28秒7月26日-6分26秒7月27日-6分24秒7月28日-6分21秒7月29日-6分17秒7月30日-6分13秒7月31日-6分8秒8月01日-6分3秒8月02日-5分57秒8月03日-5分51秒8月04日-5分44秒8月05日-5分36秒8月06日-5分28秒8月07日-5分19秒8月08日-5分10秒8月09日-5分0秒8月10日-4分50秒8月11日-4分39秒8月12日-4分27秒8月1

34、3日-4分15秒8月14日-4分2秒8月15日-3分49秒8月16日-3分36秒8月17日-3分21秒8月18日-3分7秒8月19日-2分51秒8月20日-2分36秒8月21日-2分20秒8月22日-2分3秒8月23日-1分47秒8月24日-1分29秒8月25日-1分12秒8月26日+0分54秒8月27日+0分35秒8月28日+0分17秒8月29日+0分2秒8月30日+0分21秒8月31日+0分41秒9月01日+1分0秒9月02日+1分20秒9月03日+1分40秒9月04日+2分1秒9月05日+2分21秒9月06日+2分42秒9月07日+3分3秒9月08日+3分3秒9月09日+3分24秒9月1

35、0日+3分45秒9月11日+4分6秒9月12日+4分27秒9月13日+4分48秒9月14日+5分10秒9月15日+5分31秒9月16日+5分53秒9月17日+6分14秒9月18日+6分35秒9月19日+6分57秒9月20日+7分18秒9月21日+7分39秒9月22日+8分0秒9月23日+8分21秒9月24日+8分42秒9月25日+9分2秒9月26日+9分22秒9月27日+9分42秒9月28日+10分2秒9月29日+10分21秒9月30日+10分40秒 10月01日+10分59秒10月02日+11分18秒10月03日+11分36秒10月04日+11分36秒10月05日+11分53秒10月06日+12分11秒10月07日+12分28秒10月08日+12分44秒10月09日+12分60