基于卫星云图的风矢场度量模型与自适应算法（MATLAB）

目 录

绪 论 … - 1 -
0.1 问题背景… - 1 -
0.2 概念界定… - 2 -
0.3 本文建模假设… - 3 -
0.4 主要变量说明… - 3 -
0.5 技术路线图… - 5 -
第 1 章 卫星视场坐标变换模型… - 6 -
1.1 问题分析… - 6 -
1.1.1 问题描述… - 6 -
1.1.2 坐标系构建… - 6 -
1.2 帧平面到经纬度转换模型… - 7 -
1.2.1 转换模型… - 7 -
1.2.2 模型求解… - 9 -
1.3 经纬度到帧平面转换模型… - 9 -
1.3.1 转换模型… - 9 -
1.3.2 模型求解… - 10 -
1.4 本章小结… - 12 -
第 2 章 风矢场固定尺度度量模型… - 13 -
2.1 问题分析… - 13 -
2.2 目标云定位… - 13 -
2.3 匹配云定位… - 14 -
2.4 风矢度量模型… - 15 -
2.5 模型求解… - 16 -
2.5.1 参数灵敏度… - 16 -
2.5.2 计算结果… - 18 -
2.6 本章小结… - 18 -
第 3 章 风矢场自适应度量模型… - 19 -
3.1 自适应变化风矢场算法… - 19 -
3.2 方法对比及建议… - 21 -
3.3 本章小结… - 21 -
第 4 章 风矢场等压面模型… - 22 -
4.1 问题分析… - 22 -
4.2 等压面的确定… - 22 -
第 5 章 模型改进… - 24 -
参考文献… - 25 -

绪 论

0.1问题背景
卫星云图在掌握大气环流、中长期天气预报以及灾害性天气学的研究中有重要作用。它由地球同步卫星上的红外探测仪探测地球上空的温度数据再转换成灰度数据制作而成。附件中定标数据文件 k_temp.txt 给出了灰度数据与温度数据的转换关系，k_temp.txt 内有 1024 个实型数，依次是图象灰度数据为 0 到 1023 所对应的 K 氏温度值，灰度值为 -1 时对应的是地球以外的探测点。红外探测仪扫描采样时，按步进角(南北方向)和行扫描角（东西方向）均为 140 微弧采样。
在卫星与地球中心的连线和地球表面的交点（称为星下点）处的分辨率大约是 5
公里。本题提供的卫星探测数据文件都是 2288×2288 的灰度值矩阵，矩阵的每个元素都对应地球上或地球外的一个探测点（或称采样点）。同步卫星离地球中心的高度为 42164000 米, 星下点在东经 86.5 度, 北纬 0 度，星下点对应的矩阵
元素位于矩阵的第 1145 行和第 1145 列相交处。
为解答本题，首先要确定灰度矩阵中每个元素对应的采样点在地球上的经纬度。地球可视为理想椭球，可由地球的一个经过南北极的椭圆截面绕南北极的连线旋转而得到，椭圆截面的长半轴（赤道半径）=6378136.5m, 短半轴（极半径）
=6356751.8m；据此就可以将灰度矩阵中非负元素的行列号按上北下南、左西右东的地图规则换算成地球上经纬度坐标，此结果既可用于估算各探测点之间的距离，还可用于在云图上依据海岸线经纬度坐标标出海岸线以方便看图。
观测大气环流情况的一个方法是在卫星云图上标出风矢。风矢的大小和方向由云块移动的速度决定。风矢与风的速度有所不同，如某个台风中一些区域的风速可达每秒五、六十米，而台风（看作云块）中心的移动速度可能仅每小时十多公里。没有云或云块不稳定处的风矢规定为零风矢，这种用云块的移动所定义的风矢被称为云迹风。气象部门已经有一些方法根据变化的卫星云图计算云迹风， 这类方法称为云导风方法。计算云迹风时通常将云块大小限定为 16×16 个像素，
搜索范围限定为 64×64 个像素。
本题的主要目的是希望大家充分利用卫星图像数据及其特点建立尽可能准确地描述实际风矢场的度量模型和算法。
题目提供了我国风云 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）表示前一曲线已结束，将要开始下一曲线。
具体要求解决如下问题：
（1）换算视场坐标。给出灰度矩阵元素行列号对应于经纬度坐标的换算公式，建立矩阵形式的经纬度坐标文本文件，这里矩阵的第 i 行与第 j 列，分别对应灰度矩阵的 450+i 行与 450+j 列，矩阵元素是（经度，纬度）这种形式的二维数组，给出结果的范围为：灰度矩阵中的第 451 行至 550 行，第 451 列至 550 列， 文本文件取名为 jwd.txt；而在论文中给出第 500 行和第 500、501、502 列三个交叉点处对应的经纬度。进而将卫星灰度矩阵 IR1_2100.mat 转化为卫星云图（参考附录中的卫星云图），并在这个卫星云图上添加海岸线。
（2）根据三个相邻时刻的卫星数据（ IR1_2030.mat, IR1_2100.mat,

• 1 -

IR1_2130.mat），用适当的方法建立 21：00 时刻风矢场的度量模型，给出计算方案，这里限定风矢计算中像素块匹配所选用的窗口大小为 16×16 个像素，搜索
范围限定为 64×64 个像素。
利用所给的方法计算 21：00 时刻在各个（整数）经纬度处的二维风矢。确定一个二维风矢需要 4 个指标：起始点的纬度、经度，风矢的方向（单位：从正
北方顺时针方向量得的角度）和大小（单位：米/秒）。计算范围为：南纬 40 度
至北纬 40 度，东经 46 度至 126 度。
在论文中给出上面计算结果中全部非零风矢的个数，并给出在纬度 26 度， 经度分别是 52，53，54，57，58（度）处的风矢。
给出对上述窗口大小和搜索范围的限定的理解或评价。
（3）为了提高像素块匹配质量，取消问题 2 中窗口大小和搜索范围的限定， 利用图像特征匹配等思想，设计能自适应地确定窗口大小和搜索范围的有效方法。进而建立风矢场的度量模型和算法。
利用此模型和算法计算问题 2 中在纬度 26 度，经度分别是 52，53，54，57， 58（度）处的风矢。
分析比较建立的这两种风矢计算方法的优劣，结论要有定量数据支撑。提出关于云图数据量和云导风风矢场计算研究的建议。
（4）应用附件中全球温度数值预报文件 temp3.mat 计算问题 2 中非零风矢所在的等压面（用该时刻该处大气的压强值标志，物理单位为毫巴）。并在论文中给出在纬度 26 度，经度分别是 52，53，54，57，58（度）处附近求出的风矢量所在的等压面。

function draw_f(x,y,vel,ang,pres)
% 画风矢量符号
% 输入： x,y坐标，速度，角度，压强
%%
f_length = 50;
f_patit = 7;
f_width = 12;
ang = ang/180*pi+pi;
%%
if vel<20
    num_bar = 1;
elseif vel<40
    num_bar = 2;
else
    num_bar = 3;
end

%%
if pres<=200
    color_pres = 'r';
elseif pres<=500
    color_pres = 'g';
else
    color_pres = 'b';
end
%%
for i=1:num_bar
    pstart(i,1) = x + (f_length-f_patit*(i-1))*sin(ang);
    pstart(i,2) = y - (f_length-f_patit*(i-1))*cos(ang);
    pend(i,1) = pstart(i,1) + f_width*cos(ang);
    pend(i,2) = pstart(i,2) + f_width*sin(ang);
end
%%
plot( [x pstart(1,1)] , [y pstart(1,2)] , 'LineWidth', 1, 'Color', color_pres);
for i=1:num_bar
    plot( [pstart(i,1) pend(i,1)] , [pstart(i,2) pend(i,2)] , 'LineWidth', 1, 'Color', color_pres);
end