2012年全国数学建模竞赛优秀选

题 要添加到了画出的云图中。叉相关系数法将两时刻的云图进行匹配，求解21：00时刻风矢量。最后利二中的固定窗口模型相比，辨识速度提高了48.3%。 问题重 符号说 问题的假 问题分 基于云图的风矢场（云导风）度量模型的建立和求 视场坐标系的变 帧坐标映射到模型的建 映射到帧坐标模型的建 模型的求 风矢场的度量模型的建立和求 模型1：基于灰度值相关系数的风矢场度量方 模型1的改 模型中一些问题的讨 窗口大小和搜索范围自适应算 窗口大小的调 搜索范围的调 风矢量所在等压面模 模型的建 模型的求 结果分 结果分 模型评价 模型优点 模型缺点 参考文 附 云图在掌握大气环流、中长期天气预报以及性天气学的研究中有转换关系，k_temp.txt102401023所K-1时对应的是地球以外的探测点。[注：地球是(南北方向)和行扫描角（东西方向）均为140微弧（1弧度＝ 是5公里。本题提供的探测数据文件都是2288×2288的灰度值矩阵，矩阵的每个元素都对应地球上或地球外的一个探测点（或称采样点。同步离地 米,星下点在东经86.5度,北纬0度，星下点对应的矩阵元素位于矩阵的第1145行和第1145列相交处。 .5m,短半轴 南、左西右东的地图规则换算成地球上坐标，此结果既可用于估算各探测立矩阵形式的坐标文本文件，这里矩阵的第i行与第j列，分别对应灰度件取名为jwd.txt；而在中给出第500行和第、、列三个交叉点处对应的。进而将灰度矩阵IR1_2100.mat转化为云图（参考附录中的云图，并在这个云图上添加海岸线。保存云图为PDFJPEG文件，取名为spic.pdf或spic.jpg。IR1_2130.mat16×16个像素，搜索2中窗口大小和搜索范围的限定，分析比较应用附件中全球温度数值预报文件temp3.mat计算问题2中非零风矢所在的等压面（用该时刻该处大气的压强值标志，物理单位为毫巴。并在中给W J h alphaXOY 扫描线与YOZ NJ ONXOY ON与YOZ 地球上表面某点N的X轴坐标 地球上表面某点N的Y轴坐标 地球上表面某点N的Z轴坐标 R1,R2,R3fffhd(a,bld(a,b：同一列所对应各行的相对变化率LxLy：像素块在x,y方向上的运动分量PXZ,YZ)灰度图像中点的XZ,YZ灰度值CSP10视为一个点相关系数的风矢场度量方法[1]。选取21:00时刻云图中某一的点，得64×64像素块范围内的搜索。比较的依据为目标块灰度值和合理的。对于有些云，16×16的像素块内云的灰度值分布较广，即方差较大，这样就会使得16×16的像素块不是按照一致的矢量进行移动，如果用16×1664×64个像素的范围内进行搜索，范围相对较大，这会增加计算机的计处温度与气压层之间的对应关系，所处气压层和压强的关系。由于给出的关系确定某一对应的气压层值。系，如图1所示。O

(X,Y 1以地球的中心O为原点，以地心与所在直线为Y轴，地球的长半轴为Z轴，经过原点O且与平面OYZX轴建立笛卡尔直角坐标系。其中Y轴正方向为从地心到所在点，Z轴正方向为从地心指向北极点，X轴正方向为从北极上方看，逆时针旋转与Y轴正方向成90度的方向，如图2所示NNx0,y0,z0Sd0,RL2由以上建立的坐标系可知，地心的坐标为A1(0,0,0)卫星的坐标为P(0h,0h 标(1145,1145)。地心地球可以视为理想椭球体，已知此椭球体的长半径为RL x2

矩阵中的一点MXZ,YZ)。此点与星下点Sd的坐标距离为LXYZ1145LYXZ1145。 按步进角(南北方向)和行扫描角（东西方向）均为微弧采样，所以也就可知扫描此点M时，扫描转过的角度：步进角度（

(140/106)(180/ 扫描线与平面YOZbetaLY(140/106)(180/ 当beta0时：x

(yh)k1(y

当beta0

x ztan(alpha)(yh)k(y bbb24aycyyy y

bbb24aycy当beta0时 y2

a 2k2k2R2R2ka

b2

k2k2R2)

L当beta0L

cy

k2k2R2)h2R2R2k ayk2R2b2k2

cyk2h2 可知点N的坐标为：当beta0y 22当beta0

y 22 地心连线ONXOYJX和平面YOZJZ，如图2所示。x2x2

arctan(x0

星下点所在位置为东经86.5度,北纬0度，则东经86.5度线在平面YOZ根据的定义可知：N所在纬度W

WJX J86.5JZ 以上联立就得到就可以得到帧坐标系中点M与三维坐标系中对应点N5.1.2映射到帧坐标模型的建N的经纬坐标(WJ，则ONXOYJX和平面YOZcos(Jcos(J)sin(J)2 (cos(JX)RS)2(sin(JX)RL2cos(J)cos(J) (cos(JX)RS)2(sin(JX)RL2sin(J) (cos(JX)RS)2(sin(JX)RLx0y0z

由此可知扫描线与平面XOY的夹角alpha和扫描线与平面YOZ的夹角betaalphabeta

h h

LX

LY (140/106)*(180/(140/106)*(180/ 1要求，计算出了灰度矩阵中的第451行至550行，第451列至550列的点所对应，存入了jwd.txt中。利用MATALB软件将灰度矩阵IR1_2100.mat转化为云图，并利用所oX.m给出的海岸线所对应的帧坐标值，并将其添加到了云图中，如图3所表1给定点图3标注海岸线的云的坐标变换可以得到其对应的帧坐标的位置，即对应扫描的行列值分别为ij，以(ij16×16（本文选取像素块中坐标为(i,j)的像素点，即所求对应的像素点）为中心，在21:30时刻的云图中分别向上、向下、向左、向右平行移动32个像64×64也就是说，21:00时刻的目标块经过半个小时移动到了搜索到的匹配块，则两图对应块中心的移动矢量即为所求的风矢。发生变化，在理想情况下会变为一个定值，如图4所示。

4经过对图像的像素数进行统计，发现此云图中的灰度值范围在300至p(x)1H

pr(r)

p(s)p(r) p(r) 对等式两边对r

1rsf(r) p()dr

rr0Df(D) 对于离散图像，转换

DBf(DA) max A05656将所给三个时刻的灰度值用作图，从图中可以看出，存在明显的0.5时，20:3021:0020:3021:00时刻的图中噪点个数基本不变，21:30不变，21:30时刻的图中发现了17000多个噪点。从中可以得出，同步扫描过程中的硬件误差和传感过程中的软件误差等原因，会产生噪点，21:30时刻相较少的改变原始数据，本文选择阈值为0.45，对噪点进行判断并处理。

[f(i,j)f][g(im,jn)j[[f(i,j)f]2[g(im,jn)g]2

其中，ij[1,16]分别为目标块的行与列号，mn[0,64]为搜索区域的位移改变量，f(i,j)为目标块云的灰度值，g(i,j)为搜索区域云的灰度值，f1f(i,j为目标块云的平均灰度值，g1g(im,jnN2 N2 块相对于目标块在x,y方向上分别移动了m和n个像素点。21:002050°为例，所得到的最大交叉相关系数值为0.95，在21:30时刻的云图上位于南纬20.04°东经49.54°。40度至北纬40度，东经46度至126度区域的整数处的风矢，则起始点px2/re2y2/r2 p处球的半径，且re= .5m,rp= .8m,那么地球的偏心率为：21(r/r)2 若子午线上任意一点的地心纬度为r(x2y2)1/2rp[(1tan2)/(1tan22 再计算云矢的大小和方向。示踪云目标的初始位置A(AA)F7在 中，ONOAOBOEOFR NOA

NOB90B,NEOFBA,AOBA=球面角NAB由球面三角形NABcosABr

fffr/ t为第二张图扫描到BA点的时间之差。再由球面三角形ABN的余弦定理：B)coscos(90A)sinsin(90A)cos

cosA(sinBcossinA)/(sinAcos1[(sinBcossinA)]/(sincos

fengshi.txt2652，53，54，57，58（度）处表2几个具体处的风矢（相关系数法角度（北顺速度（米/秒00的原始位置，计算预测位置与目标块之间灰度值的相关系数，并求出均值为f(x在区间上给定n1个插值节点：ax0x1xnb及其在插值节点的值f(x0),f(x1f(xn)，若函数S(x)满足S(x在每个小区间[xjxj1j0,1,nS(x)在每个小区间[abS(xj)f(xj)，j0,1,,nS(x)为f(x的三次样条插值函数。三次样条插值函数的具体求法如下三次样条函数的二阶导数S(xj)在每个小区间[xj,xj1S(x)

1

M)(xx)(x[x,

j

jS(x)

{f[x, ]

1

)h}(xx) j

(3M

j 1M(xx)2

1

M)(xx

j (x[xj,xj1],j0,1,...,n

S(xS(x)求一阶导数，分别计算每一点在帧坐标的关系数分别记为R2,R3。[fx(i,j)fx][gx(im,jn) [[fx(i,j)fx]2[gx(im,jn)gx]2 [fy(i,j)fy][gy(im,jn)

[[fy(i,j)fy]2[gy(im,jn)gy]2

f(i,jg(i,j ,, 本文取k10.4,k20.3,k30.3 R,R 3J[m,n]k1R1(m,n)k2R(m,n)kR(m, 3本文取k10.4,k20.3,k30.3计算预测块与目标块灰度值相关系数的均值为0.4757。表3几个具体处的风矢（差商法角度（北顺速度（米/秒00其中k1k2k3k4

0.4791。表4几个具体处的风矢（整体局部匹配算法角度（北顺速度（米/秒005度高于摄氏0度的位置作为无云区域，显然此点的风矢为零风矢。20:3021:0021:0021:3021:0021:3020:30本文采用计算相对变化率的方法得到窗口大小，设计思想为：对于一个果在N×N的像素块中，存在着相邻两行或者两列的灰度值变化很大，可以认为先将所要度量风矢的转换到帧坐标中，对应坐标为(x,y),将其扩展为后一列灰度值的相对变化率的绝对值。如下： b huidu(a,hd(a,b)huidu(a,b1) b huidu(a,度值的相对变化率的绝对值。如下：huidu(a1,b)huidu(a, a huidu(a,ld(a,b)huidu(a1,b) a huidu(a,列的灰度值以及第jn+1jn列的灰度值变化率最大。）8N×N绝对值最大点的坐标分别为(p,r)和(q,s)p-1rpr列的灰度值以及第q+1行第s列对于第q行第s列的灰度值变化率最大。（（9N×Nijr列和第s列所围成的区域，确定的窗口大小为(j-i)×(s-r)的像素块。表6几个具体处的风矢（自适应窗口角度（北顺速度（米/秒00BBCAD10B、C所对应的运量分别为(MV_A_x，MV_A_y)、(MV_B_x，MV_B_y)、若当前块的左边、上边、右上边的云矢未知，则先将搜索范围设定为度。计算如下：n2n2L n2n2 ,Ly n3MVxMVy分别为MVxyLxLy SxmaxMVx SymaxMVy

Smax(Sx,Sy

xx

3MVxix

3个yiyW S

WS Stemp3.txt中任意气压层对应的为二维矩阵，同层矩阵行坐标XW与纬度W的关系为： (90W)

Y

灰度图像