研究生数学建模真题解析2012D

云导风(基于卫星云图的风矢场度量模型与算法探讨)2/2/2023江惠坤南京大学机动目录上页下页返回结束——2012年全国研究生建模竞赛D题解读原题一、问题背景

一、问题背景

2012年全国研究生“华为杯”数学建模竞赛D题“基于卫星云图的风矢场（云导风）度量模型与算法探讨”，是根据大气科学和应用领域关于“云迹风”计算的一个科研问题改造而来的。初看本赛题好像比较专业，其实我们只要理解了数据来源和格式的说明，利用常识而不需要大气科学的专业知识就很容易理解赛题的问题，赛题的问题就变成了数据处理的问题，所需要的方法也只是涉及数学和电脑应用的方法。问题的重点是要探讨出比较好的适合计算云导风的方法。上页下页返回结束二、问题重述二、问题重述

本题共提供4种共6个数据文件，它们分别是：1、k.temp.txt（灰度值—温度值的对照表）2、IR1.2030.mat，IR1.2100.mat，IR1.2130.mat.（云图每个像元的灰度数据）3、coastline0.txt（海岸线的经纬度坐标）4、temp3.mat（21:00时的温度数值预报）上页下页返回结束问题1分为两个小问题：第一个小问题是视场坐标的换算，是要将灰度矩阵元素行列号（2288x2288，除去元素值为-1的那些行列号）换算成地球上的经纬度坐标，这只是一个空间解析几何的问题。第二个小问题是根据所给的数据文件（coastline0.txt）展示具有海岸线的的卫星云图，这个小问题只是数据可视化问题，是有关应用软件的一个应用。二、问题重述二(1)问题1

上页下页返回结束二、问题重述问题2是本赛题中的主要问题：题目提供了三个相邻时刻的卫星数据：IR1_2030.mat,IR1_2100.mat,IR1_2130.mat要求用适当的方法建立21:00时刻风矢场的度量模型，给出计算方案，并在限定的范围内给出计算结果。也就是要给出计算云迹风的模型、方法和部分结果。

二(2)问题2

上页下页返回结束二、问题重述问题3是要：

取消问题2中窗口大小和搜索范围的限定，利用图像特征匹配等思想，设计能自适应地确定窗口大小和搜索范围的有效方法。进而建立风矢场的度量模型和算法。二(3)问题3、4

问题4

是要：

依据全球温度数值预报的信息（temp3.mat

）去考虑云迹风所在的等压面高度（用压强表示这个高度）。（这里要用到数据文件k.temp.txt

提供的温度与灰度的对照表）上页下页返回结束三、问题1的参考解答与评述问题1-1

问题1的第一个小问题是视场坐标的换算，是要将灰度矩阵元素行列号换算成地球上的经纬度坐标，这只是一个空间解析几何的问题。这个小问题是为计算云迹风做准备的，是本赛题的基础性工作，它虽然不是本赛题的重点问题，但仍然是非常重要的工作，这里的结果直接影响到后续工作的结果，对待这样的问题需要的是认真和严谨，不能在起点上就跌倒。但本届参赛队在这个问题上得到准确结果的比例很低，还不到20%，远低于我们的预期。出现的错误五花八门，很多队的结果常常相差几十个经纬度，有的甚至分不清南纬北纬、东经西经的意义，这些带有常识性的错误完全是可以避免的。上页下页返回结束三、问题1的参考解答与评述问题1-11、视场转换的数学模型，通常有直角坐标系下方程组和球坐标系下方程组两种。直角坐标系（A）直角坐标系：其中为扫描线南北方向按步进角与xoy-平面所成夹角，为扫描线东西方向按行扫描角与yoz-平面所成夹角，l=42164000m为卫星到地心的距离，a=6378136.5m为长半轴，b=6356751.8m为短半轴；A(x,y,z)为球面上的探测点，满足方程：从这个方程组求出x,y,z.再根据求出经度和纬度上页下页返回结束三、问题1的参考解答与评述问题1-1（B）球面坐标系：Xz在上图中将zOx平面平移至原点与卫星重合，组成Oxyz以卫星为原点的直角坐标系，地心的坐标为（0，-l，0）。

在这个坐标系中建立球坐标，设r为探测点A到卫星的距离，分别表示球面坐标的和角(于是有设则探测点A(x,y,z)满足),上页下页返回结束问题1-1

三、问题1的参考解答与评述，由此解得r满足的一元二次方程，从中，经度为-3.5，纬度

下计算的结果完全一样。在指定点的结果为：解得r，又得到x，y，z

的值，再计算。计算结果与直角坐标500行500列经纬度：(46.49,32.74)500行501列经纬度：(46.57,32.74)500行502列经纬度：(46.65,32.73)上页下页返回结束三、问题1的参考解答与评述问题1-2

（2）添加海岸线这第二个小问题是根据所给的数据文件展示具有海岸线的的卫星云图。这个小问题只是数据可视化问题，是有关应用软件的一个应用，虽然不是数学建模的核心问题，但仍然是数学建模所经常要遇到的问题，特别对数学建模的完整和完美性、对展示建模的结果有着重要作用。本届参赛队对这个小问题的完成结果是令人满意的。上页下页返回结束三、问题1的参考解答与评述问题1-2

（2）添加海岸线设A(x,y,z)为球面上一点，对给定海岸线上点的经度和纬度，由解出(x,y,z),再由解出这样就可以在卫星云图上添加海岸线的点了。图示为：

上页下页返回结束问题1-2

三、问题1的参考解答与评述上页下页返回结束问题2-1

四、问题2的参考解答与评述问题2是本赛题中的主要问题。题目提供了三个相邻时刻的卫星数据（IR1_2030.mat,IR1_2100.mat,IR1_2130.mat），要求用适当的方法建立21:00时刻风矢场的度量模型，给出计算方案，并在限定的范围内给出计算结果。也就是要给出计算云迹风的模型、方法和部分结果。这个问题在大气科学和应用领域已经被广泛研究过，有些结果也正在被应用着。但由于探测手段的局限性和云系的不稳定性，现有的方法总是有不尽人意的地方。探讨出更多更好的方法是非常有必要的。上页下页返回结束四、问题2的参考解答与评述问题2-1

这个问题在本质上可以看成是图像匹配的问题，匹配的过程大致可以分为4步：第1步，对数据的去噪处理；第2步，云块匹配；第3步，风矢的计算；第4步，云检测与质量控制。下面对这4步的要点和主要问题分别给出参考解答并给予简要评述，一家之言，仅供参考。上页下页返回结束四、问题2的参考解答与评述问题2-1

先看第1步，由于云块不是刚体，它的形状会随时间不断变化，加上红外探测仪的测量数据也会有一定误差，所以对数据进行适当的去噪是必要的，去噪可以消除一些随机误差和由微小变化产生的误差，可以给后面的匹配消除一些干扰。去噪的方法有很多，包括中值滤波、高斯滤波、低通滤波等，都是很好的去噪工具。而做本题的参赛队中只有很少的队对数据做了去噪处理。当然，如果下面采用图像的频谱特征进行图像匹配，可以将滤波去噪的过程放到匹配过程中一并进行。

上页下页返回结束四、问题2的参考解答与评述问题2-2

对于第2步的“云块匹配”，首先在指定经纬度处，在21:00的卫星云图中确定一幅16x16的云块图像A，它对应灰度矩阵中的一个16x16的像素块，由这个像素块上的二维数据表征这幅图像。由于指定的经纬度都是整数，而这些整数经纬度结点一般不会对应云块图像上的某个像素点，但可以对应相邻4个像素点组成的像素块，以这个像素块为中心的一个16x16的像素块就可以作为该指定经纬度处的云块图像A.接下来的问题是要在另外两个不同时间的卫星云图中去匹配指定的图像，不妨在21:30的卫星云图中匹配的图像记为B，在20:30的卫星云图中匹配的图像记为C。而匹配图像要从两个方面入手，一是用图像的什么特征去匹配？二是用怎样的数学方法去匹配？这就是度量模型和算法问题。

上页下页返回结束四、问题2的参考解答与评述问题2-2

图像的特征有很多，如灰度特征，即全部灰度值，是图像的全部信息，其实不能算作真正意义上的特征；

梯度特征，它可以突出反映图像的纹理和边缘信息；

频谱特征，它可以区分信号在各个频段上的能量信息；甚至，灰度的均值、方差、平方和都可以作为图像的特征。至于数学方法，通常采用最大相关系数法和最小差值平方和（或平方）法。对于特征数少或只有单个值的图像特征（如均值、方差），不能采用最大相关系数法去匹配图像，通常可采用最小差值平方法去匹配图像。目前在气象工程中和大量文献中大多采用基于灰度特征或梯度特征的最大相关系数方法，而做本题的参赛队大部分也是采用这种方法，缺少自己的思考和创造性，但也有少数队用到了频谱特征和均值、方差等特征。尽管由于其他原因，其中有些队还是没有得到比较好的结果，但善于独立思考、发挥自己的创造性的精神是值得提倡的。

上页下页返回结束四、问题2的参考解答与评述问题2-3

到了第3步，要根据第2步匹配好的图像去计算风矢，风矢要求有起点、方向和大小。自然地，将图像的中心点看作云块的位置，设匹配图像A、B、C的中心点分别为a、b、c，那么a点就可以作为21:00时刻风矢的起始点，a点的经纬度坐标一般不是整数，这就是在题目的示例中，在整数经纬度处的风矢的起始点坐标都有2位小数的原因。许多参赛队都理解了这2位小数的意义，给出了正确的风矢起点的表达，但也有很大一部分参赛队没有考虑风矢起点的实际意义，不加说明就用整数经纬度表达风矢起点，这就有所欠缺。上页下页返回结束四、问题2的参考解答与评述问题2-3

至于21:00时刻风矢的方向和大小，由于风矢的精度不高，不妨把以a为起点、b为终点的平面向量作为所求风矢，即该风矢的的大小为两点之间的距离，方向的确定按题目的约定：由a点向正北方作一射线，从该射线沿顺时针方向到ab线段的角度即为所求风矢的方向。当然，将a点到b点的球面距离作为风矢的大小是更准确一点，但计算上会复杂许多。同理容易确定以c为起点、a为终点的风矢，这一风矢在下一步的质量控制中要用到。

上页下页返回结束四、问题2的参考解答与评述问题2-4

再进入第4步。由于云块随时间会变化，上面第2步的匹配很容易出现误配，为了减少误配，必须对上2步的结果进行质量控制的工作。首先要进行云检测，设置合适的阈值将无云或极少云的区域赋予零风矢，从而剔出讨论的范围。对第3步的结果进行质量控制就是利用风矢的时间连续性和空间连续性来检验其结果是否合理。应用时间连续性就是看与的差的大小与或的大小是否相差太大，与的差的大小大于阈值的就要剔除或重新计算；应用空间连续性就是检验每个风矢与其周围风矢在方向和大小上是否在一定的范围之内，例如方向控制在正负30度，大小控制在正负1.5倍。

上页下页返回结束四、问题2的参考解答与评述问题2-4

为了避免剔除的风矢过多，我们可以对第2步做修正，取消图像匹配最佳的限制，采用近似匹配，选取若干个匹配对象，例如，在21:30的卫星云图中匹配的图像选出B1,B2,B3，在20:30的卫星云图中匹配的图像选出C1,C2,C3，然后应用时间和空间连续性去选择最佳的匹配，这样做的效果会更好。要得到风矢的数值表示可参考下图：

卫星W设初始位置A

和终点位置B

是A、B的纬度和经度坐标，由上页下页返回结束四、问题2的参考解答与评述问题2-4

就容易计算出该风矢的大小和方向。在指定处计算的风矢的参考结果为：

注：1.由于技术性处理对结果有很大差别，所以允许参赛队的结果与此结果有较大的误差。2.非零风矢的个数不作为好差的依据。3.风矢的起始点的经纬度应该是非整数的上页下页返回结束四、问题2的参考解答与评述问题2-4

关于对窗口大小和搜索范围的限定的理解或评价：根据探测点的密度去估算16x16大小的窗口对应云块长宽各为一个经度或纬度稍多一点，而64x64大小的搜索范围可保证每秒60米多的风矢都在搜索范围之内。上页下页返回结束五、问题3的参考解答与评述问题3-1问题3是要：取消问题2中窗口大小和搜索范围的限定，利用图像特征匹配等思想，设计能自适应地确定窗口大小和搜索范围的有效方法。进而建立风矢场的度量模型和算法。这是一个比较开放性的问题，只要不违反常识并遵循一定的原则的方法都是被认可的。例如，确定窗口大小需遵循的原则应该是：在云块灰度值变化程度较大的地方，选择的云块窗口就可以小一些；反之，窗口就要大一些。其中，灰度值变化程度可以用合适的指标来刻画，例如，均方差、梯度能量、信息熵等。搜索范围的原则应该是：云块在半小时之内按最大风速所能到达的范围。

上页下页返回结束五、问题3的参考解答与评述问题3-2

说到最大风速，同学们在做本题时对最大风速应该要有一个大致的概念，这也是很容易做到的事情。例如，笔者在网络上搜索到的一段文字：“世界上风力资源最多的是世界第七大陆——南极洲。在那里，一般的风速是17-18米/秒，最大风速可达75米/秒以上。但是，南极洲的最大风速并不是世界的最大风速。我国海南岛琼海县，1973年4月，一次台风过境时，估计风速有81米/秒。日本的室户岬，1961年4月的一次台风过境时，风速曾达到84.5米/秒。1934年，美国华盛顿测到了103.2米/秒的风速，这才算是世界上最大的风速。”[1]这段文字就给我们提供了关于最大风速的足够的信息。加上最大风速不可能持续很久、云块的移动速度又比风速要慢，云块越大，移动速度就越慢。所以在考虑搜索范围时，将最大风速限制在70-80米/秒也就可以了。可是在这次竞赛的论文中，有些风矢大小达到每秒7、8百米，甚至还有每秒3千多米的。这种违反常识的错误实在是不应该出现啊！上页下页返回结束六、问题4的参考解答与评述问题4-1

问题4是要依据全球温度数值预报的信息去考虑云迹风所在的等压面高度（用压强表示这个高度）。（1）求风矢所在等压面基本方法：对给定的经纬度(latitude，longitude)，求出灰度矩阵中对应（接近）的行和列(i,j)，得到该点处的灰度值hui(i,j)，根据k_temp.txt将灰度值转换成温度值，与temp3.mat中对应经纬度处上方的36个预报温度值进行比较，确定云块温度介于两个预报温度值之间，于是确定该处风矢处于对应的两个等压面之间，用线性插值的方法求得该风矢所在的等压面。上页下页返回结束六、问题4的参考解答与评述问题4-2

大家在做这个问题的操作过程中很快会遇到3个方面的情况要妥善处理：一是云迹风所在经纬度不是直接对应数值预报的经纬度；二是数值预报在同一经纬度上的温度值与高度值在全局上不是线性关系，从而使得从两个途径得到的温度值在全局上不具有一一对应关系；三是云迹风处由灰度值转换来的温度值不会直接对应数值预报在相近经纬度处的某个温度值，也即该云迹风不在任数值预报中的任何一个等压面上。这些问题要由学生自己去发现和解决。上页下页返回结束六、问题4的参考解答与评述问题4-3

对第一种情况，可以在云迹风所在经纬度的周围附近将数值预报的温度值做一个平均，将这个平均值作为在该云迹风所在经纬度上修改的数值预报温度值，这样，在每个云迹风的经纬度的每一层上（共36层）都有一个修改的数值预报温度值；对第二种情况，在百度百科上分别搜索一下“云层”和“大气层”，那里提供的知识清楚地告诉我们：云层在地球上空10公里以内，而10公里以内的等压面在200毫巴至1000毫巴之间。在这个范围内，数值预报的温度值是随海拔高度单调减少，也就是随等压面的压强值（单位：毫巴）单调增加的，只要考虑在200至1000毫巴等压面之间的数值预报温度就可以了.上页下页返回结束六、问题4的参考解答与评述问题4-4

取平均值以后的温度与压强之间的关系为：1000200上页下页返回结束六、问题4的参考解答与评述问题4-5

对第三种情况，将云迹风所对应的温度值T与同一经纬度上修改的数值预报的各层（限于200-1000等压面，即就在右边那段单调区间内）的温度值进行比较，让该温度值T介于某两层的温度之间，于是就认定这个云迹风介于这两层等压面之间，用插值公式（一般用线性插值）就可计算得该云迹风所在的等压面。显然，云迹风所在的等压面一般不在数值预报中的等压面上，所以这里做插值是必须的。但遗憾的是，这有少数参赛队做了插值的工作。

按以上方法，求得纬度26度，经度分别是52,53,54,57,58（度）处的等压面的参考结果为：上页下页返回结束六、问题4的参考解答与评述问题4-6

（2）给出云导风图只要给出了像样的云导风图，不要求与下面的图一致。以下是国家气象局专业人员用本题数据所计算得到的云导风全图，供参考：（注意：本题只要求中间一部分的结果）上页下页返回结束六、问题4的参考解答与评述问题4-7

上页下页返回结束谢谢！上页下页返回结束附：2012年全国研究生数学建模竞赛D题1/11页卫星云图在掌握大气环流、中长期天气预报以及灾害性天气学的研究中有重要作用。它由地球同步卫星上的红外探测仪探测地球上空的温度数据再转换成灰度数据制作而成。附件中定标数据文件k_temp.txt给出了灰度数据与温度数据的转换关系，k_temp.txt内有1024个实型数，依次是图象灰度数据为0到1023所对应的K氏温度值，灰度值为-1时对应的是地球以外的探测点。[注：地球是被探测温度的唯一来源，如果天空无云，探测到的温度可以看成是地球表面的温度；在有云层的地方，探测到的温度相对较低，且云层越高越厚温度就越低，探测到的温度可看成云层所在区域的温度]。上页下页返回结束附：2012年全国研究生数学建模竞赛D题红外探测仪扫描采样时，按步进角(南北方向)和行扫描角（东西方向）均为140微弧（1弧度＝1000000微弧）采样。在卫星与地球中心的连线和地球表面的交点（称为星下点）处的分辨率大约是5公里。本题提供的卫星探测数据文件都是2288×2288的灰度值矩阵，矩阵的每个元素都对应地球上或地球外的一个探测点（或称采样点）。同步卫星离地球中心的高度为42164000米,星下点在东经86.5度,北纬0度，星下点对应的矩阵元素位于矩阵的第1145行和第1145列相交处。2/11页上页下页返回结束附：2012年全国研究生数学建模竞赛D题为解答本题，首先要确定灰度矩阵中每个元素对应的采样点在地球上的经纬度。地球可视为理想椭球，这个理想椭球可以由地球的一个经过南北极的椭圆截面绕南北极的连线旋转而得到，椭圆截面的长半轴（赤道半径）=6378136.5m,短半轴（极半径）=6356751.8m；据此就可以将灰度矩阵中非负元素的行列号按上北下南、左西右东的地图规则换算成地球上经纬度坐标，此结果既可用于估算各探测点之间的距离，还可用于在云图上依据海岸线经纬度坐标标出海岸线以方便看图。3/11页上页下页返回结束附：2012年全国研究生数学建模竞赛D题观测大气环流情况的一个方法是在卫星云图上标出风矢。风矢的大小和方向由云块移动的速度决定。风矢与风的速度有所不同，如某个台风中一些区域的风速可达每秒五、六十米，而台风（看作云块）中心的移动速度可能仅每小时十多公里。没有云或云块不稳定处的风矢规定为零风矢，这种用云块的移动所定义的风矢被称为云迹风。气象部门已经有一些方法根据变化的卫星云图计算云迹风，这类方法称为云导风方法。计算云迹风时通常将云块大小限定为16×16个像素，搜索范围限定为64×64个像素。4/11页上页下页返回结束附：2012年全国研究生数学建模竞赛D题本题的主要目的是希望大家充分利用卫星图像数据及其特点建立尽可能准确地描述实际风矢场的度量模型和算法。

题目提供了我国风云2号卫星获得的三个灰度矩阵，IR1_2030.mat,IR1_2100.mat,IR1_2130.mat，分别表示某天的20:30,21:00,21:30时刻红外探测仪探测到的地球上空的温度数据对应的灰度值。又给出了海岸线经纬度坐标数据文件coastline0.txt，此文件的第1列为经度(东经),第2列是纬度(北纬),每一行2个数据对应海岸线上一点，而特大数据（99999.99，99999.99）表示前一曲线已结束，将要开始下一曲线。5/11页上页下页返回结束附：2012年全国研究生数学建模竞赛D题具体要求解决如下问题：1、换算视场坐标。给出灰度矩阵元素行列号对应于经纬度坐标的换算公式，建立矩阵形式的经纬度坐标文本文件，这里矩阵的第i行与第j列，分别对应灰度矩阵的450+i行与450+j列，矩阵元素是（经度，纬度）这种形式的二维数组，给出结果的范围为：灰度矩阵中的第451行至550行，第451列至550列，文本文件取名为jwd.txt；而在论文中给出第500行和第500、501、502列三个交叉点处对应的经纬度。进而将卫星灰度矩阵IR1_2100.mat转化为卫星云图（参考附录中的卫星云图），并在这个卫星云图上添加海岸线。保存卫星云图为PDF或JPEG文件，取名为spic.pdf或spic.jpg。

6/11页上页下页返回结束附：2012年全国研究生数学建模竞赛D题2、根据三个相邻时刻的卫星数据（IR1_2030.mat,IR1_2100.mat,IR1_2130.mat），用适当的方法建立21:00时刻风矢场的度量模型，给出计算方案，这里限定风矢计算中像素块匹配所选用的窗口大小为16×16个像素，搜索范围限定为64×64个像素。利用你们的方法计算21:00时刻在各个（整数）经纬度处的二维风矢。确定一个二维风矢需要4个指标：起始点的纬度、经度，风矢的方向（单位：从正北方顺时针方向量得的角度）和大小（单位：米/秒）。计算范围为：南纬40度至北纬40度，东经46度至126度。在论文中给出上面计算结果中全部非零风矢的个数，并给出在纬度26度，经度分别是52，53，54，57，58（度）处的风矢，这些风矢的表示形式可参照下面问题4中对winsh.txt要求的格式。给出你们对上述窗口大小和搜索范围的限定的理解或评价。7/11页上页下页返回结束附：2012年全国研究生数学建模竞赛D题3、为了提高像素块匹配质量，取消问题2中窗口大小和搜索范围的限定，利用图像特征匹配等思想，设计能自适应地确定窗口大小和搜索范围的有效方法。进而建立风矢场的度量模型和算法。利用此模型和算法计算问题2中在纬度26度，经度分别是52，53，54，57，58（度）处的风矢。分析比较你们建立的这两种风矢计算方法的优劣，结论要有定量数据支撑。提出你们关于云图数据量和云导风风矢场计算研究的建议。8/11页上页下页返回结束附：2012年全国研究生数学建模竞赛D题4、应用附件