2012年全國(guó)數(shù)學(xué)建模競(jìng)賽優(yōu)秀選劉楊郭

1、 填寫）參賽（由第九屆“杯”數(shù)學(xué)建模競(jìng)賽風(fēng)計(jì)算的初步探究題 目摘要：本文針對(duì)風(fēng)的度量模型和算法進(jìn)行了深入研究，并在常規(guī)反演方法的基礎(chǔ)上提出了改進(jìn)的措施，根據(jù)題中給的數(shù)據(jù)定量計(jì)算出與常規(guī)方法比較，改進(jìn)后的效果良好。風(fēng)風(fēng)矢場(chǎng)信息，問題一：分別建立了以地球球心和同步為坐標(biāo)原點(diǎn)的兩個(gè)三維坐標(biāo)系，并計(jì)算出兩個(gè)坐標(biāo)系的變換關(guān)系矩陣 T，建立了與灰度矩陣之間的換算關(guān)系。求解出灰度矩陣中第 500 行和第 500、501、502 列三個(gè)交叉點(diǎn)處對(duì)應(yīng)的分別為：(46.3775o，32.9053o) 、(46.4623o，32.8998o)、(46.5476o，32.8939o)。矩陣形式的坐標(biāo)文件 jwd.tx

2、t 及云圖 spic.jpg 存于附件中。問題二：通過分析和梯度交叉相關(guān)系數(shù)的云圖的灰度特性和紋理特征，提出了基于最大灰度風(fēng)反演方法，利用灰度閾值云檢測(cè)方法去除地物對(duì)反演結(jié)果的影響，并且通過計(jì)算前后兩個(gè)時(shí)次的風(fēng)矢量，根據(jù)時(shí)間連續(xù)性和空間連續(xù)性準(zhǔn)則進(jìn)行質(zhì)量控制。計(jì)算得到的 21:00 時(shí)刻各附件文件 winsh.txt，五個(gè)規(guī)定點(diǎn)處的風(fēng)矢如表 1。處的二維風(fēng)矢見問題三：依據(jù)最大圖像信息熵準(zhǔn)則，自適應(yīng)調(diào)整匹配窗口大小，并提出多尺度搜索策略，有效減少計(jì)算量。采用該策略得到五個(gè)規(guī)定點(diǎn)的風(fēng)矢見表 3。問題四：利線性插值方法，計(jì)算非零風(fēng)矢量在云圖中的亮溫值和該點(diǎn)處的溫度廓線，通過亮溫在溫度廓線中的相應(yīng)位置

3、查找非零風(fēng)矢量的高度。五個(gè)規(guī)定點(diǎn)風(fēng)矢量所在的等壓面見表 4。關(guān)鍵字：風(fēng)；坐標(biāo)變換；最大交叉相關(guān)系數(shù)；信息熵；雙線性插值- 1 -一、問題重述云圖在掌握大氣環(huán)流、中長(zhǎng)期天氣預(yù)報(bào)以及性天氣學(xué)的研究中有重要作用。它由地球同步上的紅外探測(cè)儀探測(cè)地球上空的溫度數(shù)據(jù)再轉(zhuǎn)換成灰度數(shù)據(jù)制作而成。紅外探測(cè)儀掃描采樣時(shí)，按步進(jìn)角(南北方向)和行掃描角（東西方向）均為 140 微?。? 弧度1000000 微?。┎蓸?。在與地球中心的連線和地球表面的交點(diǎn)（稱為星下點(diǎn)）處的分辨率大約是 5 公里。本題提供的探測(cè)數(shù)據(jù)文件都是 22882288 的灰度值矩陣，矩陣的每個(gè)元素都對(duì)應(yīng)地球上或地球外的一個(gè)探測(cè)點(diǎn)（或稱采樣點(diǎn)）。同

4、步離地球中心的高度為42164000 米, 星下點(diǎn)在東經(jīng) 86.5 度, 北緯 0 度，星下點(diǎn)對(duì)應(yīng)的矩陣元素位于矩陣的第 1145 行和第 1145 列相交處。為解答本題，首先要確定灰度矩陣中每個(gè)元素對(duì)應(yīng)的采樣點(diǎn)在地球上的經(jīng)緯度。地球可視為理想橢球，這個(gè)理想橢球可以由地球的一個(gè)經(jīng)過南北極的橢圓截面繞南北極的連線旋轉(zhuǎn)而得到，橢圓截面的長(zhǎng)半軸（赤道半徑）=6378136.5m,短半軸（極半徑）=6356751.8m；可將灰度矩陣中非負(fù)元素的行列號(hào)由上北下南、左的地圖規(guī)則換算成地球上觀測(cè)大氣環(huán)流情況的一個(gè)方法是在坐標(biāo)。云圖上標(biāo)出風(fēng)矢。風(fēng)矢的大小和方向由云塊移動(dòng)的速度決定。風(fēng)矢與風(fēng)的速度有所不同，如某

5、個(gè)臺(tái)風(fēng)中一些區(qū)域的風(fēng)速可達(dá)每秒五、六十米，而臺(tái)風(fēng)（看作云塊）中心的移動(dòng)速度可能僅每小時(shí)十多公里。沒有云或云塊不穩(wěn)定處的風(fēng)矢規(guī)定為零風(fēng)矢，這種用云塊的移動(dòng)所定義的風(fēng)矢被稱為風(fēng)。氣象部門已經(jīng)有一些方法根據(jù)變化的云圖計(jì)算風(fēng)，這類方法稱為風(fēng)方法。計(jì)算風(fēng)時(shí)通常將云塊大小限定為1616 個(gè)像素，搜索范圍限定為 6464 個(gè)像素。本題的主要目的是希望大家充分利用準(zhǔn)確地描述實(shí)際風(fēng)矢場(chǎng)的度量模型和算法。具體要求解決如下問題：圖像數(shù)據(jù)及其特點(diǎn)建立盡可能1. 換算視場(chǎng)坐標(biāo)。給出灰度矩陣元素行列號(hào)對(duì)應(yīng)于坐標(biāo)的換算公式，建立矩陣形式的坐標(biāo)文本文件，這里矩陣的第 i 行與第 j 列，分別對(duì)應(yīng)灰度矩陣的 450+i 行與

6、450+j 列，矩陣元素是（經(jīng)度，緯度）這種形式的二維數(shù)組，給出結(jié)果的范圍為：灰度矩陣中的第 451 行至 550 行，第 451 列至 550列，文本文件取名為 jwd.txt；而在中給出第 500 行和第 500、501、502 列三個(gè)交叉點(diǎn)處對(duì)應(yīng)的。進(jìn)而將灰度矩陣 IR1_2100.mat 轉(zhuǎn)化為云圖，并在這個(gè)云圖上添加海岸線。2. 根 據(jù) 三 個(gè) 相 鄰 時(shí) 刻 的數(shù) 據(jù) （ IR1_2030.mat, IR1_2100.mat,IR1_2130.mat），用適當(dāng)?shù)姆椒ń?21：00 時(shí)刻風(fēng)矢場(chǎng)的度量模型，給出計(jì)算方案，并利用該方法計(jì)算 21：00 時(shí)刻在各個(gè)（整數(shù)）處的二維風(fēng)矢。確

7、定一個(gè)二維風(fēng)矢需要 4 個(gè)指標(biāo)：起始點(diǎn)的緯度、經(jīng)度，風(fēng)矢的方向（：從正北方順時(shí)針方向量得的角度）和大?。ǘ戎帘本?40 度，東經(jīng) 46 度至 126 度。在：米/秒）。計(jì)算范圍為：南緯 40中給出上面計(jì)算結(jié)果中全部非零風(fēng)矢的個(gè)數(shù)，并給出在緯度 26 度，經(jīng)度分別是 52，53，54，57，58（度）處的風(fēng)矢，這些風(fēng)矢的表示形式可參照下面問題 4 中對(duì) winsh.txt 要求的格式。給出對(duì)上述窗口大小和搜索范圍的限定的理解或評(píng)價(jià)。3. 為了提高像素塊匹配質(zhì)量，取消問題 2 中窗口大小和搜索范圍的限定，利用圖像特征匹配等，設(shè)計(jì)能自適應(yīng)地確定窗口大小和搜索范圍的有效方- 2 -法。進(jìn)而建立風(fēng)矢場(chǎng)的

8、度量模型和算法。并利用此模型和算法計(jì)算問題 2 中在緯度 26 度，經(jīng)度分別是 52，53，54，57，58（度）處的風(fēng)矢。定量分析分析比較建立的這兩種風(fēng)矢計(jì)算方法的優(yōu)劣，結(jié)論要有定量數(shù)據(jù)支撐。并提出關(guān)于云圖數(shù)據(jù)量和風(fēng)風(fēng)矢場(chǎng)計(jì)算研究的建議。4. 應(yīng)用附件中全球溫度數(shù)值預(yù)報(bào)文件 temp3.mat 計(jì)算問題 2 中非零風(fēng)矢所在的等壓面（用該時(shí)刻該處大氣的壓強(qiáng)值標(biāo)志，物理為毫巴）。并在中給出在緯度 26 度，經(jīng)度分別是 52，53，54，57，58（度）處附近求出的風(fēng)矢量所在的等壓面。用文本文件列表的形式給出問題 2 和本問題的計(jì)算結(jié)果，保存為電子文檔 winsh.txt. 最后在問題 1 中畫出

9、的云圖上添加所計(jì)算出的非零風(fēng)矢量結(jié)果（參照附錄中的風(fēng)圖）, 并將此標(biāo)有風(fēng)矢圖標(biāo)的云圖（即風(fēng)圖）保存為 PDF 文件或 JPEG 文件，取名為 winsp.pdf 或 winsp.jpg。二、模型假設(shè)假設(shè)地球?yàn)槔硐霗E球,可以由地球的一個(gè)經(jīng)過南北極的橢圓截面繞南北極的連線旋轉(zhuǎn)而得到，橢圓截面的長(zhǎng)半軸（赤道半徑）=6378136.5m，短半軸（極半徑）=6356751.8m。三、符號(hào)說明 p ：表示 P 點(diǎn)經(jīng)度角； p ：P 點(diǎn)緯度角；So ：星下點(diǎn)；O ：赤道上的 0 度經(jīng)度點(diǎn)； p ： p ：從星下點(diǎn)起水平掃描角；從赤道起豎直抬升角。四、模型的建立與求解問題一的解答問題分析問題一首先要給出灰度矩

10、陣中各元素行列號(hào)對(duì)應(yīng)于式，也就是根據(jù)灰度矩陣中的行列號(hào)值，來求解對(duì)應(yīng)的坐標(biāo)的換算公問題?；叶染仃囀峭缴霞t外探測(cè)儀掃描時(shí)的探測(cè)點(diǎn)，根據(jù)紅外探測(cè)儀的掃描原理可知，灰度矩陣中的元素與位置存在固定關(guān)系，同時(shí)位心和探測(cè)點(diǎn)的之間也存在固定的關(guān)系。因此，可以先通過計(jì)算它們之間逐步變換關(guān)系，最終確定灰度矩陣和之間的變換關(guān)系。問題一還需要將灰度矩陣轉(zhuǎn)化為云圖，并在云圖上添加海岸線。只要確定了坐標(biāo)與灰度矩陣之間的變換關(guān)系，那么就可以將海岸線的經(jīng)緯度轉(zhuǎn)換成灰度矩陣之后在云圖上表示出來。4.1.2 坐標(biāo)轉(zhuǎn)換模型的建立以地球球心為坐標(biāo)原點(diǎn)，赤道平面為 XOY 平面，地球球心指向同步的方向?yàn)?X 軸，X 軸沿著赤道平面

11、順時(shí)針旋轉(zhuǎn) 90o 為 Y 軸，地球球心指向北極方向?yàn)?Z 軸，建立球心坐標(biāo)系 O-XYZ。所在位置 S 點(diǎn)為坐標(biāo)原點(diǎn)，同步指向地球球心為 X 以地球同步- 3 -軸，X 軸沿著赤道平面的拓展平面順時(shí)針旋轉(zhuǎn) 90o 為Y 軸，北極方向?yàn)閆 軸建立三維直角坐標(biāo)系S-X Y Z 。球心坐標(biāo)系 O-XYZ 與坐標(biāo)系S-X Y Z 的位置關(guān)系如圖 1 所示：ZZYPphpCOps pSoX SX0pDO Y圖 1 坐標(biāo)關(guān)系及相關(guān)要素圖 1 中，P 點(diǎn)表示地球表面的某一探測(cè)點(diǎn)（灰度矩陣中第 i 行，第 j 列的元素），D 點(diǎn)是 P 在赤道平面上的投影，C 點(diǎn)是過 D 點(diǎn)向 X 軸作垂線后與 X 軸的交點(diǎn)

12、，O 點(diǎn)是赤道與零度經(jīng)度線的焦點(diǎn)，So 是星下點(diǎn)，h 表示地心與同步之間的距離 OS 42164000m 。從而可知： p POD 是 P 點(diǎn)的緯度角， p O OD是 P 點(diǎn)的經(jīng)度角。下面求解兩個(gè)坐標(biāo)系之間的變換關(guān)系：由 O-XYZ 和 p O OD 的位置關(guān)系可知，S-X Y Z 坐標(biāo)系是由 O-XYZ 坐標(biāo)系先向右平移 h，然后再繞 Z 軸方向旋轉(zhuǎn) 180o 得到。由三維圖形的幾何變換原理1：平移變換：10010Ty001Tz00100100001000 xi yi zi 1 xizi 1 xizi 1 ,yiyi（1.1）000h01Tx1繞 Z 軸旋轉(zhuǎn)變換：- 4 - coszsin

13、 z001000 sin cos0 x 1 x10zz y zy z1iiiii i00001（1.2）01 000010 00 xizi 10 yi 0 01 .假設(shè)S-X Y Z 坐標(biāo)系與 O-XYZ 坐標(biāo)系的變換關(guān)系為：xi yi zi 1 xizi 1T ,yi其中 T 為兩個(gè)坐標(biāo)系的變換矩陣，由式（1.1）、（1.2）可知：0 1011001000001 00001001 000010000 00 0T ,0h1 0 00 1 0 h0 1 1 00 所以兩個(gè)坐標(biāo)的變換關(guān)系為：101 0000100 00 xi yi zi 1 xizi 1 .yi（1.3） 0 h0 1 （ p ,

14、 p ），可以轉(zhuǎn)化為求解 P 點(diǎn)在 O-XYZ 坐標(biāo)題中需要求解的 P 點(diǎn)系下的坐標(biāo) P(xp , yp , zp ) ；根據(jù)坐標(biāo)變換關(guān)系又可以轉(zhuǎn)化為求 P 點(diǎn)在S-X Y Z 坐標(biāo)系下的坐標(biāo) P(xp , yp , zp ) ，而 P(xp , yp , zp ) 與 P 點(diǎn)在灰度矩陣中的具置P(i, j) 之間存在特定關(guān)系，根據(jù)它們之間的特定關(guān)系就可以建立方程求解出（ p , p ）情況P( xp , yp , zp )。這個(gè)過程的正向推導(dǎo)就可以在已知 P 點(diǎn)的下求解 P 點(diǎn)在灰度矩陣中的具置 P(i, j) ，而反向推導(dǎo)就可以在已知 P 點(diǎn)在灰（ p , p ）。度矩陣中的具置 P(i,

15、 j) 情況下求解 P 點(diǎn)的問題求解由灰度矩陣推導(dǎo)對(duì)于問題一，假設(shè) P 點(diǎn)的當(dāng)前位置在灰度矩陣中第 i 行和第 j 列，其中 451 i 550， 451 j 550 。根據(jù)題中紅外探測(cè)儀掃描規(guī)則和已知條件可知，P點(diǎn)可由同步 從星下點(diǎn)先水平向西旋轉(zhuǎn) p 角度，再豎直向北旋轉(zhuǎn) p 得到。其中： p ( j 1145)0 , p (1145 i)0 ,（1.4）（1.5）- 5 -0 、0 分別表示紅外探測(cè)儀行掃描角和步進(jìn)角，其值均為 140 微弧。如圖 1，假設(shè) PP(xp , yp , zp ) 。在直角三角形點(diǎn)在坐標(biāo)系S-X Y Z 中坐標(biāo)為 P(xp , yp , zp ) ，下面求解和P

16、DS 中，：DSPSsin CD y x 2 y,2（1.6）ppppsin DP z x 2 y 2 z.2（1.7）ppppp另外，題中假設(shè)地球可視為理想橢球，這個(gè)理想橢球可以由地球的一個(gè)經(jīng)過南北極的橢圓截面繞南北極的連線旋轉(zhuǎn)而得到，橢圓截面的長(zhǎng)半軸（赤道半徑）a=6378136.5m, 短半軸（極半徑）b=6356751.8m，標(biāo)系下的方程為：到地球在 O-XYZ 坐x2 y2z2 1 .a2b2通過坐標(biāo)轉(zhuǎn)換，將地球轉(zhuǎn)換到S-X Y Z 坐標(biāo)系下，得到地球在S-X Y Z 坐標(biāo)系下的方程為：(x h)2 y 2 z 2 1 ,（1.8）a2b2因?yàn)?P(xp , yp , zp ) 在地

17、球上，所以也滿足（1.8）式。聯(lián)立式子（1.4）、（1.5）、（1.6）、（1.7）、（1.8），當(dāng)給定 i 和 j 的值時(shí)，就可以計(jì)算出 xp 、 yp 和 zp 的值。在聯(lián)立求解時(shí)，根據(jù)（1.6）和（1.7）式將 yp 和zp 用 xp 表示，然后將 yp 和 zp 帶入（1.8）式中得到一個(gè)關(guān)于 xp 的一個(gè)復(fù)雜的一元二次方程， xp 值的范圍為(h SSo, h) (h a, h) 。這里采用黃金分割法求解這個(gè)復(fù)雜的一元二次方程。即首先將 xp 的一元二次方程轉(zhuǎn)化為 f (xp ) 0f (x ) 在(h a, h) 內(nèi)必然有一個(gè)零點(diǎn)， 此點(diǎn)即為 xp 的值。取的形式，那么f ( x)

18、 f ( x ),分別從h 和 h a 兩個(gè)端點(diǎn)逐步向內(nèi)采用黃金分割法求解 f (x) 的零點(diǎn)，當(dāng) f (x) 的值逐步近為 0，滿足給定的一個(gè)極小的誤差值時(shí)，此時(shí)的 x 點(diǎn)即可認(rèn)為是 xp 點(diǎn)的值。再根據(jù)前面求解的式（1.3）坐標(biāo)變換關(guān)系：xi yi zi 1 xizi 1T ,yixp yp zp 1 xp yp zp 1T.1則從而就確定了 P 點(diǎn)在 O-XYZ 坐標(biāo)系下的坐標(biāo) P(xp , yp , zp ) 。角 p 和 p 與 P 點(diǎn)的坐標(biāo) P(xp , yp , zp ) 有如下關(guān)根據(jù)圖 1 可知，P 點(diǎn)的系：- 6 -ypsin( ) ,sopx2 y2ppzpsin .px2

19、 y2 z2ppp其中，so 86.5 表示o星下點(diǎn)所在位置的經(jīng)度。進(jìn)而：yp p so arcsin( )x2 y2pp.z arcsin ppx2 y2 z2ppp至此，已知灰度矩陣中元素的行列號(hào)求導(dǎo)出了該元素的坐標(biāo)。根據(jù)前面的推導(dǎo)，求得第 500 行和第 500、501、502 列三個(gè)交叉點(diǎn)處對(duì)應(yīng)分別為：(46.3775o，32.9053o) 、(46.4623o，32.8998o)、(46.5476o，32.8939o)。的(2)由要在推導(dǎo)灰度矩陣云圖上添加海岸線，就必須要把給定的海岸線坐標(biāo)轉(zhuǎn)化到灰度矩陣中，求解出在灰度矩陣的具置。由 似。推導(dǎo)灰度矩陣是由灰度矩陣推導(dǎo)的逆過程，其推導(dǎo)原

20、理類( p ,p )，可求 P 點(diǎn)在 O-XYZ 坐標(biāo)系下的坐標(biāo)已知 P 點(diǎn)的P( xp , yp , zp )：xOP cos p cos(so p )p yOP cos sin( ) .（1.9）ppsopzp OP sin p又因?yàn)檫^ P 點(diǎn)和地球南北軸（Z 軸）形成的平面與地球球面的交線是過球心的一個(gè)橢圓，所以 OP(a cos )2 (bsin )2 ，即求解出 P( x, y , z )。ppppp然后再根據(jù)坐標(biāo)變換關(guān)系矩陣（1.3）式，可以求解出 P 點(diǎn)在S-X Y Z 坐標(biāo)系下的坐標(biāo) P(xp , yp , zp ) 。由（1.6）式和（1.7）式： p arcsin ( yp

21、 p arcsin( zp 最后由（1.4）式和（1.5）式xp y ) ,22pxp y z ) .222ppP 點(diǎn)在灰度矩陣中的行列值 P(i, j) ：i 1145 p0 . j 1145 p0- 7 -式中： p 0 和 p 0 分別表示對(duì) p 0 和 p 0 的值取整數(shù)。下面給出添加了海岸線的由如圖 2 所示?；叶染仃嘔R1_2100.mat轉(zhuǎn)化的云圖，圖 2 加入海岸線信息的云圖問題二的解答問題分析風(fēng)通常追蹤圖像塊(示蹤云)在三幅連續(xù)云圖中的移動(dòng)來估算風(fēng)速,并由示蹤云的高度指定風(fēng)場(chǎng)的高度.其原理如圖 3 所示，左側(cè)為 t0 時(shí)刻的云圖，A 區(qū)域?yàn)槭聚櫾茐K，右側(cè)為 t1 時(shí)刻云圖，在

22、搜索區(qū)域 C 中尋找示蹤云塊在 t1 時(shí)刻的位置。實(shí)際處理中主要包括四個(gè)步驟：1）數(shù)據(jù)預(yù)處理：對(duì)用于計(jì)算風(fēng)的云圖進(jìn)行精確定位；示蹤云追蹤：在目標(biāo)圖像中搜索示蹤云的位置；高度指定：指定計(jì)算風(fēng)場(chǎng)的高度；質(zhì)量控制：利用三幅云圖計(jì)算出的兩個(gè)風(fēng)場(chǎng)結(jié)果進(jìn)行連續(xù)性和空間一致- 8 -性檢驗(yàn)等質(zhì)量控制。圖 3風(fēng)計(jì)算原理示意圖4.2.2 基于灰度和梯度最大交叉相關(guān)系數(shù)的風(fēng)反演算法um示蹤云追蹤方法中目前使用最廣泛的是最大交叉相關(guān)系數(shù)（ cross-correlation coefficient, MCC）法2，即在搜索區(qū)域滑動(dòng)，計(jì)算示蹤云塊（為 1616 像素圖像塊）在目標(biāo)搜索區(qū)（為 6464 像素圖塊）里的交

23、叉相關(guān)系數(shù)的最大區(qū)域即為目標(biāo)區(qū)域，交叉相關(guān)系數(shù) R1(m, n) 為： f (i, j) f g(i m, j n) gijR (m, n) .（2.1）11/ 2 f (i, j) f 2 g(i m, j n) g 2ijij其中i, j 1,16 分別為示蹤云模板中的行與列號(hào)， m, n 1, 48為目標(biāo)搜索區(qū)位移改變量， f (i, j) 為示蹤云模板的灰度值， g(i, j) 為目標(biāo)搜索區(qū)云的灰度值，記 1 f (i, j ) 為示蹤云的平均灰度，g 1 g(i m, j n) 為目標(biāo)搜索fN 2N 2ijij區(qū)的平均灰度，其中取 N 16 。傳統(tǒng)的 MCC 方法僅利用了圖像的灰度信

24、息，而在云圖上，相鄰區(qū)域的灰度具有空間連續(xù)性，在較小的圖像塊中灰度基本相同，很難分辨出云塊的紋理結(jié)構(gòu)特征。圖像的灰度梯度反映了相鄰區(qū)域的灰度變化，能夠提取圖像中目標(biāo)物的輪廓信息，如圖 4 所示。因此，本文將圖像灰度的梯度信息加入相關(guān)系數(shù)的計(jì)算中，梯度信息的交叉相關(guān)系數(shù)為： f (i, j) fijg(i m, j n) gR (m, n) .（2.2）21/ 22 2f (i, j) fg(i m, j n) gijij其中f (i, j) 為 f (i, j)的灰度梯度， g(i, j) 為 g(i, j) 的灰度梯度， f 和 g 為平均灰度梯度。構(gòu)造函數(shù)如下：Jm, n (1 )R1(m

25、, n) R2 (m, n) ,（2.3）其中 為權(quán)重系數(shù)（本文取 0.5 ）。當(dāng) J 取得最大時(shí)對(duì)應(yīng)的區(qū)域即為示蹤云下一時(shí)刻的位置，目標(biāo)搜索區(qū)域內(nèi)相對(duì)于示蹤云在 x, y 方向上分別移動(dòng)了m 和 n 個(gè)像素點(diǎn)，建立地球模型可計(jì)算出目標(biāo)云相對(duì)于示蹤云的位移量，除以兩幅圖- 9 -像的觀測(cè)時(shí)間間隔到風(fēng)矢量。(a)灰度(b)梯度圖 4 灰度云圖及其梯度對(duì)比4.2.3 二維風(fēng)矢量計(jì)算方法假設(shè)需要求解起始點(diǎn) A 點(diǎn)風(fēng)矢的方向和大小，現(xiàn)已知 A 點(diǎn)的坐標(biāo)(A,A) ，且通過相關(guān)的計(jì)算得出了下一時(shí)刻點(diǎn) B 點(diǎn)的坐標(biāo)(B ,B ) ；那么 A 點(diǎn)風(fēng)矢的方向是從正北方向順時(shí)針旋轉(zhuǎn)到風(fēng)的來向之間的角度（正北風(fēng)的

26、風(fēng)矢方向?yàn)?0o，正南風(fēng)的風(fēng)矢方向?yàn)?180o，正東風(fēng)的風(fēng)矢方向?yàn)?90o，正西風(fēng)的風(fēng)矢方向?yàn)?270o）,A 點(diǎn)到 B 點(diǎn)的平均風(fēng)速為 A 點(diǎn)風(fēng)矢的大小。A 點(diǎn)風(fēng)矢大小的計(jì)算：vA .tABAB表示 A 點(diǎn)到 B 點(diǎn)的弧長(zhǎng)，在半個(gè)小時(shí)內(nèi)的長(zhǎng)度相對(duì)較小，可以用AB近似代替，所以：vA .tAB OB OA (x x )2 ( y y )2 (z z.)2又因?yàn)椋築ABABA再根據(jù)前面求解的（1.9）式，并將 A 和 B 兩點(diǎn)的出vA 。A 點(diǎn)風(fēng)矢方向的計(jì)算：以經(jīng)度為橫坐標(biāo)，緯度為縱坐標(biāo)建立二維坐標(biāo)系如圖 5 所示：值代入，可以求解- 10 -ABABABOAB圖 5坐標(biāo)因?yàn)轱L(fēng)向是指正北方向順時(shí)

27、針到風(fēng)的來向之間的角度，風(fēng)向的范圍是0o ,360o )，根據(jù) A 點(diǎn)和 B 點(diǎn)經(jīng)度的大小分成兩種情況計(jì)算風(fēng)向。當(dāng) 時(shí)，其風(fēng)向的范圍是0o ,180o ，如圖 5 第一象限中AB 所示，此時(shí)AB的風(fēng)向 為： arccos(A B 0o ,180o ;),( )2 ()2ABAB當(dāng) 時(shí)，其風(fēng)向的范圍是(180o ,360o ) ，如圖 5 第四象限中AB 所示，此AB時(shí)的風(fēng)向 為： 2 arccos(A B (180o , 360o ) .),( )2 ()2ABAB4.2.4 云檢測(cè)與質(zhì)量控制在進(jìn)行風(fēng)反演之前，首先要判斷該區(qū)域是否存在云團(tuán)，無云區(qū)在相鄰時(shí)刻的云圖中是保持的，因此無法反演風(fēng)矢量。

28、在云圖中，灰度比較暗的區(qū)域一般為陸地和海洋，而區(qū)域?yàn)樵葡?。由于只有紅外通道的云圖數(shù)據(jù)，因此本文采用閾值法進(jìn)行云檢測(cè)，如果示蹤云的平均灰度 f 大于某一閾值（本文將灰度閾值設(shè)為 100）則有云區(qū)域。風(fēng)的質(zhì)量控制根據(jù)天氣學(xué)中大氣運(yùn)動(dòng)的原理進(jìn)行如下控制3：時(shí)間連續(xù)性質(zhì)量控制，對(duì)前后 2 時(shí)次計(jì)算出的風(fēng)矢量差大于二者風(fēng)速均值的風(fēng)剔除；空間連續(xù)性質(zhì)量控制，對(duì)每一個(gè)風(fēng)矢量檢驗(yàn)與周圍風(fēng)是否有連續(xù)性。4.2.5 試驗(yàn)結(jié)果與分析本文根據(jù)三個(gè)相鄰時(shí)刻的 IR1_2130.mat)，利用上述方法進(jìn)行數(shù)據(jù)(IR1_2030.mat ， IR1_2100.mat ，風(fēng)反演。圖 6(a)、(b)、(c)分別為三個(gè)衛(wèi)星數(shù)據(jù)

29、經(jīng)過云檢測(cè)后的子圖，大致范圍為：北緯 23-29,東經(jīng) 52-58。通過連續(xù)播放三幅子圖可以看出，云圖中的云系整體向東移動(dòng)。- 11 -(a) IR1_2030(b) IR1_2100(c) IR1_2130圖 6 相鄰時(shí)刻云圖子圖圖 7 是圖 6 的風(fēng)結(jié)果示意圖。表 1 給出了部分整數(shù)經(jīng)緯坐標(biāo)位置風(fēng)反演的風(fēng)速和風(fēng)向。從結(jié)果可以看出， 風(fēng)的風(fēng)向和圖 6 中云團(tuán)的移動(dòng)方向基本一致，風(fēng)速大小也較為合理。但是，部分風(fēng)矢量與周圍風(fēng)的空間連續(xù)性較差，而且在部分位置沒有計(jì)算出風(fēng)矢量。根據(jù)實(shí)驗(yàn)發(fā)現(xiàn)，造成上訴現(xiàn)象的主要原因是由于像素塊匹配窗口過小引起。圖 8 是某一 1616 匹配窗口放大后的圖像，由于窗口過

30、小，且灰度分布均勻，很難從中提取特征，因此容易造成云導(dǎo)風(fēng)結(jié)果的不穩(wěn)定。圖 7風(fēng)效果圖圖 81616 匹配窗口表 1 問題二結(jié)果注：表中“-”表示數(shù)據(jù)無法提取4.3 問題三的解答從問題二的結(jié)果可知，當(dāng)匹配窗口過小時(shí)，由于缺乏云圖紋理信息，很容- 12 -經(jīng)度緯度角度（北順）速度（米/秒5226-532626420.4542627124.4572627026.25826-非零風(fēng)矢?jìng)€(gè)數(shù)1302易造成誤匹配。將窗口大小設(shè)為 3232，搜索范圍限定為 9696 時(shí)，反演結(jié)果如圖 9 和表 2 所示。增大搜索窗口后的非矢數(shù)量，以及風(fēng)矢之間的空間連續(xù)性較問題二結(jié)果相比都有了明顯改善。但是，窗口變大后，也增加

31、了計(jì)算的時(shí)間復(fù)雜度。因此，本文提出了基于圖像熵的窗口自適應(yīng)方法和多尺度搜索策略。圖 9 固定窗口 3232風(fēng)效果圖表 2 固定窗口 323 2 結(jié)果圖像的熵是一種特征的統(tǒng)計(jì)形式，它反映了圖像中平均信息量的多少，熵越大，圖像信息量越大。圖像熵的定義如下：Entroy p(i) log2 p(i) ,i其中 i 為灰度階，i=1, 2, , 255， p( i )表示灰階 i 在圖像中出現(xiàn)的頻率。窗口自適應(yīng)算法步驟如下：以 1616 為最小值，3232 為最大值，分別計(jì)算介于兩者之間所有匹配窗口的熵；比較各窗口的熵值，熵最大的窗口作為匹配窗口； (3)匹配窗口的長(zhǎng)寬各拓展 32 像素作為搜索窗口。

32、在搜索窗口較大時(shí)，可增大搜索步長(zhǎng) step (本文step 3 )進(jìn)行粗略的搜索，得到初步匹配結(jié)果。然后，在上下左右step 個(gè)像素點(diǎn)范圍內(nèi)搜索，得到最終匹配結(jié)果。圖 10 和表 3 給出了自適應(yīng)窗口方法的部分計(jì)算結(jié)果，圖 11 是得到的風(fēng)圖。通過與問題二的結(jié)果比較，自適應(yīng)窗口方法具有更強(qiáng)的空間連續(xù)性，很好反映了云團(tuán)的運(yùn)動(dòng)趨勢(shì)。- 13 -經(jīng)度緯度角度（北順）速度（米/秒）52262548.8532627225.2542627825.8572627521.5582627317.3圖 10 自適應(yīng)窗口風(fēng)示意圖表 3 自適應(yīng)窗口 3636 結(jié)果建議：通過頻域變換或者小波方法提取圖像中的高頻成分，根

33、據(jù)紋理特征進(jìn)行匹配；研究非線性的搜索策略，通過分析在各方向上匹配相關(guān)系數(shù)的變化，選擇下一個(gè)搜索方向，從而降低計(jì)算復(fù)雜度。- 14 -經(jīng)度緯度角度（北順）速度（米/秒）522627321.2532628021.0542627825.8572627521.56582627317.3非零風(fēng)矢?jìng)€(gè)數(shù)1921圖 11風(fēng)圖問題四的解答問題分析在問題二的解答過程中，已經(jīng)建立了 21:00 時(shí)刻風(fēng)矢場(chǎng)的度量模型，并計(jì)算出了該時(shí)刻在各個(gè)（整數(shù)）處的全部非零二維風(fēng)矢。本問題要求運(yùn)用全球溫度數(shù)值預(yù)報(bào)產(chǎn)品，計(jì)算上述非零風(fēng)矢所在的等壓面。由于地球各緯度上空各處的溫度值與相應(yīng)點(diǎn)大氣的氣壓值有著一定的溫度氣壓廓線關(guān)系，因此要

34、計(jì)算非零風(fēng)矢所在的等壓面，只需要計(jì)算得到各處非零風(fēng)矢所在位置的溫度值即可。4.4.2 問題求解通過對(duì)該問題的分析，制定了求解該問題的基本思路：由追蹤云團(tuán)的移動(dòng)計(jì)算得到每一個(gè)風(fēng)矢量（問題二中已完成），再由云團(tuán)代表的風(fēng)矢溫度，通過一定的轉(zhuǎn)換獲取風(fēng)場(chǎng)的高度。按照上述思路，設(shè)定解答該問題的步驟如下：第一步：導(dǎo)出問題二中求得的 21:00 時(shí)刻在各個(gè)整數(shù)處的二維非零- 15 -風(fēng)矢；通過問題一中給出的坐標(biāo)變換關(guān)系，求出各個(gè)非零風(fēng)矢在云圖上對(duì)應(yīng)的精確位置（用坐標(biāo)表示）。第二步：結(jié)合題中所給灰度數(shù)據(jù)與溫度數(shù)據(jù)的轉(zhuǎn)換關(guān)系（k_temp.txt），采線性插值法，求出每個(gè)精確位置處的溫度值 tb。第三步：對(duì)于同一非零風(fēng)矢的精確位置，同樣采線性插值法，求出其在全球溫度數(shù)值預(yù)報(bào)文件 temp3.mat 中對(duì)應(yīng)的各層的溫度值 t（共 36 個(gè)）。第四步：將 tb 與 36 個(gè) t 做對(duì)比，取與 tb 最接近的兩個(gè)值，采用線性插值法，結(jié)合文件 temp3.mat，求出 t