To be continue … …

主要参考文献:[2015年国赛高教杯奖A题]电子科技大学-太阳影子定位的多目标优化模型

最近用java写了个求影长的程序,原地爆炸

文章目录

  • 一、摘要
  • 二、问题重述
    • 1. 问题背景
    • 2. 需要解决的问题
  • 三、问题分析
    • 3.1 问题一的分析
    • 3.2 问题二的分析
    • 3.3 问题三的分析
    • 3.4 问题四的分析
  • 四、模型假设
  • 五、定义与符号说明
    • 5.1 符号说明
    • 5.2 名词解释
  • 六、模型的建立与求解
    • 6.1 模型一:基于太阳赤纬与时差计算的结构方程模型——问题一的求解
      • 6.1.1 建模准备——影子长度计算公式
      • 6.1.2 影子长度关于各个参数变化规律
        • (1) 纬度φ\varphiφ与日数NNN一定,影子长度ddd随时间t′t't′变化规律
        • (2) 纬度φ\varphiφ与时间t′t't′一定,影子长度ddd随天数N的变化规律
        • (3) 日数NNN与时间t′t't′一定,影子长度ddd随天数NNN变化规律
        • (4) 影子长度在指定日期时间下的天安门广场9:00-15:00的变化曲线图
    • 6.2 模型二:基与最小二乘法模型——问题二求解
      • 6.2.1 建模准备
      • 6.2.2 数据的预处理
      • 6.2.3 问题二的求解——基于半搜索算法
        • 1. 模型阶段一:数据点的选取
        • 2. 模型阶段二:半搜索算法思路
        • 3. 方法二:循环遍历搜索法
      • 6.2.4 模型检验与结果分析
    • 6.3 模型三: 基与分治思想的循环遍历法——问题三的求解
      • 6.3.1 建模准备
      • 6.3.2 多目标优化模型的建立
      • 6.3.3 模型的求解
        • 1. 模型阶段一:数据预处理
        • 2. 模型阶段二:基于分治思想的循环遍历法
        • 3. 模型阶段三:求解结果
    • 6.4 模型四:基与坐标变换的图形射影模型——问题四的求解
      • 6.4.1 建模准备
      • 6.4.2 第 1 小问——拍摄日期已知
      • 6.4.3 第 2 小问——拍摄日期未知
  • 七、模型的评价与推广
    • 7.1 模型的评价
      • 7.1.1 模型的优点
      • 7.1.2 模型的缺点与不足
    • 7.2 模型的改进
  • 参考文献
  • 附录

一、摘要

本文是一个对太阳影子定位的研究模型。通过对太阳赤纬角、时区分划的计算、研究以及直杆在不同时角下影长的变化,分别建立了基于太阳影子定位的结构方程模型、最小二乘模型和基于分治思想的循环搜索、多目标优化等模型。

首先,本文针对问题一的指定地理位置和日期的直杆影长的变化问题,对影子定位方程进行了简化,利用几何知识建立了太阳影子定位的结构方程。得到影长与纬度日数时间三个参数之间的解析关系式。并以天安门广场的直杆为例,在另两个参数一定的情况下,利用MATLAB软件得到影长随参数的变化规律。影长关于时间的变化为先减小后增大,并在正午12点达到最小值;影长关于日数的变化为先减小后增大,夏季影长最短,冬季影长最长;影长关于纬度的变化为在0∘∼60∘(N/S)0^{\circ}\sim60^{\circ}(N/S)0∘∼60∘(N/S)之间随纬度的变化呈正相关。

其次,针对问题二的日期时角已知而地理位置未知的直杆影长的变化问题,利用最小二乘法建立半搜索模型。通过最小二乘曲线拟合,得到影长随时间变化的多项式函数,利用MATLAB求导得到影长最短的时刻为当地正午12点(北京时间12:35),并根据两地的时差得到当地经度。基于各个时刻影子长度,建立了单目标优化模型,将全球(经纬度)上任一点影长与实际测量值之差的绝对值之和最小为目标函数,约束条件为纬度和杆长的范围,利用半搜索算法,减小搜索范围,得到直杆可能长度,可能的地点(另一种方法是多目标优化,同模型三)。针对问题三,需要由影子顶点坐标确定直杆所处的地理位置和拍摄日期,基于各时刻影长和所测时段内影子方向的变化,基于多目标优化模型建立了循环搜索法,在问题二的基础上,将所测时段内地球上任一点太阳方位角的变化值与实际影子方向角的变化值之差的绝对值之和最小为第二目标函数,约束条件为杆长经度纬度日数范围,利用循环搜索法分别得到可能的杆长、经纬度和日期,残差分析结果表明模型的参数精确度较高。

最后,针对问题四的拍摄视频,先将其转化为AVI格式,利用MATLABPS等软件以3分钟为步长,将视频中直杆影子的顶点坐标转换为实际顶点坐标,通过计算机图像处理相关的透视和投影变化的知识进行分析。并利用日期已知和未知的情况,将问题四的求解分别转化为问题二、三,利用搜索算法得到视频中可能的直杆长度、所处的经纬度及拍摄日期。

关键词: 太阳影子定位;最小二乘;二分查找;遍历算法;多目标优化

二、问题重述

1. 问题背景

太阳影子定位技术是指通过分析实际或视频中物体的影子长度变化确定视频拍摄的地点和拍摄日期的一种方法。

2. 需要解决的问题

问题1:建立影子长度变化的数学模型,分析影子长度关于各个参数的变化规律,并应用建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场(北纬39度54分26秒,东经116度23分29秒)、3米高的直杆的太阳影子长度的变化曲线。

问题2:根据某固定直杆在水平地面上的太阳影子顶点坐标数据,建立数学模型确定直杆所处的地点。并将模型应用于附件1的影子顶点坐标数据,给出若干个可能的地点。

问题3:根据某固定直杆在水平地面上的太阳影子顶点坐标数据,建立数学模型确定直杆所处的地点和日期。将模型分别应用于附件2和附件3的影子顶点坐标数据,给出若干个可能的地点与日期。

问题4:附件4为一根长3米的直杆在太阳下的影子变化的视频,建立数学模型求解出若干个可能的视频拍摄地点。若拍摄日期未知,能否根据视频确定出拍摄地点与日期?

三、问题分析

问题的研究对象是测量的直杆,研究内容为其影子长度的变化情况。该问题描述了在不同测量的地理位置、时间和日期下影子长度的变化,并在知道影子顶点坐标的情况下求解可能的直杆长度及测量地的经纬度和日期。

3.1 问题一的分析

第一小问:问题中说明影子长度与一些参数有关,这就要求我们寻找影子长度涉及哪些相关参数。首先通过相关名词定义,定义了太阳高度角,并通过太阳高度角找到直杆长度和影长的函数关系。其次根据太阳高度角的计算公式,结合直杆长度和影子长度的解析式,得到影子长度随参数变化的函数表达式。结合表达式研究在某个参数值变化,其而它参数值不变的情况下,得到影子长度的变化规律。

第二小问:问题要求我们在已知一些参数的情况下,求影子长度变化曲线。我们在上述中已经求出影子长度函数表达式,题目要求我们研究时间参数在指定时间内变化的情况下,得到影子长度的变化规律。

3.2 问题二的分析

根据直杆在水平地面上影子的顶点坐标求出其所处的若干个可能地点。由问题一知,影子长度随时间变化曲线图具有二次函数的特点,根据附件1运用最小二乘法拟合得到时间与影子长度的曲线方程,从而得到当地影子长度最短的时刻点(一般地定义当地正午12:00影子长度为最短),根据两地的时差以及北京经度(东八区,约为120°E),即可粗略地计算出当地的经度。

经度已知而杆长和纬度未知,可将问题视为一个优化问题,每个时间点的影子坐标都对应一个影长,从而确定关于影长的单目标函数,使每个地点的影子长度与实际测量值之间的累计绝对值之差达到最小,可考虑通过半搜索方法来减小搜索范围并逐步找到最优解(同模型三)。

3.3 问题三的分析

问题三相对于问题二而言,直杆影长的拍摄日期也是未知的。同样将其视为一个多目标优化问题,确定关于影子长度、影子方向的变化角度为另一目标函数,使每个地点的影子长度、影子方向变化角度与实际测量值间的累计绝对差值达到最小,可通过循环搜索算法来实现,第一次搜索出可能的拍摄地点和日期的大范围,第二次搜索减小经纬度、杆长、日期的搜索范围以及步长(提高精度),逐步减小搜索范围找到最优解。

3.4 问题四的分析

视频中的直杆影子变化是在平面上的像,需要通过计算机图像处理中的透视投影成像的相关知识进行分析,将视频图像中直杆影子的顶点坐标转换为实际的顶点坐标。以3分钟为取样步长,取出视频中直杆影子的顶点坐标,并利用MATLABPS等软件转化为实际顶点坐标,问题四就转化为问题二、三的求解,同样地考虑在问题二、三中建立的优化模型,基于分治思想利用计算机遍历算法,通过多次搜索来控制搜索范围从而得到最优解。

四、模型假设

  1. 假设海拔对影子长度的影响可忽略;
  2. 问题四中相机拍摄的角度对影子实际长度的测量影响不大;
  3. 题中所给的直杆顶点坐标数据是实时实地测量得的,是真实可靠的;
  4. 忽略大气折射的影响;
  5. 地球自转为匀角速度;
  6. 太阳光线为平行光.

五、定义与符号说明

5.1 符号说明

符号 符号说明
hhh 太阳高度角
LLL 直杆长度
ddd 影子长度
φ\varphiφ 所测地点的地理纬度
α\alphaα 所测地点的地理经度
δ\deltaδ 赤纬角
ω\omegaω 时角
t′t't′ 所测地点时间
ttt 北京时间
N 日数(以1月1日为第一天,依次类推)

5.2 名词解释

太阳高度角:对于地球上的某个地点,太阳高度角是指某地太阳光线与通过该地与地心相连的地表切线的夹角,即太阳光的入射方向和地平面之间的夹角。

赤纬角:赤纬角又称太阳赤纬,是地球赤道平面与太阳和地球中心的连线之间的夹角。

时角:一个天体的时角定义为天子午圈与天体的赤经圈在北极所成的球面角,或在天赤道上所夹的弧度。

六、模型的建立与求解

6.1 模型一:基于太阳赤纬与时差计算的结构方程模型——问题一的求解

6.1.1 建模准备——影子长度计算公式

根据太阳高度角定义,太阳角为对于地球上的某个地点,太阳高度角是指太阳光的入射方向和地平面之间的夹角hhh。根据该定义画出某一点的太阳角、直杆和影子示意图,如图1所示,即可得到太阳高度角与直杆长度和影子长度关系。
tan⁡h=Ld(6.1)\tan h=\frac{L}{d}\tag{6.1} tanh=dL​(6.1)

图1

太阳高度角计算[1]公式
sin⁡h=sin⁡φsin⁡δ+cos⁡φcos⁡δcos⁡ω(6.2)\sin h=\sin \varphi \sin\delta +\cos \varphi \cos \delta \cos \omega\tag{6.2} sinh=sinφsinδ+cosφcosδcosω(6.2)赤纬角计算[1]公式
sin⁡δ=0.39795⋅cos⁡[0.9853⋅(N−173)](6.3)\sin \delta=0.39795\cdot \cos[0.9853\cdot (N-173)]\tag{6.3} sinδ=0.39795⋅cos[0.9853⋅(N−173)](6.3)即
δ=−23.44∘⋅cos⁡(360365⋅(N+10))\delta=-23.44^{\circ}\cdot \cos\big (\frac{360}{365}\cdot(N+10)\big)\\ δ=−23.44∘⋅cos(365360​⋅(N+10))时角是以正午12点为0°开始算,且时角每一小时增加15°,即
ω=(t′−12)⋅15∘(6.4)\omega=(t'-12)\cdot15^{\circ}\tag{6.4} ω=(t′−12)⋅15∘(6.4)在不同纬度时间不同,计算所测地点时间利用北京时间转化。经度相差1°,时差相差1小时,可得所测地点时间计算公式[2]如下
t′=t−8+α15∘(6.5)t'=t-8+\frac{\alpha}{15^{\circ}}\tag{6.5} t′=t−8+15∘α​(6.5)联立式(6.1)~(6.5),得到影子长度的计算公式
{d=L⋅cot⁡[arcsin⁡(sin⁡φsin⁡δ+cos⁡φcos⁡δcos⁡ω)]δ=−23.44∘⋅cos⁡(360365⋅(N+10))ω=(t′−12)⋅15∘t′=t−8+α15∘(6.6)\begin{aligned} \begin{cases} d=L\cdot\cot\big[\arcsin(\sin\varphi\sin\delta+\cos\varphi\cos\delta\cos\omega)\big]\\ \delta=-23.44^{\circ}\cdot\cos(\frac{360}{365}\cdot(N+10))\\ \omega=(t'-12)\cdot15^{\circ}\\ t'=t-8+\frac{\alpha}{15^{\circ}} \end{cases} \end{aligned}\tag{6.6} ⎩⎪⎪⎪⎨⎪⎪⎪⎧​d=L⋅cot[arcsin(sinφsinδ+cosφcosδcosω)]δ=−23.44∘⋅cos(365360​⋅(N+10))ω=(t′−12)⋅15∘t′=t−8+15∘α​​​(6.6)

6.1.2 影子长度关于各个参数变化规律

通过以上分析,我们已经给出影子长度ddd与纬度φ\varphiφ、日数NNN、时间t′t't′三个参数(由于经度影响时区的确定,进而影响当地时刻的确定,故经度对影子长度的影响可转化为区时对影子的影响)之间的关系式,如式(6.6)所示。式(6.6)中影子长度ddd与各个参数间的解析关系式较为复杂,基于此,我们先确定其中两个参数的值,然后分析影子长度ddd随余下一个参数的变化规律。故我们列出一下三种情况(其中直杆的高度为3米)。

(1) 纬度φ\varphiφ与日数NNN一定,影子长度ddd随时间t′t't′变化规律

我们假设地点为北京天安门广场(北纬39度54分26秒,东经116度23分29秒),日期为2015年10月22日,直杆高度为3米。由于北京时间为东经120度下的区时,故我们先将其转化为东经116度23分29秒位置下的时间。基于(6.4)式,具体的转化公式为:
t′=t−8+116+23/60+29/360015∘(6.7)t'=t-8+\frac{116+23/60+29/3600}{15^{\circ}}\tag{6.7} t′=t−8+15∘116+23/60+29/3600​(6.7)式中,ttt为北京时间。

图2

观察图2,在北京时间9:00-15:00的时间段内影子长度呈先减小后增大的趋势,正午12点时刻达到最小值3.782米,9:00和15:00影子长度分别为7.573 米和6.207米。

clear
clc
%问题一
H=3;%杆长3m
%计算2015年10月22日北京赤纬-11.5637度
sigma=asind(0.39795*cosd(0.98563*(295-173)));
phi=dms2degrees([39,54,26]);%纬度:度分秒转化为度
alpha=dms2degrees([116,23,29]);%经度:度分秒转化为度
t=9:0.25:15;%每间隔15min测量一次
omega=15*t+alpha-300;
L=H*cotd(asind(sind(phi)*sind(sigma)+cosd(phi)*cosd(sigma).*cosd(omega)));
num=length(L); %数据的点数是数列的长度
N0=datenum([0 0 0 9 0 0]); %起始时间9:00
dN=datenum([0 0 0 0 15 0]);
N=N0+(0:num-1)*dN; %根据数据点数,产生15分钟间隔时间点
plot(N,L,'s-');
datetick(gca,'x','HH:MM');
xlabel("时间")
ylabel("影长/米")
(2) 纬度φ\varphiφ与时间t′t't′一定,影子长度ddd随天数N的变化规律

北京天安门广场(北纬39度54分26秒,东经116度23分29秒)为测量地点、正午12时为测量时间、3米为测量直杆长度,天数以15天为间距记录影子长度,作出一年的影子长度随天数变化曲线图。如图下3所示。

图3

通过上图3,可以观察到从第一天起(1月1号为第一天),影子长度随着天数增加而减少,到了第173天(此时为6月22日夏至日),影子长度最短。这天之后,影子长度变为随着天数增加而增加的变化规律。

clear
clc
%问题一第2小问灵敏度分析2,此时改变天数N
H=3;%杆长3m
%计算随天数N改变时北京的赤纬
N=1:15:365;
sigma=asind(0.39795.*cosd(0.98563.*(N-173)));
phi=dms2degrees([39,54,26]);%纬度:度分秒转化为度
alpha=dms2degrees([116,23,29]);%经度:度分秒转化为度
t=12;%以正午12点为例进行灵敏度分析
omega=15*t+alpha-300;
L=H*cotd(asind(sind(phi).*sind(sigma)+cosd(phi).*cosd(sigma)*cosd(omega)));
m=min(L);
d=find(L==min(L));
plot(N,L,'s-');
xlabel("天数")
ylabel("影长/米")
(3) 日数NNN与时间t′t't′一定,影子长度ddd随天数NNN变化规律

以东经116度23分29秒为测量地点经度、2015年10月22日正午12时为测量时间、测量直杆长度为3米,改变纬度参数,作出一年的影子长度随纬度(0°~60°)变化曲线图。

图4

观察图4,我们取纬度的范围,绘制影长纬度的变化曲线。2015年10月22日北京时间正午12:00下影子长度与纬度呈正相关,即直杆的影子长度随着纬度的增大而增大。

clear
clc
%问题一第3小问灵敏度分析3,此时改变维度phi
H=3;%杆长3m
%计算随维度phi改变时北京的赤纬
N=275;
sigma=asind(0.39795*cosd(0.98563.*(N-173)));
phi=[0:3:60];%纬度
%phi=dms2degrees([phi,0,0]);
alpha=dms2degrees([116,23,29]);%经度:度分秒转化为度
t=12;%以正午12点为例进行灵敏度分析
omega=15*t+alpha-300;
L=H*cotd(asind(sind(phi)*sind(sigma)+cosd(phi)*cosd(sigma)*cosd(omega)));
plot(phi,L,'rd-');
xlabel("纬度")
ylabel("影长/米")
(4) 影子长度在指定日期时间下的天安门广场9:00-15:00的变化曲线图

在给定地理位置(北京天安门广场,即北纬39度54分26秒,东经116度23分29秒)、所测日期(2015年10月22日)、区时时间段(北京时间9:00-15:00)、直杆高度(3 米)的情况下,绘制出直杆的太阳影子长度的变化曲线,即是特定参数下太阳影长的变化。这与6.1.2中所给出的情况1(纬度与日数一定,影子长度随区时的变化规律)如出一辙,该条件下直杆的太阳影子长度变化曲线即如图2所示。在北京时间9:00-15:00的时段内影子长度呈先减小后增大的趋势,正午12点达到最小值3.782米,9:00和15:00影子长度分别为7.573米和6.207米。

6.2 模型二:基与最小二乘法模型——问题二求解

6.2.1 建模准备

根据某固定直杆在水平地面上的太阳影子顶点坐标数据,建立数学模型确定直杆所处的地点.将模型应用于表1的影子顶点坐标数据,其中坐标系以直杆底端为原点,水平地面xyxyxy为平面,直杆垂直于地面,测量日期为2015年4月18日),求出若干个可能的地点。

将式(6.1)变形得
arccot(L/d)=h(6.8)arccot(L/d)=h \tag{6.8} arccot(L/d)=h(6.8)两边用sin⁡\sinsin函数作用,得
sin⁡(arccot(L/d))=sin⁡φsin⁡δ+cos⁡φcos⁡δcos⁡ω(6.9)\sin(arccot(L/d) ) =\sin\varphi\sin\delta+\cos\varphi\cos\delta\cos\omega \tag{6.9} sin(arccot(L/d))=sinφsinδ+cosφcosδcosω(6.9)整理得
L=d⋅[1(sin⁡φsin⁡δ+cos⁡φcos⁡δcos⁡ω)2]1/2(6.10)L=d\cdot [\frac{1}{(\sin\varphi\sin\delta+\cos\varphi\cos\delta\cos\omega)^{2}}]^{1/2}\tag{6.10} L=d⋅[(sinφsinδ+cosφcosδcosω)21​]1/2(6.10)式(5.10)中含有 3 个未知数HHH、φ\varphiφ、ω\omegaω,通过方程组求解未知数, 过程繁琐复杂。 因此,引入最小二乘法的曲线拟合方法。给定平面上的点(xi,yi)(x_{i},y_{i})(xi​,yi​),i=1,2,⋯,ni=1,2,\cdots,ni=1,2,⋯,n,用最小二乘法求f(x)f(x)f(x),使δ=∑δi2\delta=\sum\delta_{i}^{2}δ=∑δi2​达到最小。其中δi\delta_{i}δi​为点(xi,yi)(x_{i},y_{i})(xi​,yi​)与曲线y=f(x)y=f(x)y=f(x)的距离。曲线拟合的实际含义是寻求一个函数y=f(x)y=f(x)y=f(x),使在某种准则下与所有数据点最为接近,即曲线拟合得最好。最小二乘准则就是使所有散点到曲线的距离平方和最小,拟合时选用一些简单的“基函数”(例如幂函数,三角函数等等)φ0(x)\varphi_{0}(x)φ0​(x),φ1(x)\varphi_{1}(x)φ1​(x),⋯\cdots⋯,φm(x)\varphi_{m}(x)φm​(x)来线性表示
f(x)=c0φ0+c1φ1+⋯+cmφm(6.11)f(x)=c_{0}\varphi_{0}+c_{1}\varphi_{1}+\cdots+c_{m}\varphi_{m}\tag{6.11} f(x)=c0​φ0​+c1​φ1​+⋯+cm​φm​(6.11)式(6.10)为依据太阳影子长度ddd含杆长LLL 、纬度φ\varphiφ的理论公式, 建立的最小二乘法模型。

图5

6.2.2 数据的预处理

测量日期为 2015 年 4 月 18 日, 根据(6.3)式, 计算出δ=10.0227∘\delta=10.0227^{\circ}δ=10.0227∘。 以直杆底端为原点, 水平地面为xyxyxy平面, 则太阳影子长度d=xi2+yi2d=\sqrt{x_{i}^{2}+y_{i}^{2}}d=xi2​+yi2​​,xix_{i}xi​、yiy_{i}yi​分别是影子顶点xxx轴、yyy轴坐标。利用附件1的数据, 将时间转化为时角,计算出测量时间所对应的影子长度, 如表1所示。

表1

北京时间 影子长度 北京时间 影子长度 北京时间 影子长度
14.7 1.1496 15.05 1.3894 15.4 1.6613
14.75 1.1822 15.1 1.4262 15.45 1.7033
14.8 1.2135 15.15 1.4645 15.5 1.7462
14.85 1.2491 15.2 1.5015 15.55 1.7901
14.9 1.2832 15.25 1.5402 15.6 1.8350
14.95 1.3180 15.3 1.5799 15.65 1.8809
15.0 1.3534 15.35 1.6201 15.7 1.9279

利用MATLAB进行最小二乘法的拟合需要利用未知数(H,φ,ω)(H,\varphi,\omega)(H,φ,ω)的一组估计值作为初值。 对经度进行估计, 根据自然规律: 一天内影子的长度会经历先减小再增大的过程, 而影子长度达到最短的时刻就是测量当地12:00, 即正午太阳高度。 应用 MATLAB 中的cftool工具箱对时角、 及其对应的太阳影子长度进行曲线拟合, 分别假设拟合函数为指数函数、 抛物线、 三阶多项式等并对置信区间、 误差和函数的整体分布特点进行检验, 最终选择二次多项式作为拟合函数, 由表 1 数据得到曲线函数的拟合结果为
y=0.1489t2−3.752t+24.13(6.12)y=0.1489t^{2}-3.752t+24.13\tag{6.12} y=0.1489t2−3.752t+24.13(6.12)式 (6.12) 中 , 拟 合 精 度 指 标 SSE 为1.669×10−51.669\times10^{-5}1.669×10−5, R-square 为1,Adjuster-square 为1, RMSE 为9.63×10−49.63\times10^{-4}9.63×10−4。SSE 及 RMSE 越接近0, R-square和 Adjuster-square 越接近1, 表明拟合效果越好。 式(5.12)为抛物线, 较好地反映出物体影子长度随时角变化的趋势, 并用 MATLAB 绘出式(6.12)函数曲线, 如图 6 所示。

图6

求解直杆所在地经度: 初步判断二次函数对称轴为t≈12.59t\approx12.59t≈12.59, 为了将t=12.59t=12.59t=12.59包含图形中且使图形具有较为明显的规律性, 运用MATLAB软件可得到(6.12)在t∈[9,16]t\in[9,16]t∈[9,16]上的图像, 结果如图 7 所示. 因为北京时间为t=12.59t=12.59t=12.59,即北京时间为 12:35 分左右时, 直杆所在地位于当地正午, 根据两地之间的时差以及北京经度即可计算出直杆所在地经度[3]为111∘E111^{\circ}E111∘E。

图7

clear
clc
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
];%附件1的x坐标
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
];%附件1的y坐标
d=sqrt(x.^2+y.^2);%附件1的影长
t1=14.7:0.05:15.7;
plot(t1,d,'bo');hold on;
t=9:0.25:16;
y=0.1489*t.^2-3.752*t+24.13;
plot(t,y,'k.-.');
xlabel("时间(h)");ylabel("影长(m)");
legend("真实值","拟合曲线");

6.2.3 问题二的求解——基于半搜索算法

1. 模型阶段一:数据点的选取

已知测量日期为 2015 年 4 月 18 日, 即N=108N=108N=108的情形。 根据(5.3)式, 求得此时太阳的赤纬角δ=10.0227∘\delta=10.0227^{\circ}δ=10.0227∘, 对应时角从14.7∼15.714.7\sim15.714.7∼15.7, 时间间隔为 0.05,取当地正午太阳高度, 即北京时间 12:35, 此时直杆影长为最短 0.4934 m。 同时,选取正午的另一个好处是cos⁡ω=cos⁡0=1\cos\omega=\cos0=1cosω=cos0=1, 便于后续计算。

2. 模型阶段二:半搜索算法思路

设置杆长 L的搜索范围为1∼6m1\sim6m1∼6m, 步长为0.05; 纬度的搜索范围为−90∘S∼90∘N-90^{\circ}S\sim90^{\circ}N−90∘S∼90∘N, 步长0.1。 基于前述的处理, 问题二的半搜索算法[4]的时间复杂度为O(n2)O(n^{2})O(n2), 大大缩短了搜索的时间并提高了精度。 最后, 根据半搜索的结果得到(三组)可能的地点列表如下。

表2

组别 杆长(m) 经度(度) 纬度(度)
第一组 1.8 109E -19.3S
第二组 1.9 109E -1.3S
第三章 2 109E 16.7N

经检验, 第一组参数与题给数据严重不符, 所以舍去, 对于二、 三组参数,根据作图及统计分析, 发现能够很好地描述所给日期时间段内直杆影长的变化。

clear
clc
%第二问方法一
L=1:0.01:3;%设置杆长的搜索范围
phi=-20:0.01:30;%设置搜索的纬度范围
N=31+28+31+18;%计算该日期下的日数
delta=asind(0.39795*cosd(0.98563*(N-173)));%该日下的赤纬角
A=sind(delta);B=cosd(delta);
d=0.4934;K={};k=0;%当地正午12点对应的杆影长度
for i=1:length(L)for  j=1:length(phi)D=L(i)*cotd(asind(A*sind(phi(j))+B*cosd(phi(j))));if abs(d-D)<0.000008k=k+1;list=[L(i),phi(j)];K{k}=list;list=[];elsecontinue;endend
end
3. 方法二:循环遍历搜索法

这是根据循环遍历搜索法做出的拟合度较高的曲线(思路见模型三)。

图8

此法同样适用于第三、四问的求解(可改变步长得到更高的精度)。

clear
clc
%第2问方法二遍历搜索
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坐标
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坐标
d=sqrt(x.^2+y.^2);%杆影长
K={};k=0;%KK=[];KKK=[];
alpha=0:1:180;%经度
phi=-90:1:90;%纬度
L=1:0.1:3;%杆长
N=31+28+31+18;%计算该日期下的日数
delta=asind(0.39795*cosd(0.98563*(N-173)));%该日下的赤纬角
t=14.7:0.05:15.7;%拍摄时的北京时间
A=sind(delta);B=cosd(delta);
for i=1:length(alpha)for j=1:length(phi)for m=1:length(L)flag1=0;flag2=0;%标记%计算目标函数1影长的误差omega=15*t+alpha(i)-300;%根据北京时间计算的时角%计算各地(经纬)不同北京时间的影长D=L(m)*cotd(asind(A*sind(phi(j))+B*cosd(phi(j)).*cosd(omega)));%KK1=sum(abs(d-D));KK=[KK,KK1];if (sum(abs(d-D)))<0.5flag1=1;%目标函数1满足要求elsecontinueend%计算目标函数2方向角的误差theta1=atand(y(1)/x(1));%实测第1个北京时刻与x轴的夹角theta2=atand(y(end)/x(end));%实测从第2个北京时刻开始与x轴的夹角theta0=abs(theta2-theta1);%实测20个时段内与第1个时刻的方位角的%首先求非实测的thetas1h1=asind(A*sind(phi(j))+B*cosd(phi(j))*cosd(omega(1)));thetas1=acosd((sind(h1)*sind(phi(j))-sind(delta))/(cosd(h1)*cosd(phi(j))));%再求从第2项开始的thetaskh=asind(A*sind(phi(j))+B*cosd(phi(j))*cosd(omega(end)));thetask=acosd((sind(h)*sind(phi(j))-sind(delta))/(cosd(h)*cosd(phi(j))));thetaalpha=abs(thetask-thetas1);%非实测20个时段内与第1个时刻的方位角的%KK2=sum(abs(thetaalpha-theta0));KKK=[KKK,KK2];if (sum(abs(thetaalpha-theta0)))<0.1146flag2=1;elsecontinueendif flag1==1 && flag2==1k=k+1;list=[alpha(i) phi(j) L(m)];K{k}=list;list=[];endendend
end

6.2.4 模型检验与结果分析

图9

我们对计算出的地点(经纬度)和直杆长度检验,代入计算的参数值计算14:42 到 15:42 时段内的影长与实际影长做比较并计算残差值。计算出题中所给附件 1 中各个时刻(时角)的影子长度。对于地理位置(东经109度,北纬16.7度),在北京时间2015年4月18日14:42到15:42时段内, 利用式(6.9)计算出影长。利用MATLAB软件作出真实值计算值图像,如图8、9所示,通过图8、9的观察,可以看出真实值与计算值相差不大,两条折线基本吻合,模型计算出结果较可靠。作出了两组数据的残差值,并画出含有 95%CI 的残差图,发现所有直线均包含了0,没有异常点的出现,所以说明我们的模型和参数都是正确的,见图10。

图10

综合上面的检验和分析, 我们可以得出结论: 所构建模型的误差值极小, 故建立的模型及求解结果具有较高的可靠性与精确度(下面的程序可用于检验其他标准)。

clear
clc
%第二问残差检验
L=1.2;%杆长
alpha=119.5;%经度
phi=-10;%纬度
%杆长1.8m时,纬度19N,经度112.5E;杆长1.2m时,纬度-10S,经度119.5E
N=31+28+31+18;%计算该日期下的日数
delta=asind(0.39795*cosd(0.98563*(N-173)));%该日下的赤纬角
A=sind(delta);B=cosd(delta);
t=14.7:0.05:15.7;%设置当天时间范围
omega=15*t+alpha-300;%转化为时角
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坐标
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坐标
d=sqrt(x.^2+y.^2);%当地时角的杆影长
D=L*cotd(asind(A*sind(phi)+B*cosd(phi).*cosd(omega)));%计算的杆影长
figure(1);
plot(D,'s-');hold on;plot(d,'o-');d=d';D=D';
figure(2);
[b,bint,r,rint,stats]=regress(D,d);
rcoplot(r,rint);

6.3 模型三: 基与分治思想的循环遍历法——问题三的求解

6.3.1 建模准备

附件2、附件3中的影子长度的测量时间段分别是北京时间12:41-13:41、北京时间13:09-14.09,每隔3分钟测一次,共21次(i=1,2,⋯,21)(i=1,2,\cdots,21)(i=1,2,⋯,21)。

表3

符号 符号说明
d′(i)d'(i)d′(i) 附件2中每个时刻实测影子长度
d′′(i)d''(i)d′′(i) 附件3中每个时刻实测影子长度
d0(i)d_{0}(i)d0​(i) 第iii个时刻实测影子长度(i=1,2,⋯,21)(i=1,2,\cdots,21)(i=1,2,⋯,21)
dαφ(i)d_{\alpha\varphi}(i)dαφ​(i) 全球范围内任一点(α\alphaα表示该未知的经度坐标,φ\varphiφ表示该位置的纬度坐标)第iii个时刻的影长

另一方面,在实测北京时间段内,直杆影子方向的变化角度可通过几何知识得到。建立直角坐标系,以直杆底端为端点,水平地面为xyxyxy平面,坐标数值单位为米。其中符号说明如下

表4

符号 符号说明
θαφ(i)\theta_{\alpha\varphi}(i)θαφ​(i) 全球范围内任一点(α\alphaα表示经度,φ\varphiφ表示纬度)第个时刻影子方向与第1个时刻影子方向之间的夹角(i=1,2,⋯,20)(i=1,2,\cdots,20)(i=1,2,⋯,20)
θi\theta_{i}θi​ 附件2中第iii个时刻影子方向与xxx轴的夹角(第iii个时刻12:41)
θi′′\theta_{i}''θi′′​ 附件3中第iii个时刻影子方向与xxx轴的夹角(第iii个时刻13:09)
θ0(i)\theta_{0}(i)θ0​(i) 第iii个时刻影子方向与第1个是时刻影子方向的夹角

根据上述关系式可得如下模型,其中
θ0=∣θ1−θi∣(6.13)\theta_{0}=|\theta_{1}-\theta_{i}|\tag{6.13} θ0​=∣θ1​−θi​∣(6.13)根据几何知识,有
tan⁡θ1=y1x1,tan⁡θi=yixi(6.14)\tan\theta_{1}=\frac{y_{1}}{x_{1}},\tan\theta_{i}=\frac{y_{i}}{x_{i}}\tag{6.14} tanθ1​=x1​y1​​,tanθi​=xi​yi​​(6.14)

6.3.2 多目标优化模型的建立

根据建模前的准备,建立以下多目标优化模型:

实际影长与计算影长绝对值之差之和最小,有
min⁡∑i=121∣dαφ(i)−d0(i)∣(6.15)\min \sum_{i=1}^{21}|d_{\alpha\varphi}(i)-d_{0}(i)|\tag{6.15} mini=1∑21​∣dαφ​(i)−d0​(i)∣(6.15)实际方向角的变化与计算方向角的变化之差之和最小,有
min⁡∑i=120∣θαφ(i)−θ0(i)∣(6.16)\min \sum_{i=1}^{20}|\theta_{\alpha\varphi}(i)-\theta_{0}(i)|\tag{6.16} mini=1∑20​∣θαφ​(i)−θ0​(i)∣(6.16)循环搜索的范围,得如下约束条件
s.t.{−180∘≤α≤180∘−90∘≤φ≤90∘0<L<60<N<365,N为整数(6.17)\begin{aligned}s.t. \begin{cases} -180^{\circ}\leq\alpha\leq180^{\circ}\\ -90^{\circ}\leq\varphi\leq90^{\circ}\\ 0<L<6\\ 0<N<365,N为整数 \end{cases} \end{aligned}\tag{6.17} s.t.⎩⎪⎪⎪⎨⎪⎪⎪⎧​−180∘≤α≤180∘−90∘≤φ≤90∘0<L<60<N<365,N为整数​​(6.17)其中,α\alphaα表示经度,φ\varphiφ表示纬度,LLL表示直杆的实际长度,NNN表示日期。联立上述(6.15)~(6.17)式,即得循环搜索算法的目标函数与约束条件。

在优化目标中,dαφ(i)d_{\alpha\varphi}(i)dαφ​(i)的计算公式为:
{dαφ(i)=L⋅cot⁡[arcsin⁡(sin⁡φsin⁡δ+cos⁡φcos⁡δcos⁡ω)]δ=arcsin⁡{0.39795⋅cos⁡[0.98563⋅(N−173)]}ti=(t0i−12)⋅15∘(6.18)\begin{aligned} \begin{cases} d_{\alpha\varphi}(i)=L\cdot \cot[\arcsin(\sin\varphi\sin\delta+\cos\varphi\cos\delta\cos\omega)]\\ \delta=\arcsin\left\{ 0.39795\cdot\cos[0.98563\cdot(N-173)] \right\}\\ t_{i}=(t_{0i}-12)\cdot15^{\circ} \end{cases} \end{aligned}\tag{6.18} ⎩⎪⎨⎪⎧​dαφ​(i)=L⋅cot[arcsin(sinφsinδ+cosφcosδcosω)]δ=arcsin{0.39795⋅cos[0.98563⋅(N−173)]}ti​=(t0i​−12)⋅15∘​​(6.18)其中݀d0(i)=xj2+yj2d_{0}(i)=\sqrt{x_{j}^{2}+y_{j}^{2}}d0​(i)=xj2​+yj2​​,将太阳方位角记为θs\theta_{s}θs​,计算公式为:
cos⁡θs=sin⁡(h)sin⁡φ−sin⁡δcos⁡(h)⋅cos⁡φ(6.19)\cos\theta_{s}=\frac{\sin(h)\sin\varphi-\sin\delta}{\cos(h)\cdot\cos\varphi}\tag{6.19} cosθs​=cos(h)⋅cosφsin(h)sinφ−sinδ​(6.19)θαφ(i)\theta_{\alpha\varphi}(i)θαφ​(i)的计算公式为:
θαφ(i)=∣θs1−θsi∣(6.20)\theta_{\alpha\varphi}(i)=|\theta_{s1}-\theta_{si}|\tag{6.20} θαφ​(i)=∣θs1​−θsi​∣(6.20)

6.3.3 模型的求解

1. 模型阶段一:数据预处理

首先,我们对附件 2、附件3中的数据处理,利用它们的顶点坐标计算出它们实测影长。我们以݀d0(i)d_{0}(i)d0​(i)作为第i个时刻实测影长。对于太阳方位角,我们通过附件2、附件3中的顶点坐标根据几何知识求出方位角。再将第i个时刻的方位角依次与第1个时刻的方位角作差的绝对值。我们以θ0(i)\theta_{0}(i)θ0​(i)作为附件 2、附件3第i 个时刻方位角与第1个时刻方位角差的绝对值。

其次,我们基于分治思想的循环遍历法思想,对经纬度、杆长、天数参数值进行遍历,将这些参数值代入式(6.19)和式(6.23)中,计算出时间段内的影长和方位角,再依次将第i个时刻后的方位角与第1个时刻方位角作差取绝对值。

2. 模型阶段二:基于分治思想的循环遍历法

首先,设置杆长 L的遍历范围为1∼6m1\sim6m1∼6m,步长为0.1; 纬度φ\varphiφ的遍历范围为−90∘S∼90∘N-90^{\circ}S\sim90^{\circ}N−90∘S∼90∘N,步长为1; 经度α\alphaα的遍历范围为−180∘W∼180∘E-180^{\circ}W\sim180^{\circ}E−180∘W∼180∘E,步长为1;天数NNN的遍历范围为1到365。 我们将以上计算出来的θ0(i)\theta_{0}(i)θ0​(i) 、θ0′(i)\theta _{0} '(i)θ0′​(i)与θ0′′(i)\theta_{0}''(i)θ0′′​(i)作差后绝对值进行累加,若不满足目标函数min⁡∑i=121∣dαφ(i)−d0(i)∣<0.004\min \sum_{i=1}^{21}|d_{\alpha\varphi}(i)-d_{0}(i)|<0.004min∑i=121​∣dαφ​(i)−d0​(i)∣<0.004,则进行下一次遍历循环;反之,继续筛选。若满足目标函数min⁡∑i=020∣θαφ(i)−θ0(i)∣<0.004\min \sum_{i=0}^{20}|\theta_{\alpha\varphi}(i)-\theta_{0}(i)|<0.004min∑i=020​∣θαφ​(i)−θ0​(i)∣<0.004,将该次遍历的参数值记录下来。继续遍历,直到遍历结束(可以减小步长)。

3. 模型阶段三:求解结果

我们遍历得到附件2、附件3的两组可能结果。

表5

序号 日期(天) 杆长(m) 经度(度) 纬度(度)
1 197 1.90 79 39
2 145 1.90 79 39

分析表3, 根据附件2中所给出的直杆影子顶点坐标,最终找出了直杆长度和所处地点(经纬度)。结果为:杆长1.90米;日期为7月16日对应地点为(东经79度,北纬39度)附近。 日期为5月25日对应地点为(东经79度,北纬39度)附近。

表6

序号 日期(天) 杆长(m) 经度(度) 纬度(度)
1 225 3 127 69
2 145 3 109 -27

据表 4,根据附件2中所给出的直杆影子顶点坐标,最终找出了直杆长度和所处地点(经纬度)。结果为:杆长3米;日期为 8月13日对应地点为(东经127度,北纬69度)附近。日期为5月25日对应地点为(东经109度,南纬27度)附近。

图11

电脑性能好的经度可以设为-180~180

clear
clc
%第3问附件2第1次搜索
x=[-1.2352,-1.2081,-1.1813,-1.1546,-1.1281,-1.1018,...-1.0756,-1.0496,-1.0237,-0.9980,-0.9724,-0.9470,...-0.9217,-0.8965,-0.8714,-0.8464,-0.8215,-0.7967,...-0.7719,-0.74730,-0.7227];%杆影的x坐标
y=[0.1730,0.1890,0.2048,0.2203,0.2356,0.2505,...0.2653,0.2798,0.2940,0.3080,0.3218,0.3354,...0.3488,0.3619,0.3748,0.3876,0.4001,0.4124,...0.4246,0.4366,0.4484];%杆影的y坐标
d=sqrt(x.^2+y.^2);%杆影长
K={};k=0;%KK=[];KKK=[];
alpha=0:1:90;%经度
phi=-50:1:50;%纬度
L=1:0.1:3;%杆长
N=1:1:250;%日数
t=12.6833:0.05:13.6833;%拍摄时的北京时间
for i=1:length(alpha)for j=1:length(phi)for m=1:length(L)for n=1:length(N)flag1=0;flag2=0;%标记%计算目标函数1影长的误差delta=asind(0.39795*cosd(0.98563*(N(n)-173)));%赤纬角A=sind(delta);B=cosd(delta);omega=15*t+alpha(i)-300;%根据北京时间计算的时角%计算各地(经纬)不同北京时间的影长D=L(m)*cotd(asind(A*sind(phi(j))+B*cosd(phi(j)).*cosd(omega)));%KK1=sum(abs(d-D));KK=[KK,KK1];if (sum(abs(d-D)))<0.1flag1=1;%目标函数1满足要求elsecontinueend%计算目标函数2方向角的误差theta1=atand(y(1)/x(1));%实测第1个北京时刻与x轴的夹角theta2=atand(y(2:end)./x(2:end));%实测从第2个北京时刻开始与x轴的夹角theta0=abs(theta2-theta1);%实测20个时段内与第1个时刻的方位角的变化%首先求非实测的thetas1h1=asind(A*sind(phi(j))+B*cosd(phi(j))*cosd(omega(1)));thetas1=acosd((sind(h1)*sind(phi(j))-sind(delta))/(cosd(h1)*cosd(phi(j))));%再求从第2项开始的thetaskh=asind(A*sind(phi(j))+B*cosd(phi(j)).*cosd(omega(2:end)));thetask=acosd((sind(h).*sind(phi(j))-sind(delta))./(cosd(h).*cosd(phi(j))));thetaalpha=abs(thetask-thetas1);%非实测20个时段内与第1个时刻的方位角的变化%KK2=sum(abs(thetaalpha-theta0));KKK=[KKK,KK2];if (sum(abs(thetaalpha-theta0)))<0.6flag2=1;elsecontinueendif flag1==1 && flag2==1k=k+1;list=[alpha(i) phi(j) L(m) N(n)];K{k}=list;list=[];endendendend
end
%修改步长进行第二次小范围搜索

误差如下图

图12

clear
clc
%第三问附件2残差检验
L=2;%杆长
alpha=80;%经度
phi=38;%纬度
N=209;%日数
delta=asind(0.39795*cosd(0.98563*(N-173)));%该日下的赤纬角
A=sind(delta);B=cosd(delta);
t=12.6833:0.05:13.6833;%设置当天时间范围
omega=15*t+alpha-300;%转化为时角
x=[-1.2352,-1.2081,-1.1813,-1.1546,-1.1281,-1.1018,...-1.0756,-1.0496,-1.0237,-0.9980,-0.9724,-0.9470,...-0.9217,-0.8965,-0.8714,-0.8464,-0.8215,-0.7967,...-0.7719,-0.74730,-0.7227];%杆影的x坐标
y=[0.1730,0.1890,0.2048,0.2203,0.2356,0.2505,...0.2653,0.2798,0.2940,0.3080,0.3218,0.3354,...0.3488,0.3619,0.3748,0.3876,0.4001,0.4124,...0.4246,0.4366,0.4484];%杆影的y坐标
d=sqrt(x.^2+y.^2);%当地时角的杆影长
D=L*cotd(asind(A*sind(phi)+B*cosd(phi).*cosd(omega)));%计算的杆影长
figure(1);
plot(t,D,'s-');hold on;plot(t,d,'o-');d=d';D=D';set(gca,'ygrid','on');
legend("预测值","真实值");xlabel("time/h");ylabel("影长/m");
figure(2);
[b,bint,r,rint,stats]=regress(D,d);
rcoplot(r,rint);

6.4 模型四:基与坐标变换的图形射影模型——问题四的求解

6.4.1 建模准备

首先将拍摄到的视频转换成 AVI 格式,用MATLAB以3分钟为间隔进行信息提取,共提取到14张图片,每张图片记为f(i)f(i)f(i),(i=1,2,⋯,14)(i=1,2,\cdots,14)(i=1,2,⋯,14),然后将f(i)f(i)f(i)转换成灰度图݂f′(i)f'(i)f′(i),如图7所示(此时以提取的第8张图片为例)。并通过PhotoshopCS5等软件读取到直杆底端和影子顶点的坐标数据,得到影长。最后根据太阳影子定位模型进行求解。

图13

通过实验得出影长的像素相对值,通过 PS 软件的处理,分别14个影子对应比例的实际影子长度。

6.4.2 第 1 小问——拍摄日期已知

第1小问的求解方法同第二问的曲线拟合的最小二乘法,对于表5中所给出的13个时刻下的影长实际值,我们用曲线拟合的方法绘制出影长随时间的变化曲线。根据5.1.2中情况(1)的结论:一天中影长随时间ttt呈先减小后增大的趋势,且在当地正午时刻达到最小值,故影长随时间的变化曲线为二次函数的形式。我们考虑用 MATLAB 软件进行拟合,最后选定当地正午太阳高度的北京时间作为参考点进行半遍历搜索,此法遍历次数少精度高。

在拟合图像中,曲线最低点对应的横坐标时间值为12.6953h12.6953h12.6953h,即为北京时间12:41,同时也对应于当地时间正午 12 点。基于式(5.3)的转化关系,便可计算出视频拍摄地点的经度值为东经109.65 度。式(5.4)中, t为北京时间12 点 41 分,t′t't′为当地时间的正午12点。

6.4.3 第 2 小问——拍摄日期未知

第2小问的求解方法同第三问的循环遍历搜索算法,首先第一次搜索运用小精度在大范围内查找满足的参数范围,第二次搜索提高精度减小步长在小范围内搜索,已知杆长为2m的情况下,给出多目标优化的搜索模型。

表5略

在表 5 中,给出了视频拍摄的可能地点和日期。日期为1月18日时,位置坐标为(东经 111.65 度,南纬41度)附近;日期为1月18日时,位置坐标为(东经114.65度, 北纬40度);日期为3月9日时,位置坐标为(东经114.65度,北纬40度)附近;日期为6月21日时,位置坐标为(东经109.65度,北纬42度)附近;日期为7月14日时,位置坐标为(东经 110.65度,北纬41度) 附近;日期为8月5日时,位置坐标为(东经114.65度,北纬40 度) 附近;日期为11月24日时,位置坐标为(东经111.65度,南纬41度)附近;日期为12月12日时,位置坐标为(东经109.65度,南纬42度)附近.

七、模型的评价与推广

7.1 模型的评价

7.1.1 模型的优点

7.1.2 模型的缺点与不足

7.2 模型的改进

参考文献

[1] 王昌明.可照时数和太阳高度角计算公式的简化证明[J].山东气象,1989(02):46-48.
[2] 国安,米鸿涛,邓天宏,李亚男,李兰霞.太阳高度角和日出日落时刻太阳方位角一年变化范围的计算[J].气象与环境科学,2007(S1):161-164.
[3] 仝青山,温石刚,王轩,李明玉.基于数学建模的太阳影子定位[J].南京工程学院学报(自然科学版),2017,15(01):33-36.
[4] 姜启源.数学模型[M].北京:高等教育出版社,2011.1.

附录

2015年数模A题太阳影子定位学习笔记相关推荐

  1. 全国大学生数学建模竞赛-2015-A题-太阳影子定位

    2015高教社杯全国大学生数学建模竞赛题目 A题 太阳影子定位 如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面,太阳影子定位技术就是通过分析视频中物体的太阳影子变化,确定视频拍摄的地点和日期 ...

  2. 2021年研究生数模B题论文记录

    2021年研究生数模B题论文记录 1.常见数据处理方法: 2.相关性系数选择 3.聚类算法 4.一种数据降维方式 5.预测模型 文章来源 2021年全国大学生研究生数学建模竞赛优秀论文集合,B题,文章 ...

  3. 08年A题数码相机定位学习笔记

    To Be Continue- 文章目录 数码相机定位学习笔记 一.摘要 二.问题分析 三.刚体变换 3.1世界坐标系→\rightarrow→相机坐标系 3.2相机坐标系→\rightarrow→图 ...

  4. 2022国赛数模A题思路以及解析(附源码 可供学习训练使用)

    需要全部源码或者论文请点赞关注收藏后评论区留言 前言 发布已获得创作队伍的同意,论文最终斩获省一等奖,写的十分优秀,可供后面的数模比赛训练参考 摘要 基于摇荡模型的波浪能装置最大输出功率设计问题研究 ...

  5. 2022高教杯数学建模E思路 超详细文字内容 数模E题

    E 题 小批量物料的生产安排 某电子产品制造企业面临以下问题:在多品种小批量的物料生产中,事先无法知道物料的实际需求量.企业希望运用数学方法,分析已有的历史数据,建立数学模型,帮助企业合理地安排物料生 ...

  6. 2023年MathorCup数模D题赛题解题思路

    MathorCup俗称妈杯,是除了美赛国赛外参赛人数首屈一指的比赛,而我们的妈杯今天也如期开赛.今年的妈杯难度,至少在我看来应该是2023年截至目前来讲最难的一场比赛.问题的设置.背景的选取等各个方面 ...

  7. 2023年MathorCup数模A题赛题详细思路

    MathorCup俗称妈杯,是除了美赛国赛外参赛人数首屈一指的比赛,而我们的妈杯今天也如期开赛.今年的妈杯难度,至少在我看来应该是2023年截至目前来讲最难的一场比赛.问题的设置.背景的选取等各个方面 ...

  8. 2023年MathorCup数模A题赛题

    A 题 量子计算机在信用评分卡组合优化中的应用 在银行信用卡或相关的贷款等业务中,对客户授信之前,需要先通过 各种审核规则对客户的信用等级进行评定,通过评定后的客户才能获得信 用或贷款资格.规则审核过 ...

  9. 2023年MathorCup数模B题赛题

    B 题 城市轨道交通列车时刻表优化问题 列车时刻表优化问题是轨道交通领域行车组织方式的经典问题之一. 列车时刻表规定了列车在每个车站的到达和出发(或通过)时刻,其在实 际运用过程中,通常用列车运行图来 ...

  10. 2023五一数模b题思路分享

    B题:快递需求分析问题 网络购物作为一种重要的消费方式,带动着快递服务需求飞速增长,为我国经济发展做出了重要贡献.准确地预测快递运输需求数量对于快递公司布局仓库站点.节约存储成本.规划运输线路等具有重 ...

最新文章

  1. sql server面试题
  2. 计算机网络解决数据包丢失,数据包丢失时网络控制系统的稳定性分析及设计
  3. 我害怕接入IM云的开发者
  4. Javascript Array和String的互转换。
  5. [kubernetes] 资源管理 --- 资源预留实践
  6. 11尺寸长宽 iphone_新手必知LED显示屏尺寸规格及计算方法
  7. 什么是地址译码 理解二进制编码
  8. SAP中的默认帐户与密码
  9. 硬盘分区后的逻辑结构
  10. j2me联网_与J2ME联网
  11. 完美简单的集成高德地图导航和语音播报功能
  12. Python——实现防止微信撤回消息
  13. 亚马逊aws认证是什么?亚马逊aws认证证书含金量怎么样?
  14. 第14章 Python网络爬虫
  15. 模电—初探MOSFET
  16. 改变exe文件图标的方法
  17. 支付宝花呗额度一直不涨?阿里老员工说出原因,亲测有效
  18. Pycharm导入同级目录模块解决办法汇总
  19. 办理广东林业调查规划设计资质最新申报标准
  20. CityMaker学习教程14 水面图层的创建

热门文章

  1. 埃及金字塔之谜最完美的解释
  2. (深度原创)华为基于LTC主流程的组织销售能力提升,含相关工具模板方法!
  3. 苹果iPad手机如何无线投屏电脑使用教程
  4. Fgui: Glist 实现无限滑动 虚拟列表
  5. hbuild html5打包apk,使用HBuilder打包5+App
  6. 于仕琪libfacedetection WIN10 VS2015
  7. Awvs 12.x安装教程及常见问题
  8. 关于XP系统远程桌面的一点点记录
  9. 数学建模论文格式要求汇总
  10. 智慧教育教学案例分析