2012年全國數(shù)學建模競賽優(yōu)秀選劉楊郭

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

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

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

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

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

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

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

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

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

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

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

12、，O 點是赤道與零度經(jīng)度線的焦點，So 是星下點，h 表示地心與同步之間的距離 OS 42164000m 。從而可知： p POD 是 P 點的緯度角， p O OD是 P 點的經(jīng)度角。下面求解兩個坐標系之間的變換關(guān)系：由 O-XYZ 和 p O OD 的位置關(guān)系可知，S-X Y Z 坐標系是由 O-XYZ 坐標系先向右平移 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 坐標系與 O-XYZ 坐標系的變換關(guān)系為：xi yi zi 1 xizi 1T ,yi其中 T 為兩個坐標系的變換矩陣，由式（1.1）、（1.2）可知：0 1011001000001 00001001 000010000 00 0T ,0h1 0 00 1 0 h0 1 1 00 所以兩個坐標的變換關(guān)系為：101 0000100 00 xi yi zi 1 xizi 1 .yi（1.3） 0 h0 1 （ p ,

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

15、 j) 情況下求解 P 點的問題求解由灰度矩陣推導對于問題一，假設(shè) P 點的當前位置在灰度矩陣中第 i 行和第 j 列，其中 451 i 550， 451 j 550 。根據(jù)題中紅外探測儀掃描規(guī)則和已知條件可知，P點可由同步 從星下點先水平向西旋轉(zhuǎn) p 角度，再豎直向北旋轉(zhuǎn) p 得到。其中： p ( j 1145)0 , p (1145 i)0 ,（1.4）（1.5）- 5 -0 、0 分別表示紅外探測儀行掃描角和步進角，其值均為 140 微弧。如圖 1，假設(shè) PP(xp , yp , zp ) 。在直角三角形點在坐標系S-X Y Z 中坐標為 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è)地球可視為理想橢球，這個理想橢球可以由地球的一個經(jīng)過南北極的橢圓截面繞南北極的連線旋轉(zhuǎn)而得到，橢圓截面的長半軸（赤道半徑）a=6378136.5m, 短半軸（極半徑）b=6356751.8m，標系下的方程為：到地球在 O-XYZ 坐x2 y2z2 1 .a2b2通過坐標轉(zhuǎn)換，將地球轉(zhuǎn)換到S-X Y Z 坐標系下，得到地球在S-X Y Z 坐標系下的方程為：(x h)2 y 2 z 2 1 ,（1.8）a2b2因為 P(xp , yp , zp ) 在地

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

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

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

20、理類( p ,p )，可求 P 點在 O-XYZ 坐標系下的坐標已知 P 點的P( xp , yp , zp )：xOP cos p cos(so p )p yOP cos sin( ) .（1.9）ppsopzp OP sin p又因為過 P 點和地球南北軸（Z 軸）形成的平面與地球球面的交線是過球心的一個橢圓，所以 OP(a cos )2 (bsin )2 ，即求解出 P( x, y , z )。ppppp然后再根據(jù)坐標變換關(guān)系矩陣（1.3）式，可以求解出 P 點在S-X Y Z 坐標系下的坐標 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 點在灰度矩陣中的行列值 P(i, j) ：i 1145 p0 . j 1145 p0- 7 -式中： p 0 和 p 0 分別表示對 p 0 和 p 0 的值取整數(shù)。下面給出添加了海岸線的由如圖 2 所示?；叶染仃嘔R1_2100.mat轉(zhuǎn)化的云圖，圖 2 加入海岸線信息的云圖問題二的解答問題分析風通常追蹤圖像塊(示蹤云)在三幅連續(xù)云圖中的移動來估算風速,并由示蹤云的高度指定風場的高度.其原理如圖 3 所示，左側(cè)為 t0 時刻的云圖，A 區(qū)域為示蹤云塊，右側(cè)為 t1 時刻云圖，在

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

23、叉相關(guān)系數(shù)的最大區(qū)域即為目標區(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 分別為示蹤云模板中的行與列號， m, n 1, 48為目標搜索區(qū)位移改變量， f (i, j) 為示蹤云模板的灰度值， g(i, j) 為目標搜索區(qū)云的灰度值，記 1 f (i, j ) 為示蹤云的平均灰度，g 1 g(i m, j n) 為目標搜索fN 2N 2ijij區(qū)的平均灰度，其中取 N 16 。傳統(tǒng)的 MCC 方法僅利用了圖像的灰度信

24、息，而在云圖上，相鄰區(qū)域的灰度具有空間連續(xù)性，在較小的圖像塊中灰度基本相同，很難分辨出云塊的紋理結(jié)構(gòu)特征。圖像的灰度梯度反映了相鄰區(qū)域的灰度變化，能夠提取圖像中目標物的輪廓信息，如圖 4 所示。因此，本文將圖像灰度的梯度信息加入相關(guān)系數(shù)的計算中，梯度信息的交叉相關(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 ）。當 J 取得最大時對應(yīng)的區(qū)域即為示蹤云下一時刻的位置，目標搜索區(qū)域內(nèi)相對于示蹤云在 x, y 方向上分別移動了m 和 n 個像素點，建立地球模型可計算出目標云相對于示蹤云的位移量，除以兩幅圖- 9 -像的觀測時間間隔到風矢量。(a)灰度(b)梯度圖 4 灰度云圖及其梯度對比4.2.3 二維風矢量計算方法假設(shè)需要求解起始點 A 點風矢的方向和大小，現(xiàn)已知 A 點的坐標(A,A) ，且通過相關(guān)的計算得出了下一時刻點 B 點的坐標(B ,B ) ；那么 A 點風矢的方向是從正北方向順時針旋轉(zhuǎn)到風的來向之間的角度（正北風的

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

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

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

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

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

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

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

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

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