2012年全國研究生數(shù)學建模競賽優(yōu)秀選

1、參賽（由 填寫）第九屆“杯”研究生數(shù)學建模競賽學校國防科學技術大學參賽隊號900020161. 徐海洋2.隊員3.參賽（由 填寫）第九屆“杯”研究生數(shù)學建模競賽題 目基于網(wǎng)絡流相關的云圖技術摘云圖計算要：的過程進行建模、求解和優(yōu)化；根據(jù)傳統(tǒng)本文氣象法存在的缺陷，提出了新的計算方法網(wǎng)絡流相關法并建立了像素窗口和搜索范圍的自適應調整模型。通過云圖實現(xiàn)的計算，需要建立灰度矩陣與的換算。云圖的灰度本文利用射影幾何基本原理，建立了基于幾何的簡單模型將矩陣到坐標的模型，并實現(xiàn)了在云圖上繪制海岸線等地理信息。本文首先用相關法建立了的基本模型，選取中心矩相關系數(shù)作為匹配的相關指標。由于傳統(tǒng)相關法利用局部點對點

2、的匹配方式，會造成一定的凌亂風矢出現(xiàn)，并且未對云圖信息進行充分利用。本文在相關法基礎上，提出了網(wǎng)絡流相關法，直接利用三幅云圖信息，采用周圍輔助采樣點對中心采樣點進行修正的方式，有效避免了傳統(tǒng)相關法缺陷，提高了計算精度。此外，本文通過實例計算與分析，得出了像素窗口大小和搜索范圍變化的對算法影響一般規(guī)律，并在此基礎上提出了可以自適應調整窗口大小和搜索范圍的方法，提高了匹配質量。同時本文還基于垂直大氣輪廓線方程，通過最小二乘擬合的方式建立了在衛(wèi)圖灰度與所在等壓面的模型，實現(xiàn)了對等壓面的估計。獲得了根據(jù)已建立的換算視場坐標模型，本文利用風云 2 號所拍攝的相鄰的云圖，對相關法以及本文提出的網(wǎng)絡流相關法

3、和自適應調整法進行了實例計算。其中相關法、網(wǎng)絡流相關法以及自適應調整法所得到非零數(shù)分別為 4491、4532 和 5376。并且通過比較和分析發(fā)現(xiàn)，本文提出的網(wǎng)絡流聯(lián)合相關法以及自適應調整法的風場連續(xù)性與一致性均優(yōu)于相關法，且 16×16 窗口1大小對于間隔 30 分鐘的云圖序列是較為合理，而 64×64 的搜索范圍則過大。最后，本文中的模式識別、旋轉匹配、背景灰度修正和輔助數(shù)據(jù)的應用等模型的拓展進行了初步的討論。：；相關法；最大費用最大流；模糊模式識別2目錄一、二、三、四、4.1問題重述4模型假設4符號說明5問題分析6技術64.2 視場坐標換算6的方法探究74.3.1 相

4、關法84.3相位分析法8法94.3.24.3.3 TCFM4.3.4 算法討論94.4 窗口大小與搜索范圍的.10量的等壓面高度10模型的建立與求解124.5五、5.1 換算視場坐標模型125.1.1 地球球面上的物點投影到視場平面上像坐標的換算135.1.2 灰度矩陣的像素行列號與視場平面上對應點的關系145.1.3 地球曲面上采樣點的坐標確定16模型175.2.1 相關系數(shù)法175.25.2.2 網(wǎng)絡流5.3 自適應相關法20.245.3.1 窗口大小與搜索范圍的進一步分析245.3.2 自適應匹配算法描述27等壓面的估計31結果分析325.4六、6.1 網(wǎng)絡流相關法與相關系數(shù)法326.2

5、 自適應匹配法與相關系數(shù)法35七、八、模型的評價38模型的拓展388.1 模式識別388.2 旋轉匹配398.3 背景灰度修正398.4 輔助數(shù)據(jù)的應用40九、參考文獻403一、 問題重述我國幅員遼闊、南北緯度差異大，在同一時間不同地點具有不同的氣象特征。氣象的變化對交通、工農(nóng)業(yè)生產(chǎn)等許多部門有著重要的影響。地球上平均有1/3到1/2的區(qū)域被云層覆蓋，云不僅對氣候的變化起著舉足輕重的作用，也是人類判別天氣情況的重要依據(jù)云圖是云層分布與運動的重要，并以其直觀性、準確性和信息的多樣性受到人們廣泛關注。利用氣象云圖獲得云跡風，也稱為，是根據(jù)對連續(xù)幾幅云圖中同一云塊的跟蹤，從云的運動估算出該區(qū)域的風速

6、、風向與所在高度。近幾十年來，利用氣象資料反演的已廣泛應用在臺風范圍、熱帶氣旋移向預報、暴雨和洪澇災害分析等方面，并越來越顯示其重要的作用。因此，利用是極為重要的。云圖對云跡風的問題 1：將云圖的視場坐標換算為，建立灰度矩陣與坐標的對應關系，將灰度矩陣繪制海岸線輪廓。圖，并根據(jù)海岸線的在云圖上添加問題 2：根據(jù)給出的三個相鄰時刻的云圖數(shù)據(jù)，建立該時刻場的度量模型，給出計算方案，并在給定的窗口大小與搜索范圍條件下計算整數(shù)處的二維，對采用的窗口大小和搜索范圍進行理解與評價。問題 3：分析調整窗口大小和搜索范圍對計算結果的影響，設計自適應調整窗口大小和搜索范圍的有效方法，建立場的度量模型和算法。并利

7、用獲得的。分析比較不同計算方法的優(yōu)劣，提出關模型和算法計算數(shù)于云圖數(shù)據(jù)量和處的二維場計算研究的建議。問題 4：根據(jù)全球溫度數(shù)值預報數(shù)據(jù)估算出用來推算出移動距離的云塊的環(huán)境溫度，由此溫度推算出量的等壓面高度。再據(jù)此高度進行量分層。根據(jù)計算結果，在云圖上添加。二、 模型假設假設 1：所接受到的紅外探測光線沿直線，并且在過程中不存在吸收、反射或折射。這樣云圖的拍攝過程可用幾何光學簡化模型建立和計算，且并影響計算結果精度。假設 2：地球為理想橢球，可以由地球的一個經(jīng)過南北極的橢圓截面繞南北極的連線旋轉而得到，且忽略地球表面的起伏對云圖的影響。假設 3：氣象相對于地球中心的位置保持不變。假設 4: 不考

8、慮地球和天空的其它熱源對云層溫度的影響。在測定時間范圍內(nèi)，地球表面的背景溫度保持不變，且對云層溫度無影響。事實上，在后面拓展部分可以看到地物以及海洋和大陸對本文模型中不考慮這些因素。云圖是有所影響，但由于相關地理數(shù)據(jù)不足，假設 5：不考慮云圖的失真、畸變，認為是完美的光學成像。實際上由于衛(wèi)4星距離地球較近，在拍攝過程中，不是遠場成像，所得圖片存在一定的畸變。但是由于是用兩幅或兩副以上圖片進行比較，像差與畸變會在一定程度消除。故可以認為不考慮這些成像誤差。假設 6：本文認為從任意間隔 1 小時內(nèi)的兩幅云圖可以用于的計算。假設 7：云圖中等壓面與溫度以及灰度的關系符合已知的統(tǒng)計數(shù)據(jù)，即可以通過已有

9、的統(tǒng)計數(shù)據(jù)估計示蹤云的氣壓高。假設 8：假設文中所討論的地區(qū)的風速不超過 120 km/h，也就是搜索窗口固定為64´ 64 像素，由于在氣象學中超過 117 km/h 定義為 12 級颶風，為出現(xiàn)概率較低，因此這個假設是較為合理的。天氣狀況，三、 符號說明符號符號說明(𝑥 𝑖,𝑗 , 0, 𝑧 𝑖,𝑗 )（x, y, z）𝜑 𝑖,𝑗 𝜃 𝑖,𝑗 f GGSS灰度矩陣中行列號為(i,j)的像素所

10、對應的球面物點的空間點坐標視場平面坐標灰度矩陣中行列號為(i,j)的像素所對應的緯度坐標灰度矩陣中行列號為(i,j)的像素所對應的云圖灰度矩陣元素行列號與對應云圖的灰度矩陣坐標坐標的相鄰時刻云圖的灰度矩陣像素的窗口矩陣相鄰時刻云圖像素的窗口矩陣云圖上目標模塊內(nèi)的像素灰度矩陣(i, j )GI ,J 云圖上匹配模塊內(nèi)的像素灰度矩陣G'(i, j)M ,N 目標模塊平均灰度GI , J 匹配模塊平均灰度G'M , N T P(p,q) OC TOCGf J溫度矩陣等壓面輔助采樣點的行列號交叉原點矩相關系數(shù)原點矩相關度總和 最大費用流模型背景修正矩陣5四、 問題分析4.1技術獲得云跡

11、指用連續(xù)幾幅（間隔一般小于 1 小時）靜止利用靜止氣象氣象圖像追蹤圖像上示蹤云模塊的位移，根據(jù)云圖的灰度估算出用來推算移動距離的云塊的環(huán)境溫度，由此溫度并通過已有統(tǒng)計數(shù)據(jù)或其它輔助計算出示蹤云模塊所代表的水汽特征所在的高度層次，以獲得這些層次上風的估計值，即定量的計算出云的移動速度及云所在的等壓高度。反演風場是風廓線儀等云圖應用的一個重要方面。盡管當前可采用無線電探空儀、來獲得風場，但觀測頻次和探測范圍仍有較大不足，特別在海洋上，風場的觀測缺乏，目前已成為一種重要的用云圖來觀測和反演風場就顯得尤為重要。云跡風1。它作為陸地觀測網(wǎng)常規(guī)風測量的補充資料，對于數(shù)值預報初始場的構造、臺風、高空急流及水

12、汽輸送的研究等都具有非常重要的意義。但必須指出的是，由于只能掃描到以及沒有被覆蓋的中低云，。而當大氣運動的則失去了意義。對因此，在客觀上不能得到被覆蓋的云層的作用力沒有使云產(chǎn)生運動或者云自身沒有運動時，于全球風場資料的獲得風場互相補充。獲得的風場需要與各種其它探測工具獲得的目前所能得出的主要為（包括風行與風速）以及所在的氣壓高，為了直觀與觀察方便，一般以標有圖標的云圖的方式給出云跡風信息。由于云圖上的一個像素代表地球上的幾公里或者十幾公里，因此的精度是決定強風的量級和影響范圍、暴雨的落區(qū)等的預報準確度的關鍵因素。而由于天氣分析與的實時性的范圍內(nèi)的。方法的計算復雜度也是需要在可接受在實現(xiàn)的過程需

13、要主要涉及幾個問題：視場坐標與換算、具體算法的實現(xiàn)的與相關參數(shù)的設定以及高度的確定，以下對這幾個問題進行初步的分析。4.2 視場坐標換算為了實現(xiàn)的計算，首先需要實現(xiàn)云圖灰度矩陣與的換算。而云圖的拍攝過程就是對地球的表面云塊進行射影變換（即中心投影變換），變換中心為拍攝云圖的氣象，投影面為的任意視場平面，如下圖，其中ABC 是地球表面的云塊，則 A' B 'C ' 即為視場平面上云塊的對應的像，其中視場平面垂直于與地球中心的連線的。6SBCACBA圖 1對地球表面的中心投影變換在此基礎上，就可以通過簡單的幾何知識建立與云圖上的坐標（即灰度矩陣的行列號）對應關系。而具體的計

14、算的方式可以分為兩種不同的方式：將云圖的灰度矩陣的像素行列號到坐標，或者將經(jīng)緯度坐標這兩種到像素矩陣行列號方法。雖然在理論上互為逆，本質上是等價的，但在實際計算中存在很多差異。前一種方法，每一個像素都能到一個坐標，對于在云圖上繪制海岸線等地理信息在理論上不存在像素誤差，但需要進行搜素匹配，降低了計算效率。而后一種方法的推導過程更加復雜，雖然對于云圖繪制地理信息不需要匹配搜索，計算更為簡便，但結果可能會由于計算精度問題而導致像素誤差。本文中所考慮的采樣點較為少，因此為了保證計算精度選用將灰度矩陣的像素行列號到坐標的方法進行坐標的換算。4.3的方法探究資料反演云跡從上世紀70年應用靜止氣象始的，早

15、期的云跡風方法是采用類似的基反演是通過人工識別追蹤云團的移動，而目前主要的于局部云塊匹配的思想，基本原理見圖。以云圖上的一個采樣點為中心選取一個指定的（一般為方形）的目標模塊（搜索窗口），即示蹤云，在下一時刻云圖的同樣位置附近指定的搜索范圍（虛線框）內(nèi)尋找與目標模塊相匹配的模塊（匹配模塊），然后通過兩個模塊的匹配信息計算風速與風向。而對于零的，中心的標準一致，均是將氣壓高度大于950毫巴的物象以及風速小于4m/s的像元認為是地物2。7搜索區(qū)域匹配模塊目標模塊圖 2 局部云塊匹配原理的計算當前，云圖反演風場的具體方法有很多種，主要包括相關法、相位分析法與TCFM法。4.3.1 相關法相關法是目前

16、開發(fā)和使用最為成方法，利用相關法處理常規(guī)云圖序列計算云跡風在美國、用。、歐洲、澳大利亞等和地區(qū)得到了廣泛的應相關法的原理是通過計算目標模塊的灰度矩陣與搜索范圍內(nèi)所有模塊的灰度矩陣的相關度或者相似程度，然后選取相關度最大的模塊作為匹配模塊，通過兩個模塊間的像素差（即相對位置）計算。按相關度的計算方式的不同分為絕對值法、序列相似性檢驗法和交叉相關系數(shù)法等方法。由于這些方法僅僅是計算相關度的公式不同，因此在計算速度和精度上沒有實質的上差異。其叉系數(shù)法由于數(shù)學理論依據(jù)較為充分，而被各國氣象部門廣泛應用。由于相關法計算風速時是整像素數(shù)，對于位移不足一個像素的云團，或者在連續(xù)的間隔小于5分鐘的快速區(qū)域云圖

17、上，目標模塊的位移一般小于目前靜止紅外云圖上一個像素的尺度，就會出現(xiàn)“亞像素尺度位移”的問題3,4。由于相對于常規(guī)云圖而言，快速區(qū)域云圖僅提高了云圖的時間分辨率，而云圖的空間分辨率并沒有提高，因此，此時傳統(tǒng)的相關法在風速計算上會產(chǎn)生很大的誤差。雖然可以通過擬合的來估計匹配度的峰值，但這種方法不僅在相關法計算量已經(jīng)很大的基礎上，進一步增加計算量，且并不能有效地解決問題。4.3.2相位分析法5教授在研究1問題時，創(chuàng)新性地提出了分鐘時間分辨率區(qū)域云圖中的“亞像素尺度位移”相位分析技術，對示蹤云進行頻域波譜分析，由諧波的相位變化計算波速，最后由平均波速計算云跡風，美國GOES1分鐘間隔云圖進行探索性試

18、驗表明，這一方法不僅可以追蹤位移小于一個像素的8示蹤云塊，獲得云跡，而且可以同時獲得速度譜及其方差。這種方法更適合于短時間間隔的快速區(qū)域云圖，成為技術上對傳統(tǒng)的相關技術的一種補充。該系統(tǒng)生成的云跡風場可以滿足實際業(yè)務和科研工作的需要。相位法實質是目標模塊與下一時刻的相同位置模塊做頻譜匹配，通過相位差計算，因此，當云塊位移大于兩個像素時，計算過程中將出現(xiàn)相位重疊現(xiàn)象，計算的云塊的相速度和真實速度不一致。所以對于時間間隔較長的云圖序列，尤其是30分鐘至60分鐘間隔的常規(guī)云圖，相位法不適用。4.3.3 TCFM法由于時間間隔較長的云圖序列對于相位分析法則不適用，而相關法計算位移又會由于像素分辨率問題

19、導致風速計算存在誤差。為了綜合利用相6with位分析和最大相關兩種方法的優(yōu)點，提高普通間隔云圖序列的精度方法的概念相結合，提出 TCFM（Tracking cloud等研究者們將兩種Combined Fourier phase analysis andum correlation）方法。該方法首先用傳統(tǒng)相關法進行模塊匹配，然后通過匹配模塊與目標模塊的匹配度是否發(fā)生云塊的亞像素位移，若發(fā)生亞像素位移，則進一步利用相位分析法對目標模塊和匹配模塊進行頻譜分析，根據(jù)相位差計算亞象素位移，計算得到亞像素對傳統(tǒng)匹配法得到的進行修正。該方法有效解決了長時間間隔云圖亞像素位移的問題，提高了云跡風的精度，但由于

20、在相關法的基礎上加入了快速變換，計算復雜度較高。4.3.4 算法討論由于計算復雜度較高，為了縮小匹配搜索范圍，一般采用二次匹配法：一般情況下先將張紅外云圖原分辨率減半，計算云圖上給定模塊的位移，稱為粗匹配；以此粗匹配作為云跡風初值來縮小搜索區(qū)，用原分辨率云圖重新搜索，稱為精匹配，從而得到比粗匹配精度高的云跡少，僅為 6400 個，因此二次匹配法不是必要的。由于本文所選取的示蹤云較目前大部分的均是采用三幅相鄰時刻的云圖進行風計算的。事實上，兩幅相鄰時間的云圖提供信息已經(jīng)足夠用于計算，如何有效的利用三幅云圖的信息是值得討論的。目前各國氣象部門主要是采用前兩圖得到的云跡對后兩圖云跡進行質量的方法。歐

21、洲 ESA 也采用二次匹配，但使用紅外和可見光兩種云圖，即用紅外圖匹配得到的云跡風場（稱為粗匹配風場）作為初始風場，進行可見光圖匹配（稱為精匹配）。而實際上，這樣做并沒有有效地利用三幅圖的灰度信息，而且目前大多數(shù)方法是基于單點的模塊匹配。因此當計算出現(xiàn)錯誤時，無法對錯誤的方式進行修正。因此提出能夠有效利用三幅圖，并建立與周圍采樣點之間的精度有很大幫助。的方法，對提94.4 窗口大小與搜索范圍的搜索窗口大?。K大小）的選取是影響算法速度與精度的一個重要因素。如果選取的搜索窗口太小，雖然會有效地提高計算速度，但是由于過小的示蹤云缺乏足夠特征，會導致所得雜亂無章，出現(xiàn)較大誤差。而如果選取的搜索窗口

22、過大，所得的精度會有所提高，但由于示蹤云的特征太復雜，使得極有可能無法在下一時刻云圖上進行像素匹配或匹配度過低，而且過大的搜索窗口會極大的影響計算速度。Lunnon 和 Lowe(1992)在用指出相關法處理 云圖時選取模塊大小為 16×16 像素時，得到的風場要比由 32×32 的像元矩陣得到的風場效果更好。而且，在計算速度上選擇 16×16 像素的目標模塊要比選擇 32×32 像素的目標模塊快兩倍。目前主要是采用 16×16 的固定窗口大小，而自適應的根據(jù)云塊特征或者區(qū)分度來調整窗口大小。同樣，搜索范圍的大小對模型計算也會產(chǎn)生重要的影響。顯

23、然搜索范圍設越小，計算越快，反之范圍設置越大會更精確，但是延長了計算時間，但是由于風速與云圖間隔時間的限制，過大的搜索范圍反而導致模型對模塊匹配的錯誤 從而降低精度。題目中給出的搜索范圍限定為 64×64 個像素，云塊大小限定為 16×16 個像素，可以通過計算得知相應的最大風速限制為 120 km/h。一般的都是根據(jù)最大風速使用固定大小的搜索范圍，但由于云圖中大部分云塊運動速度遠小于最大風速，因此使用相同的搜索范圍是沒有必要的。而如果能給出根據(jù)搜索過程中的匹配度進行自適應的調整搜索范圍的方法，毫無疑問在保證精度的前提下會極大地提高計算效率。此外，采樣點的選取也十分重要。采

24、樣點太少了毫無疑問會降低云跡的信息量，而太多會必然導致搜索范圍的重疊過多，對計算帶來不必要的麻煩，而且信息冗余。一般采用隔 1 到 2 個取一個點進行匹配，這里本文采用每個整數(shù)作為采樣點進行計算。4.5量的等壓面高度除了給出過等壓面高度的大小、方向，還必須推算出高度。量的等壓面高度，通示蹤云等壓面高度在氣象學中的指定方法有很多，如近似法、氣候學法、暗影法、立體法、窗區(qū)通道溫度法、多通道測云高度法、多通道法和中的雙通道法等。本文采用窗區(qū)通道溫度法確定示蹤云高度。利用所給數(shù)據(jù)擬合的灰度與溫度關系圖，可以看出溫度與灰度不是成簡單的線性關系，由文獻可知符合云圖的灰度對應的是紅外線輻射量，而地球表面輻射

25、量與溫度關系近似黑體輻射公式，因此在一般情況下，需要根據(jù)灰體輻射規(guī)律進行灰度與溫度的擬合于插值。但由于本題已經(jīng)給出每個灰度對應的是紅外輻射溫度值，因此在本題中不需要對灰度和對應溫度進行擬合。10圖 3 灰度與溫度的關系曲線而溫度與氣壓高度同樣不是呈現(xiàn)一般的線性關系，根據(jù)文獻資料表示，當氣壓高于 200 HPa 時為對流層，熱量絕大部分來自地面，氣溫隨著高度的增加而遞減。溫度與等壓面成對數(shù)關系：T = a + b ln P而利用所給的溫度與氣壓高度數(shù)據(jù)，可以畫出溫度與氣壓高度關系曲線。下圖為任意選取了十個不同獲得的溫度與氣壓的關系曲線。圖中可以看到，當氣壓高于 200 HPa 時基本符合上式。當

26、氣壓低于 200 HPa 時，溫度與氣壓的關系不符合上述的指數(shù)關系，而在此范圍內(nèi)通過溫度推導壓強也會得到兩個結果， 這顯然是不符合實際的。而且氣壓低于 200 HPa 時，處于平流層、中間層和熱成層，這幾個區(qū)域中幾乎不存在云層，因此可以利用溫度與等壓面的對數(shù)關系對風矢的等壓面進行估計7。圖 4 溫度與氣壓的關系曲線由于統(tǒng)計數(shù)據(jù)中的溫度與等壓面是離散的，對于中間溫度需要插值。溫度與等壓面的關系式為對數(shù)關系，因此需要作對數(shù)最小二乘擬合，而且等壓面小于200 HPa 時的溫度統(tǒng)計數(shù)據(jù)不符合對數(shù)規(guī)律，因此僅使用等壓面高于 200 HPa 的數(shù)據(jù)。11五、 模型的建立與求解5.1 換算視場坐標模型為實現(xiàn)

27、根據(jù)云圖計算云跡風，首先要確定灰度矩陣像素中每個元素對應的采樣在地球上的建立從灰度矩陣像素，本文根據(jù)之前分析部分的討論，通過空間的換算模型。的方法云圖對應的灰度矩陣為G = Gij ，其中i, j Î1, 2288。首先建立空定義間直角坐標系，為計算簡潔，以地球中心為原點，地球中心與連線為 y 軸，南北極軸線為 z 軸，赤道平面為 x-y 平面，建立三維直角坐標系。記在該坐，其中 D = 42164000 ，為同步標系下的坐標為(0, D, 0)離地球中心的高度，如圖。圖 5 坐標系圖像所拍攝的圖像投射到的平面，也就是說任意垂直于 y視場平面是指軸（即與地心的連線）平面均可作為的視場

28、平面，其中靠近地球的這一側的平面均為虛像面，另一側均為實像面，這里為了計算方便取 x-z 平面（為虛像面）。125.1.1 地球球面上的物點投影到視場平面上像坐標的換算視場（x-z 平面）的一個像點坐標如何對應到地球曲面上相坐標。首先討論將應物點所對應的記視場平面上的某個虛像點 Q 的坐標為(x, 0, z)( i , j )(i, j )對應的地球曲面上物點 P 坐標為(x,y,z)。則經(jīng)過下形式：與像點 Q 的直線方程如x - x(i, j )= z - z(i, j )yD=(1)0 - x(i, j )0 - z(i, j )化簡可以得到xyz1-= 1-D(2)x(i, j )z(i

29、, j )地球可視為一個理想橢球，這個理想橢球可以由地球的一個經(jīng)過南北極的橢圓截面繞南北極的軸線旋轉得到，則地球球面的橢球方程可以表述為x2 + y2 + z2= 1(3)a2a2c2其中a = 6378136.5為橢圓截面的長半軸（赤道半徑）， c = 6356751.8為短半軸（極半徑），與上述直線方程聯(lián)立可求解得到與視場平面像點的連線在橢球面的上交點，其中一個即地球表面上對應的物點。ì x222+ y + z= 1ï a2a2c2(4)íxyzï1-= 1-Dïîx(i, j )z(i, j )消去變量 x 和 z 可得x22z

30、y2y(i, j )+)(1-)2 += 1( i , j )(5)a2c2a2D為后面表述與計算方便，這里記22xz(i, j )(i, j )K = (+)a2c2上述的二次方程存在兩個解，也就說理論上球球面上兩個點，y 坐標如下在一個方向上可以拍攝到地K / D ±K / D2 +1/ a2 - K / a2K / D2 +1/ a2y± =13但是注意到由于地球球面的遮擋，實際上所能拍攝到的云圖只有距離衛(wèi)星更近那一點。所以，所求對應物的點只能取 y 值更大的解。如圖中只能看到坐標為 𝑥+, 𝑦+,𝑧+ 的點，點

31、19909;, 𝑦,𝑧 被遮擋。圖 6 理論上可拍攝到的地球表面上的兩個點位置示意圖顯然有 y+ ³ y- ，則取K / D +K / D2 +1/ a2 - K / a2y+ =(6)K / D2 +1/ a25.1.2 灰度矩陣的像素行列號與視場平面上對應點的關系為得到灰度矩陣中元素對應的地球，這里首先確定灰度圖上每個像素在該坐標系下視場平面（即 x-z 平面）對應的坐標，再根據(jù)視場平面與地球曲面點坐標的關系，推導出像素行列號與坐標的關系。但是僅依靠空間坐標、視場平面坐標以及像素行列號無法求解出相應物點在地球曲面上的實際經(jīng)緯度坐標，必須根據(jù)所給的灰

32、度上某已知參照物（如地球）的相關像素信空息，確定像素行列號與視場平面 x-z 上坐標對應關系。本文通過視場平面間坐標，根據(jù)灰度矩陣計算出地球投影在視場平面上的橢圓實際尺度。地球在視場平面上的整體射影范圍，顯然當方程判別式為零時K - K1+= 0(7)D2a2a2直線與橢球相切，此時D2K0 = D2 - a2(8)故當且僅當視場平面上的像點的坐標滿足22xzD2(i, j )+(i, j )£ K=(9)0D2 - a2a2c214與視場平面的連線與地球曲面存在交點，也就是說，在該范圍內(nèi)的像點對應的物點是在地球表面上的。對上述不等式做如下變形有22xz(i, j )+(i, j )

33、£ 1(10)( K a)2( K c)200上述不等式所表明所能觀測到的地球曲面上的點所對應的視場平面像點，都必須落在 x-z 平面上的滿足不等式(10)的橢圓區(qū)域內(nèi)，也就是說該不等式描述的區(qū)域就是地球在視場平面上的像，如下圖的陰影區(qū)域。陰影橢圓邊界為與地球橢球的切點組成。其中2K0 a 和2K0 c 分別為視場平面上地球對應的像橢圓的在 x 方向上的長軸與 z 方向上的短軸。zyx22xz(i, j )(i, j )+£ 1( K a)2( K c)200圖 7所能觀測到的地球曲面上的點在視場平面上的投影區(qū)域由于所得到的灰度圖中的宇宙背景灰度為-1，故當像素的灰度值為

34、0 或者正數(shù)時，表明該像素位于地球對應的像橢圓（即圖中的陰影）上。則根據(jù)灰度的分布范圍可以知道地球對應的像素分布情況，從而可知道像橢圓長軸短軸在灰度圖像上的對應的像素長度，根據(jù)所給數(shù)據(jù)分別為 2169 像素與 2163 像素。由于像素大小都是相同的，因此由射影幾何可以知道像素大小對應的該坐標系下的尺度都是相等的。建立像素與選定視場平面的坐標尺度的關系：ìï2K0 a = 2169´ width(11)íïî2 K0 c = 2163´ height其中 width 與 height 分別為灰度圖上像素在該坐標系下視場平面15

35、（即 x-z 平面）對應的x 方向和z 方向的坐標尺度。計算可得ì2width=K0 aïï2169í(12) K0 aïheight= 2ïî2169則有灰度矩陣的元素與該坐標系下視場平面（即 x-z 平面）對應的坐標ìïx(i, j) =(1145 - j) × width(13)íz=(1145 - i) ×heightïî (i, j )5.1.3 地球曲面上采樣點的坐標確定根據(jù)上面求得的空間坐標關系，通過下列反三角函數(shù)即可將地球球面上點的與緯度分

36、別為q(i, j) 和空間坐標換算成地球上坐標，記像素(i,j)對應的j(i, j ) 。ìï y =K / D +K / D2 +1 / a2 - K / a2K / D2 +1 / a2ïï z= ( 1 - 1 )z(14)í y(i, j )yDïï x = ( 1 - 1 )xï y(i, j )yDîìjzz / y= arctan= arctanï(i, j )x2 + y2(x / y)2 +1ï(15)í= arccot x + dïq&

37、#239;î(i, j )y其中， d 是由d = m - 90o ， m 為不為 0 且處于東經(jīng)位置時所帶來的修正，有位置的，在本題目條件下 m = 86.5o 。當本模型應用于所處的其它位置時，由于東西經(jīng)的不同，需進行不同的修正。由函數(shù)組(12)-(15)，即可建立由灰度矩陣像素到f : (i, j) ® (q(i, j),j(i, j) )的(16)：本文利用上面得到的502 列三個交叉點處對應的以得到在灰度矩陣中第 500 行和第 500、501、，其他詳細信息見附件 jwd.txt 中。表格 1 通過視場坐標換算得到的數(shù)據(jù)16灰度矩陣行列號（500,500）（50

38、0,501）（500,502）(46.197, 32.719)(46.282, 32.713)(46.366, 32.708)此外，本文還通過上述和簡單的搜索匹配，實現(xiàn)了在云圖上繪制海岸線等地理信息，效果見下圖。可以從圖中看到我國國土輪廓線，澳大利亞輪廓線，輪廓線，以及長江黃河等明確的地理特征，說明以上算法精確可靠。圖 8繪制了海岸線的云圖5.2模型文章前面已經(jīng)確立了射 f。在此基礎上，本文對模。云圖灰度矩陣元素行列號與對應坐標的映進行數(shù)學建模，首先采用傳統(tǒng)的相關法進行建5.2.1 相關系數(shù)法(q ,j ) 的采樣點，利用第 1 問計算得到的對于給定的關系 f: (I, J ) ®

39、(q,j )f(16)通過簡單的搜索匹配，可得到灰度矩陣中相應像素的位置(I , J ) 。定義大小為1716×16 像素的窗口矩陣 S ，S 是灰度矩陣G 的子矩陣，且S9,9 為窗口的像素中心，S9,9 = GI ,J 。228816162288模塊矩陣 S云圖灰度矩陣 G圖 9 云圖灰度矩陣與模塊矩陣的關系以及模塊矩陣中心點的確定如上圖中所示，一方格為一個像素，中間為 16×16 像素大小的搜索窗口，并將它定義為矩陣 S，將第 9 行第 9 列的像素作為窗口的中心。S 為 2288×2288 像素大小的云圖灰度矩陣 G 的子矩陣。在相鄰一時刻的相鄰云圖，獲得

40、云圖的灰度矩陣G' ，G' 中選取的匹配中心位置為(M , N )相應的窗口矩陣為S ' ，同理S ' 是灰度矩陣G' 的子矩陣，S '為窗口，9,9的像素中心， S '= G'。9,9M ,N 取 搜 索 范 圍 為 64×64像 素 ， 從 圖 中 可 以 看 到 匹 配 模 塊 滿 足 ，M ÎI - 24, I + 24， N ÎJ - 24, J + 24。18642×24166416(I,J)目標模塊搜索區(qū)域圖 10 相關法中搜索區(qū)域的確定本文這里選用交叉原點矩相關系數(shù) OC（C

41、ross-Correlation Coefficient of Original Moments） 進行目標模塊匹配。定義如下：ååG()(i, j )i, j G'OC(I , J , M , N ) =I ,JM ,N(17)ååéG()2åå()2éùi, j ù'Gi, jëûëûI ,JM ,N(i, j )和G'(i, j) 分別表示間隔時刻獲得的其中，G云圖上目標模塊I ,J M ,N 和匹配模塊內(nèi)的像素灰度矩陣 S、

42、S '；OC(I, J , M , N) 為匹配位置在(I , J ) 處的i, ji, j模塊 S與在(M , N ) 的模塊 S ' 的原點矩相關系數(shù)。i, ji, j另外也可用中心矩相關系數(shù) CC（Cross-Correlation Coefficient of Central Moments）進行目標模塊的匹配。定義如下：åå éëG(i, j ) - Gé(i, j ) - G'ùù'GûëM ,N ûI ,JI ,JM ,NCC(I , J , M ,

43、 N ) =(18)ùù22ééååéG(i, j ) - GI ,J û ååM ,N(i, j ) - G'ù'GëêëëM ,N ûúûI ,J其中， G和G'為模塊內(nèi)像素的平均灰度，但CC(I, J , M , N) 僅反映模I , J M , N 塊內(nèi)圖像強度分布的表面起伏特征的相似程度，并不反映平均強度的變化。OC(I, J , M , N) 的大小，既含有圖像強度表面起伏幅度及其

44、時間變化的信息，又與平均強度的變化有關。需要指出的是，OC(I, J , M , N) 和CC(I, J , M , N) 在作用上并無優(yōu)劣之分，只是從不同角度反映了模塊內(nèi)的圖像變化信息。但是由于 OC 和 CC 的計算中涉及到除法，在實際計算中會出現(xiàn)分母為零的情況。例如，采用中心矩相關系數(shù)法進行相關計算，若匹配模塊的十六個灰度是完全相同，則分母為零。而原點矩相關系數(shù)法同樣存在這個問題，只要模塊的灰度全部為 0，則也出現(xiàn)分母為零的情況。但是在實際測試中，出現(xiàn)灰度為零的概率很小，出現(xiàn)一個模塊中十六個灰度相同的情況是有可能的。因此，本文選擇交19(I,J)2×24叉原點矩相關系數(shù) OC

45、進行目標模塊匹配。關于風速的計算采用普遍的近似算法，用目標模塊和匹配模塊間的像素數(shù)目乘以像素分辨率除以云圖間隔時間計算得到。方向以云圖上的正北方向進行標定。而由之前的分析部分提高，根據(jù)者等壓面大于 950HPa 時，認為是地物，取氣象中心的標準，風速小于 4m/s 或為零。在本題中可以近似認為像的分辨率 5 公里/像素。事實上素分辨率為矢，所以云圖中靠云圖主要用來測量橫向風部分的誤差較小，所以的分辨率在變化較小的范圍內(nèi)是可以接收的。5.2.2 網(wǎng)絡流相關法文中上述的相關法是采用點對點的匹配方式來進行計算，但是有的時候計算結果中會存在較大的誤差。如下圖中所示，顯然中間點的風向與的周圍點的風向存在

46、明顯的差別。由于風的變化一般都是連續(xù)的，中間點的風向顯然是不合理的，為凌亂。事實上，這主要是由于傳統(tǒng)相關法只是局部的點與點匹配，對于缺乏特征的云塊，很容易在搜索范圍內(nèi)找到大量的具有相同特征的高匹配度云塊，從而大大影響獲得的的正確性。而該方法沒無法利用點周圍的相關信息，使得相鄰點之間的匹配沒有相互影響。因此沒有辦法利用周圍的點對這種相似度匹配的錯誤進行修正，通過這種點與點的相關度的匹配，必然會出現(xiàn)匹配錯誤。圖 11 亂流示意圖目前對于上述相關法出現(xiàn)的缺陷，各國氣象部門常采用后兩幅云圖得到來調整或的對前兩幅得到的進行質量，即利用客觀規(guī)律或刪除認為不合理的，最終生成高質量的云跡風風場。這樣可以在一定

47、程度上進行修正，但是實際上還是沒有有效地利用到三幅圖，且還是沒有利用到周圍點的特征，質量只是將相鄰時刻的進行對比從而估計每個點所得的度，并不能對錯誤進行正確的排除或修正。因此，提出如何有效利用三幅云圖信息，并且能有效利用周圍點相互關系對傳統(tǒng)相關法缺陷進行修正的算法是具有現(xiàn)實意義和利用價值的。考慮到尋常的測量中會通過重復多次測量結果取平均值，以此為測量到的量，減小測量誤差。在此處，為了減少測量誤差，我們同樣采取多次取樣測量的方式，而由于到云塊的變化以及的連續(xù)性，可以在采樣點周圍選取若干個輔助的采樣點，通過對輔助采樣點的匹配信息和對所選采樣點的進行修正。本文法進行闡述。這種缺陷提出了網(wǎng)絡流相關法，

48、下面對具體模型和實現(xiàn)方20對于采樣點(I,J)，選取與(I,J)相鄰的八個采樣點(I-1,J-1)，(I-1,J)，(I-1,J+1)，(I,J-1)，(I,J)，(I,J+1)，(I+1,J-1)，(I+1,J)，(I+1,J+1)作為輔助采樣點。分布如下圖所示。圖 12 中心采樣點與周圍輔助采樣點的分布示意圖八個輔助采樣點分別以這幾個采樣點為中心，按與之前相同的規(guī)則選取目標窗口，大小均為 16×16 。對應的灰度矩陣記為Sp,q ，其中p ÎI -1, I +1, I, q ÎJ -1, J +1, J在相鄰一時刻的相鄰云圖G' ，G' 中匹配

49、中心位置為(M , N)相應的窗口矩pq，陣為S '，其中M Î p - 24, p + 24， N Î+ 24p,q pq。設上述的九個采樣點中任意兩個采樣點( p, q) 和( p', q' ) 在下一時刻對應的中心位置不能重合，即(Mp , Nq ) ¹ (Mp' , Nq' ) 。根據(jù)文章前面內(nèi)容可知，原點矩相關系數(shù)的計算公式為ååG()(i, j )'i, j Gp,qM p ,NqOC( p, q, M , N ) =pq2ååéëGp,q

50、29;åM p ,Nq(i, j )ùû2(i, j )ùé'Gëû為了充分利用云圖信息和周圍輔助點的信息，文章不采用單個采樣點的相關系數(shù)法，而是采用九個輔助采樣點的原點矩相關度。定義總原點矩為TOC，定義如下TOC(I , J ) =åOC( p, q, M p , Nq )(19)pÎI -1,I ,I +1qÎJ -1,J ,J +1且(M p , Nq ) ¹ (M p' , Nq' )TOC 即為九個采樣點獲得的原點矩相關系數(shù)的總和。TOC 越大，表明

51、目標21GI-1,J-1GI-1,JGI-1,J+1GI,J-1GI,JGI,J+1GI+1,J-1GI+1,JGI+1,J+采樣點與周圍輔助采樣點在下一時刻云圖上的相對位置越接近，所獲得的越。下相近。取 TOC 最大的情況下所獲得的九個的平均矢量為采樣點的圖方法中如何利用八個輔助采樣點進行修正的示意圖。111213111213212223212223313233313233算法無誤差221112131112132122232123313233323331實際情況13 相關系數(shù)法計算圖如果匹配算法中可能出現(xiàn)的情況出現(xiàn)凌亂，則文章利用九個采樣點得到的與傳統(tǒng)的采用一個采樣點得到的結果相同。但是如果

52、出現(xiàn)圖中的實際情況，如果采用傳統(tǒng)的算法，得到的錯誤為n。但是如果采用本文的網(wǎng)絡流相關為利用八個輔助采樣點修正后的n ' ，顯然該最終的法，最終與真實的更為接近。因此，采用本文的方法可以得到更為精確的實際于傳統(tǒng)的算法。，結果明顯優(yōu)求這 9 個采樣點的顯然是一個最優(yōu)多點匹配問題，而處理該方法的典型思路即把它轉換為最小費用最大流模型。在有三幅不同時間的云圖（記為圖A，B，C）的情況下，由于幾乎在短時間內(nèi)發(fā)生突然反向的變化，所以可以同時參考多幅云圖構建如下圖論模型：考慮一個邊上帶權的有向圖G =VAB ,Vc, E，VAB，Vc 分別對應網(wǎng)絡流中抽象的點，其中 Va 表示灰度矩陣 B 中所選取的九個采樣點與下面搜索范圍內(nèi)= V (M , N , p, q)。V = V (M ' , N ' ) 表示 C 圖上選取所有點配成的點對，VABpqcpq的與 B 圖中選取的匹配點(M ' , N ' ) 。pqE= V (M , N , p, q),V (M ' , N ' ),W ( p, q, M ' , N ' )，ABCpqpqBCpq其中W ( p, q, M ' , N ' ) = OC( p, q, M ' , N ' ) ，為 BC 云圖中( p, q) 與(M