集-国赛获奖1992-年2015优秀

/定位参数分析影长的变化规律，通过中间变量高度角、赤纬角、时角确立影长变化模型。其次利用影长变化模型求解出在900-1500影长变化曲线类似一12:123.84m。最后分析响程度，得到影长变化的相对误差限为0.23%。针对问题二，1到两个稳定的直杆位置为109012E,19012N和95024E25012N，利用从文献之差与对应方位角之差相减，建立误差分析模型，对结果进行分析，得到109012E,19012'N准确性更高进行求解，得到附件2，附件3的直杆所处位置与日期见下表。位置坐标7105.4910-11251.092510-使用附件1的数据检验问题三模型，求出日期结果相对误差为0.05%，证明模针对问题四：首先对中提取的影长进行灰度化，二值化，小连通域处理，运用Canny边缘检测算法提取预处理后的轮廓，通过变换提取出杆影顶端坐标。建立单目标优化模型，进行全局搜索得到拍摄地点为' ：拟迭代法人工鱼群算法多目标优化模型Canny边缘检测算问题重述问题的背景确定拍摄的地点和日期有待研究。问题的相关信息根据附件提供的相关信息，以坐标系以直杆底端为原点，水平地面为xy平需解决的问题问题一：以9:0015:00为时间，天安门广场（北纬39542，东经）为地点，建立描述长度变化的数学模型，分设若干参数分析高3米的直杆长度的变化规律；问题二：为了确定在水平地面上固定直杆的拍摄地点，以直杆的 问题三：为了确定在水平地面上固定直杆的拍摄地点和日期，以直杆问题四：（1）根据附件4所给，建立确定 解若干个可能的拍摄地点；（2）将附件4所给 模型假设与符号说明模型的假设假设2： 假设3： 模型的符号说明赤纬角，的直射纬Ad，，NE时差，真时与平时的差CMx,x,问题分析问题一的分析影长与上述五个参数之间的关系。查阅相关资料得到可以高度角、赤玮广场的5个参数带入影长变化模型，可得到杆影的变化曲线，分析长问题二的分析1的顶点坐标数据求解直杆位置。可视为已知影长坐标、日期和时刻，求影长所在的地点的问题。首先应根据影长坐标计算影长，本文拟将附件1中的影长、时刻及日期代入问题一中的影长变化函数中，可得到含有21个方程式的值最小为单目标的最优化模型对附件1进行求解其中约束条件可考虑的范问题三的分析日期的问题。首先应将附件2、3中的21组影长及时刻代入问题一中的影长变化影长差值最小、影长的曲线斜率差最小为目标的双目标最优化模型对附件2、3问题四的分析针对问题四首先需对中长度进行提取，通过对图像预处理，将彩色图像灰度化，再对小连通域处理。拟采用一种Canny边缘检测算法提取预处理后的数据处理实际影长附件实际影长（详细见附录1）如下表：………………角度变化角度见附录2）如下图所示：32ZZ10

角度变化

1X

21节气日期划分我国一年365天，共24节气[1]，由于地球在黄经轨道上围绕公转，造至、秋分、冬至将一年365天划分为四季如下表：2问题一的解答模型的准备背景知识高度角hs：相对于地面的高度角sinhssinsincoscos赤纬角：的直射纬度，根据Cooper近似方程可得赤纬角23.45sin360284n 方位角A：从北方沿着地平线顺时针度角的角sinAcoscos时角（1）时、时真太阳时T1：连续两次经过上中天的时间间隔平时：日常生活中所使用的时间ET1E9.87sin2B7.53cosB式中B360n81n也要考虑地球轨道的偏心率和地球的轴向倾斜影响，所以也要加入时差E。TC4LgT3 式中T150 所以当真时可以通过使用前面两个修正因子调整本地时间 中得到TT4LgT3 影长出现时间段由于直杆的影长在出现的时候才有，所以应确定每天出现的时间天内出现的时间。地球赤纬角的表达式[2]为23.44sin360284 日出、日没时处于地平线位置，高度角hs0，上式变为0=sinsincoscoscossinsincoscos-tan=arccos-tantan从上式可见，日出时1=-，日没时2=，日出、日没时角和当地纬度和太阳赤纬角有关，因为-1cos1，所以当=0时，北半球在正午出现在正现极昼，即大于24h的白昼。当高度角hs0时，此时为日出、日没的方位角cosAcos N=2arccostantan式中N影长变化模型的建立(x,y)f,,n,t,时间四个变量有关，由此可建立以高度角、方位角、时角、赤纬角、杆的高度为变量与端点位置坐标的数学模型，得到影长关于五个变量以直杆为坐标原点O，正南方向为X轴，正东方向为Y轴，直杆为Z轴建立平面直角坐标系。假设的端点为Qxy，故长度为落影点Q点Otanh

x2 /而方位角方程满足：yxtan联立上面两式可得端点坐标yxtanhcotx 11tan hcot hcot htanAcoths1tan2 1tan2Lhsarcsinsinsincoscoscos cossinA=arcsin cos

284n23.45sin360 模型的求解与分析影长变化曲线的求解根据日出、日没的太阳时角确定天安门广场每日直杆影长在632am1728pm区间内用软件（源程序见附录3）求解，得到2015年10月22日时间9:00-15:00之间广场3米高的直杆影长及高度角如下表：表3广场9:00-15:00所求数时间高度角（度）杆影长广场9点到15 影长变化曲广广场9点到15 影长变化曲87的的 654

时

0影长变影长变化 时间图 图 5.3.2影长变化曲线的分析由图1可得，基于五个变量的影长变化规律符合 最短影长l3.84(m)位于时间12点12分。由于 时间指东经120°，即的位置，故最短影长应出现在时间12点之后。，由图2可得广场影长在9:00-15:00时间段内变化的趋势呈先减，各个参数与影长的变化规律根据上述建立的影长变化模型，长度与经度、纬度、日期、时间以及杆日期参数11月31之间变化。利用软件编程（源程序见附录4）求解出日期变化对太阳影长与日期的关系6度4度的长2的长率率影长的变化0影长的变化

日期序号日期对影长的影响 日期的变化图4日期变化对影长的关表4四节气对长度变化的影节日变化率节日变化率320321--922923621622-9.64E-12211222--经纬度、时刻及杆长参数变某一个变量得到这一变量与影长的变化规律，利用软件（源程序见附录4）分别求解出影长与经度、纬度、时刻及杆长的变D化规律如下图：A影长与杆长度A

太阳影长与度的BB的长的长度0 直杆的长C太阳影长与度的C的长的长度10 纬度：30°S-

8的长的长度46080100120140160经度:60°E-D太阳影长与间的D的长的长度64 天内的图5影长与杆长、及时刻的变化规使杆长在0m10m之间变化，随着直杆长度的增加，影长与杆长呈线性相关。由图B可得，控制直杆所在纬度、日期、时刻及杆长不变，使直杆的经度在600W1800W上变化，当经度达到120°时，影长达到最小值。由图C可得，控制直杆所在经度、日期、时刻及杆长不变，使直杆的纬度在300S300N上变化，确定广场每日直杆影长6:32am17:28pm区间内变化，通过观察，发现12点之前影长逐步减少，12误差分析由N

d

赤玮角与天数有关，平均误差为1.71%高度角与纬度、天数、经度有关，平均误差为0.57%。最后求解得到影长变化误差限为0.23%。问题二的解答虑的范围；最后通过穷举搜索，得到直杆所在的地点。单目标优化模型的建立确定目标函数minll 确定约束条件

180023.50 综上所述，直杆所处 minblilbi180023.50lb x2模型的求

ltf,,d,t,拟牛顿算法基本思想拟牛顿迭代法条件用B表示对矩阵H的近似近似矩阵，用D表示对矩阵逆矩阵H1的B DH设经过k1此迭代后得到xk1，此时将目标函数fx在k1附近作 fxf f x 1x T2f x k k k k k k目标函数两边同时作用一个梯度算子fxfxk1Hk1xxk1对上式fxfxk1Hk1xxk1xxk，得gk1gkHk1xk1xkSkxk1xkykgk1gkgk1gkHk1xk1xky SSH1 k1 k1拟条件在迭代过程中的矩阵Hk1进行约束。因此对Bk1近似的矩 。则y S yk算法步

k k k1 k1Step1：给定初值,,h和精度阀值D0I，k=0Step2：确定搜索的方向dkDkgk。Step3：利用kargminlxkdkSkdkxk1xkSkStep4gk1ykgk1gkDk1Step5：令kk1Step2。求解结果根据上述步骤，利用软件编程（源程序见附录5）求解，得到附件1中直杆的位置可能为以下两个地点：51经度位置坐 纬将表格中两个地点位置坐标代入地图[6]，得到两个地点的具体方位如下图图6位置1具体方 图7位置2具体方偏角之差偏角误差分析，通过附件一的的X，Y坐标，拟合出一条曲线，算 arctan 将问题二中求得直杆地点、杆高、以及附件1期、时刻，计算高度角，sinAcos cos s Aarcsin

AA tt

cos 最后以附件1中相邻时刻影长变化曲线夹角之差与对应方位角之差相将,A求得值拟合成一条曲线，计算出两条曲线的相关性如下图-

-- -- -角度变化率角度变化率-角度变化率角度变化率- -- -- - 15.215.4

15.215.4图8位置1角度变化率与时间的关 图9位置2角度变化率与时间的关实际影长与求解影长的对比与分2的影长相对比如下表：61、2与实际影长的相关性，得到相关6.5.模型的检本文将文献[7]中的数据（详细请见附录6）检验问题三的模型。将建筑模型检验结果表示：7118046'32003' '104606.6结果分析表格，运用软件编程（源程序见附录7）求解对单目标优化模型求得的影长与附件1中实际值的差值对比曲线如下：2

5x10x10A地所长度的x10B地所长度的43差差差差210

0 时 时10由上图可得，误差精度为103。位置1中所计算的影长与实际影长的差值最大为2.4103，最小为1.4104，其差值的规律呈先减小，后增加的规律；位置2中差值最大为4.8103，最小为2104，差值的变化规律先减小后增加。问题三的解答双目标优化模型的建立确定目标函数minlll 确定约束条件

minbibi

180023.50综上所述，直杆所处地点双目标最优化模型如下：

minbtb lif,,dt,t,lbi

x2y2s.ttanb 模型的求

180023.50人工鱼群算法流程图开开YYNNiYgenNY结确定最优iiXi 为，得到Xnext2,Ynext2Xi聚群行为，得到Xnext1,Ynext1在给定范围内初始话鱼群X1,X XN,gengenXiXnextXi11求解结果运用上述建立的双目标最优化模型，并依据题中附件2、3中顶点坐标数据，利用软件（源程序见附录9）求得附件2、3中直杆所处位置及82、3杆长（米 71011252、3中直杆具体所处地理位置，本文根据上述表格中两图12附录2具体方 图13附录3具体方结果分析东经78.9°北纬41°点影长与时间的关的长度的长度112.712.812 13.113.213.313.413.5时间14误差分析式中w1

w2—实际

ew1-件2、3中实际影长的相对误差值如下表：95.4910-1.092510-由上表可得，双目标优化模型影长与附件2实际影长相对比，误差为5.4910-43实际影长相对比，误差为1.092510-4，问题四的解答模型的准备图像预处理（1） 像进行灰度化处 Grayi,j0.299Ri,j0.587Gi,j0.144Bi,j（2）图像进行二值化处理图像呈现出明显的黑白效果。定义一个二值化函数fx',y'，为阀值。（3）小连通域处理

fx,y

fx',y'fx',y'Canny边缘检测算法Step1：对图像进 滤 hx,y,

x2y22gx,y为平滑后图像，用hx,y,对图像fx,y的平滑可表示为gx,yhx,y,fx,y式中Step2：计算梯度幅值和方向已平gx,y的梯度，可以使用22一阶偏导有限差分近似计算x和Y两个偏导数阵fxyf'xy，即：x f'x,y f'x,y

fx1,yfx,yfx1,y1fx,y2fx,y1fx,yfx1,y1fx1,y2 x, 2 x, 22x x,yx,yarctanGxx,yGx,yyMxy—图像的边缘强度x,y—边缘的方向。使Mx,y取得局部最大值的方向x,y，反映了边缘的方向。yStep3：非极大值抑将梯度角离散为圆周的四个扇区之一，以便用33的窗口做抑制运算，4个扇区的标号为0到3，对应33临域的四种可能组合，在每一点上，临域的中心Mx,y与沿着梯度线的两个像素相比Mx,y的梯度值不比梯度线的两个相邻像素梯度值大，则令Mx,y0。Step4：双阀值算法检测及连接边首先对非极大值抑制图像作用两个阈值k1和k2，两者关系k10.4k2把梯度值小于k10，1k2的像素灰度值设022连接边缘步骤对图2pxy时，px,为开始点的轮廓线，直到轮廓线终点qx,y1与图2qxy点位置对应sx,y8邻近区域。如sx,y点的8sx,y存在，则将其包括到图像2中，作为rxy点。从rxy12当完成对包含px,y变换C c,C 23C—变换矩阵，u,v为原始的左边，对应得到变换后的坐标xyx',yx

c,cC可以看做由4

c,c用于平移，c,cT产 变换31

c, 22 cucvx

u

v

cucvy

' c13uc23vAD为Y4A0, Bb, Cc, D0,d，记4个点对 变换的坐标A'xy B'xy C'xy Dxy GGFADy x图15中直杆底 y0 x ccc y ccc1 13

1

13 x2

c11+c21c31c13x2c23 y2

c12c22c32c23y2c23y2x ccc y cc3 23 3 23 计算变换矩阵 1

x1 x

y1

y1 y 2 2 x0x1x2x3 y0y1y2y3c11 x1x0c12x1x3y3x1x3y3x1y1x1y1x1y1

12231 c

yycy c22

yycy 233 利用软件编程（源程序见附录10）对处理过程如下图彩色图片灰度图片二值图片二值降噪图片16模型的建立（1）确定目标函数minll 综上所述，建立的单 minlblilbilif,,d,t,s.t180023.50模型的求解拍摄地点的求解Step2:使用MATL的fmincon的内置函数工具箱，函数形如lfmincon,h3米，为39054'26''N,116023'29''E作为初始值差，差值l最小时该组可能为直杆所在地点。根据上述步骤，利用软件编程（源程序见附录）求解，得到附结果分析2度2度的长2的长2

影长 时间的关1

9 时间17的地理位置处于（东经109°，北纬40.1°），利用地图对查询得到具体地理位置处于内附近。拍摄地点与日期模型建立 法得出拍摄地点为11006'E,39036'N和日期为2015年7月19日。杆高敏感度分析直杆高度变化对的影98差7差6的误5的误4321 直杆的高度差18在此，使杆高在1.6mh2.3m变化，通过观察上图可得直杆位置变化较明模型的评价、改进及推广模型的评价模型的优点问题一中所建立的长度变化的模型，选取直杆所在经度、纬度、日模型的缺点模型的改进模型的推广参考文献 节气太阳垂直照射地球的纬度变化.能热发电技术.中国电.能热发电技术.中国电.20135月.拟法及其收敛性[D].湖南大学地 .利用视坐标的计算进行物高测量与定位[J].林业大学学(自然科学版),1991,03:89-史峰 智能算法.航空航天大学.2011年7 数字图像处理与识别系统的开发[D].东南大学加泰罗尼亚.Canny边缘检测算法原理及其VC实现详解)..2015913日何援军.和投影变换——论图形变换和投影的若干问题之三[J].计附t=%%日期的处理及赤纬角的计算 %n为一年期的序%%时间的处理及时角的 %t0为当地时间 %%求解高度角%%求出最低点的坐标%%画出变化曲xlabel('时间');ylabel('的长度');holdonaxis([9153误差检%plot(x,y-0.4,'r %画出数 %画出数据% 计算得到的曲线','误差域%%%画出变化率曲%%%%xlabel('时间%%%%%%求出函数表达t= %t为时fai %fai为纬seita date= %date为一年期的序h= h为杆%时间的处理及时角的计 %t0为当地时间 %日期的处理及赤纬角的计算%求解高度globalL;%%求出影长与日期的的关t= %t为时fai %fai为纬seita seitadate= %date为一年期的序h= hxlabel('日期序号');ylabel('的长度xlabel('日期的变化');fori=1:size(d,1) %n为一年期的序%%求出影长与直杆长度的关t= %t为时fai %fai为纬seita seitadate= %date为一年期的序h= h为杆xlabel('直杆的长度');ylabel('的长度%直杆长度的变化对影长的影响xlabel('直杆的长度的变化长度');ylabel('影长的变化率%%求出影长与经度的的关t= %t为时fai %fai为纬seita date= %date为一年期的序h= hxlabel经度:60°Eylabel('的长度%经度对影长的影响l3=diff(L3,seita);xlabel('经度的变化');ylabel %%求 影长与纬度的的关t= %t为时fai %fai为纬seita seitadate= %date为一年期的序h= hxlabel纬度：30°S30°N');ylabel('的长度');%纬度对影长的影响l4=diff(L4,fai);xlabel('纬度的变化');ylabel %%求 影长与时间的的关t= %t为时fai=%fai为纬seita=date=date为一h=h为杆title('影长与时间的关系xlabel('一天内的时刻');ylabel('的长度');%一天中时间变化对影长的影响l5=diff(L5,t);xlabel('一天中时间变化');ylabel('影长的变化率[a,b,c]=xlsread('附件1-3.xls','附件[~,~,n]=xlsread('附件1-3.xls','附件%%根据数据拟合出p=polyfit(x,y,2);y1=polyval(p,x);%%计算出影长fori=1:length(c)%%日期的处理 %n为一年期的序%%寻找forforforcal=(h.*(1-(sin((469.*pi.*sin((2.*pi.*(date+cos((pi.*(15.*abs(seita./15+t-8)-180))./180).*...

error=sum(abs(cal-l)./l)/length(l);if

检验所求得下影长与实际影长的关global%第一个fori=1:21xlabel('时间');ylabel差值%%title('东经109.2°北纬19.2°点影长与时间的关系 %hold%%legend('所求得地点的影长与时间的关系曲线','题目数据拟合出的影长与时间%%%%%第二个fori=1:21xlabel('时间');ylabel差值%%title('东经95.4°北纬25.2°点影长与时间的关系 %hold%%legend('所求得地点的影长与时间的关系曲线','题目数据拟合出的影长与时间曲线',0);%%%%长间度高=- %n为一年期的序%%寻找forforforcal=(h.*(1-(sin((469.*pi.*sin((2.*pi.*(date+cos((pi.*(15.*abs(seita./15+t-8)-180))./180).*...

error=sum(abs(cal-l)./l)/length(l);if

[a,b,c]=xlsread('附件1-3.xls','附件%%根据数据拟合出%附件y1=polyval(p,x1);plot(x1,y1,'*r',x1,y1,'-legend('题目所给影端的位置','拟合得到的影端曲线title('影长变化曲线xlabel('X(米计算出影长与时间fori=1:length(c)%时间与影长的拟z1=eval(sym('0.0982*t1^2-%求出最低点的坐标forforfor cos((pi.*(15.*abs(seita./15+t-8)-180))./180).*... 284))./365))./3600).*sin((pi.*fai)./180)+...cos((pi.*fai)./180).*cos((pi.*(15.*abs(seita./15+t-8)

error=sum(abs(cal-l1)./l1)/length(l1);if

xlabel('时间');ylabel('的长度');holdon[d,e,f]=xlsread('附件1-3.xls','附件%%根据数据拟合出%附件y2=polyval(p,x2)