氣象遙感圖像及格點場重采樣插值方法

葉金印，邱旭敏，黃 勇，張春莉

1.淮河流域氣象中心，安徽 蚌埠 233040

2.安徽省氣象科學研究所，合肥 230031

氣象遙感圖像及格點場重采樣插值方法

葉金印1，邱旭敏1，黃 勇2，張春莉1

1.淮河流域氣象中心，安徽 蚌埠 233040

2.安徽省氣象科學研究所，合肥 230031

1 引言

氣象遙感圖像產品和氣象格點場資料在天氣預報業務和科研中發揮著不可替代的作用[1-2]，尤其是在氣象客觀自動化預報中備受關注[3-4]。氣象業務數據按表現形式分為兩類：一類是非規則的離散類數據，如氣象站點觀測數據；一類是規則的格點（柵格）類數據，如氣象遙感圖像及數值天氣預報格點場。在氣象業務尤其是科研中，常常需要對各類資料進行融合，以便于氣象成因分析。原始氣象遙感圖像和氣象格點場的分辨率和地圖投影方式往往不同[5]，若要對它們進行綜合分析，需要進行重采樣插值處理。對規則的格點（柵格）類數據進行重采樣插值是首要環節，也是重要的環節之一。氣象遙感圖像產品都是以某種地圖投影方式表現的，與在經緯網格坐標系中表現的氣象格點場無法直接匹配，需要將氣象遙感圖像重采樣插值到經緯網格坐標系。同時，以經緯網格坐標系表示的數值天氣預報產品（如降水預報格點場），也因數值預報模式不同而具有不同的分辨率和區域范圍[6-7]；除數值天氣預報格點場外，其他氣象格點（柵格）場也存在使用不同分辨率表示的問題[4]；為便于預報自動化和科學研究，也需要將不同分辨率的氣象格點場插值到相同分辨率的經緯網格坐標系中。

從不同研究目的和應用需求出發，對于圖像重采樣和插值方法研究有不同的側重，如王茂新等人針對氣象衛星NOAA/AVHRR數據的特點，提出了三次卷積法、雙線性內插法等幾種常規采樣算法[5]。樓繡林等人又針對重采樣速度慢的問題提出了一種重采樣的快速算法，利用“塊操作”快速確定重采樣圖像中任意像元點在原始圖像中的位置，同時提出了鄰點權重法代替其他插值方法[8]。一般圖像處理中的重采樣插值主要用于對圖像本身的放大（縮小）和幾何校正；對于氣象業務中使用的氣象遙感圖像，重采樣插值主要通過地圖投影坐標轉換來提取以及填補氣象信息。鄰點權重法是圖像產品插值方法中最為實用的一種[8]，而最近鄰點法[9]是在計算方法上與鄰點權重法最為接近的一種。因此，本文對不同分辨率的氣象遙感圖像產品，分別采用鄰點權重法及最近鄰點法進行重采樣插值，并對插值結果進行分析。

雙線性插值方法在圖像和格點場插值中的應用比較廣泛，且效果較好[10]。李得勤等人采用雙線性插值方法對數值預報模式產品進行降尺度處理，發現插值結果與實際觀測的均方根誤差以及相關系數和插值前相差甚小[11]。王岳山首次將貝塞爾插值方法應用在數值預報格點場插值中，指出采用貝塞爾方法進行插值不僅速度快而且精度也非常高[12]。因此，本文對不同分辨率的氣象格點場資料分別采用貝塞爾插值及雙線性插值方法插值，并對插值結果進行分析。

2 氣象遙感圖像重采樣插值方法

氣象遙感圖像重采樣插值的目的是將一個坐標系中的點陣信息在另一個坐標系中按一定格式“最逼真地”表示出來。基本方法有直接法和間接法，如圖1所示。直接法，是從某一種地圖投影坐標的原始圖像上的像元點坐標出發，按照地圖投影坐標轉換公式求出目標坐標系中像元點對應的坐標，然后將原始圖像上像元點(x，y)處的數據直接賦予或經過插值賦予目標坐標系(X，Y)處的像元點。間接法，是從目標坐標系的像元點坐標出發，按照地圖投影坐標轉換公式求出原始坐標系中像元點的坐標，然后將原始坐標系中像元點(x，y)處的數據直接賦予或經過插值賦予目標坐標系中(X，Y)處的像元點。

圖1 直接法和間接法示意圖

氣象遙感圖像重采樣插值關鍵環節有2個：地圖投影坐標向目標坐標系的轉換以及對像元點數據進行插值。本文以蘭勃托投影衛星云圖[13]為例，描述將其重采樣插值到經緯網格坐標系中的方法。

2.1 地圖投影坐標轉換方法

地圖投影坐標轉換[14-15]有兩種方式：（1）由目標坐標系（經緯網格坐標系）坐標值求出地圖投影坐標系（蘭勃托投影）中坐標值，稱之為正解，即由經緯度計算投影坐標：(φ，λ)→(X，Y)。

（2）從地圖投影坐標（蘭勃托投影）向目標坐標系（經緯網格坐標系）轉換，稱之為反解。即由投影坐標計算經緯度：(X，Y)→(φ，λ)。

2.2 像元點數據插值方法

（1）最近鄰點法[9]：這種方法以距離被計算點最近的一個像元的灰度值作為輸出像元的灰度值。僅通過距離比較，就近取值，僅有一個像元參與計算，因而計算效率高，但重采樣插值效果較差，灰度的連續性受到一定程度的破壞。

（2）鄰點權重法，其示意圖如圖2所示。

圖2 鄰點權重法示意圖

插值點p周圍4個原始像元點的灰度值都參與計算，根據距離d的不同，各個像元點灰度值的貢獻不同，距離小的像元貢獻大，距離大的像元貢獻小。由于“鄰點權重法”在決定插值點灰度值時距離起決定作用，能夠盡可能多地保持原圖像的信息。

兩種插值方法相比，鄰點權重法的計算量顯然要大于最近鄰點法的計算量。

2.3 氣象遙感圖像重采樣插值方法的實現

直接法和間接法不同之處在于投影變換中源坐標系和目標坐標系的選擇不同。坐標轉換后，對像元點的插值方法可以根據需要選取合適的方法，如最近鄰點法及鄰點權重法。

直接法實現步驟：首先建立一個經緯網格坐標系，利用地圖投影坐標轉換反解公式解析出原始圖像每個像元點的經緯度(φ，λ)，進而對經緯度坐標系上每一個網格點(φi，λi)進行插值。在插值方法上，最近鄰點法需要對原始圖像中所有像元點的經緯度進行判斷，來找出距離所求網格點(φi，λi)最近的像元點，將像元點的灰度值賦予該點；鄰點權重法則需要找出距離該點最近的周圍4個像元點，判斷過程復雜，計算量大。

間接法實現步驟：首先建立一個經緯網格坐標系，再從該經緯網格坐標系出發，利用地圖投影坐標轉換正解公式計算出每個網格點(φi，λi)在原始圖像中的坐標點(X，Y)，計算出坐標點(X，Y)周圍的4個像元點坐標，即原始圖像橫軸方向第i和第i+1及縱軸方向第j和第j+1的4個交叉點。進而對該坐標點(X，Y)進行插值，并將值賦予網格點(φi，λi)。在插值方法上，最近鄰點法僅需找到距離坐標點(X，Y)最近的像元點，并將值賦予網格點(φi，λi)。鄰點權重法需要找到距離坐標點(X，Y)周圍的4個像元點，并按照距離權重對灰度值進行插值計算。

比較直接法和間接法的計算過程，間接法重采樣插值的計算量遠小于直接法。

3 氣象格點場插值方法

與氣象站觀測數據（不規則離散數據）不同的是，數值天氣預報模式輸出的客觀分析場和要素預報場均以一定經緯網格距的格點場表示，不同數值天氣預報模式輸出產品具有不同的經緯范圍和分辨率。針對業務和科研的需要，常常要把不同分辨率的數值預報模式輸出產品進行標準化，形成一定經緯度范圍內分辨率相同的氣象格點場。本文以數值天氣預報模式輸出的規則格點場為例，比較分析雙線性插值方法和貝塞爾插值方法的插值效果。

3.1 雙線性插值方法

雙線性插值，又稱為雙線性內插。在數學上，雙線性插值是有兩個變量的插值函數的線性插值擴展，其核心思想是在兩個方向分別進行一次線性插值。假如想得到未知函數f在點P=(x，y)的值，假設已知函數f在Q11=(x1，y1)、Q12=(x1，y2)、Q21=(x2，y1)以及Q22=(x2，y2)四個點的值。首先在x方向進行線性插值，得到R1和R2，然后在y方向進行線性插值，得到P。雙線性插值示意圖如圖3所示。

圖3 雙線性插值示意圖

3.2 貝塞爾插值方法

將貝塞爾公式[12]以中央差分格式寫成插值公式，只保留到第二階差分：

其中，t=(x-x0)/h，h是格距。

y-1、y0、y1、y2分別為同一條直線上相鄰的4個點x-1、x0、x1、x2的值，y為位于x0與x1之間的點x的值。將貝塞爾插值方法應用到氣象格點場插值中，需要在兩個方向上分別進行一次插值：先在x方向上進行貝塞爾插值，再在y方向上進行貝塞爾插值，得到最終結果。

4 應用分析

4.1 氣象遙感圖像重采樣應用分析

本文從氣象業務科研需求出發，對直接法和間接法重采樣插值的效果進行分析，重點分析評估易于實現的間接法中最近鄰點法和鄰點權重法的插值效果。

以氣象業務中應用的FY2E和FY2D蘭勃托投影紅外衛星云圖為例，進行重采樣插值應用分析。FY2E紅外云圖像元點陣為512×512，像素點間距為13.0 km，投影中心點經緯度為：（110°E，30°N），左下角經緯度為：（86.59°E，0.64°S）；FY2D紅外云圖像元點陣為1 200×1 200，像素點間距為4.9 km，投影中心點經緯度為：（100°E，35°N），左下角經緯度為：（77.32°E，6.59°N）。

（1）直接法重采樣插值應用分析：圖4（a）及圖4（b）分別為FY2E蘭勃托投影紅外云圖（2011年12月7日8時，GMT+ 8）原始圖像和采用直接法進行地圖投影坐標轉換后的圖像（Fortran SGL繪圖，窗口分辨率為1 200×700，下同）。FY2E紅外云圖經過直接法地圖投影坐標轉換后共有512× 512個像元點。由圖4可以看出，在高分辨率目標坐標系中，低分辨率的FY2E紅外云圖經過地圖投影轉換后，像元點之間有很明顯的間隙（灰度值為空）。

圖4 （a）FY2E紅外云圖

圖4 （b）經過直接法投影轉換后的圖像

圖5（a）及圖5（b）分別為FY2D（2012年8月9日8時15分，GMT+8）蘭勃托投影紅外云圖原始圖像和采用直接法進行地圖投影坐標轉換后的圖像（Fortran SGL繪圖）。FY2D紅外云圖經過直接法地圖投影坐標轉換后共有1 200×1 200個像元點。由圖5可以看出，目標坐標系分辨率與原始圖像分辨率基本一致的情況下，高分辨率的FY2D紅外云圖經過地圖投影坐標轉換后的圖像并沒有顯現出間隙。

比較圖4（b）與圖5（b）發現，對于高分辨率目標坐標系而言，低分辨率的圖像相對于高分辨率的圖像，采用直接法重采樣時存在更多的需要插值的“間隙”。低分辨率圖像因坐標轉換后的像元點間隙明顯，最近鄰點法的插值會有較大誤差；而鄰點權重法可以在一定程度上減小這種誤差。高分辨率圖像因坐標轉換后的像元點間隙不明顯，最近鄰點法的插值誤差與鄰點權重法的誤差比較接近。由此可見，隨著圖像分辨率的提高，最近鄰點法在重采樣插值效果上會更接近于鄰點權重法，且因為計算簡單，在計算量上會小于鄰點權重法。

圖5 （a）FY2D紅外云圖

圖5 （b）經過直接法投影轉換后的圖像

（2）間接法重采樣插值應用分析：根據FY2E和FY2D紅外云圖分辨率及投影位置分別建立與之對應的重采樣經緯網格坐標系：FY2E重采樣區域選為10°N～50°N，90°E～130°E，網格間距0.1°；FY2D重采樣區域選為10°N～50°N，80°E～120°E，網格間距0.05°。分別采用最近鄰點法及鄰點權重法對其進行重采樣插值。

圖6為FY2E紅外云圖最近鄰點插值與鄰點權重插值結果在Fortran SGL中的顯示。其中，圖6（c）和圖6（d）分別為圖6（a）和圖6（b）局部特征放大圖像。比較圖6（c）和圖6（d）可以看出，低分辨率的FY2E紅外云圖采用兩種插值方法得到的圖像在紋理和清晰度上存在明顯差異。

圖6 FY2E紅外云圖插值結果在Fortran SGL中的顯示

圖7為FY2D紅外云圖最近鄰點插值與鄰點權重插值結果在Fortran SGL中的顯示。其中，圖7（c）和圖7（d）分別為圖6（a）和圖6（b）局部特征放大圖像。比較圖7（c）和圖7（d）可以看出，高分辨率的FY2E紅外云圖采用兩種插值方法得到的圖像在紋理和清晰度上未表現出差異。

圖8為FY2E和FY2D紅外云圖的最近鄰點法插值與鄰點權重法插值結果在MICAPS（氣象信息綜合分析處理系統）[16]中的疊加顯示。其中，圖8（c）和圖8（d）分別為圖8（a）和圖8（b）局部特征放大圖像。可以看出：低分辨率的FY2E紅外云圖的最近鄰點插值與鄰點權重插值結果等值線有明顯不吻合，鄰點權重法插值結果與原始圖像更為接近；而高分辨率的FY2E紅外云圖的最近鄰點插值與鄰點權重插值結果等值線幾乎吻合，且與原始圖像完全吻合。

圖7 FY2D紅外云圖插值結果在Fortran SGL中的顯示

圖8 最近鄰點法（紅線）與鄰點權重法（黑線）插值結果在MICAPS中的疊加顯示

為定量評估兩種插值方法的效果，引入清晰度進行分析。清晰度[17]一般使用平均梯度衡量，公式為：

式（4）中F為平均梯度，G(x，y)為影像函數，m，n為影像的行列數。根據相對誤差來分析方法之間效果的差異，公式為：

式（5）中f為相對誤差，v1為比較圖（數據）參數，v2為原始圖（數據）參數。

表1為最近鄰點法及鄰點權重法清晰度統計分析結果。對不同分辨率的圖像，最近鄰點法相對于鄰點權重法的清晰度相對誤差分別為14.171%（FY2E）和1.849%（FY2D）。對于低分辨率圖像（FY2E），最近鄰點法在重采樣插值效果上與鄰點權重法有較大差距；而對于高分辨率圖像（FY2D），兩者差距不大。

表1 最近鄰點法及鄰點權重法清晰度

由此可見，隨著圖像分辨率的提高，最近鄰點法在重采樣插值效果上會更加接近鄰點權重法，且最近鄰點法計算方法簡單、計算量小。在客觀預報自動化業務中，對于高頻次、高分辨率遙感圖像，重采樣插值采用最近鄰點法更為實用。

4.2 氣象格點場插值結果分析

以氣象業務中日常應用的歐洲中期天氣預報中心（ECMWF）降水預報場（2012年8月13日20時，GMT+8）為例進行插值方法分析。歐洲中期天氣預報中心降水預報場是間距為0.25°的等經緯距格點場，區域范圍為0°E～180°E、0°N～90°N。

分別采用貝塞爾插值和雙線性插值將其插值為網格距為0.2°的等經緯距格點場，區域范圍取在114°E～120°E、29°N～35°N。圖9為兩種插值結果與原降水預報場在MICAPS系統中疊加顯示圖。貝塞爾插值結果與原始降水場等值線基本吻合（圖9（a）），而雙線性插值結果與原始降水場等值線吻合度不如貝塞爾插值（圖9（b））。

圖9 ECMWF降水預報場（紅線）與兩種插值結果（黑線）疊加顯示圖

表2為降水預報原始場以及采用貝塞爾和雙線性插值的統計分析結果，貝塞爾插值結果相對于降水原始場，標準差相對誤差0.640%，平均值相對誤差0.226%；而雙線性插值標準差相對誤差3.025%，平均值相對誤差0.605%。可以看出貝塞爾插值結果與雙線性插值結果均較為接近原始場，貝塞爾插值結果優于雙線性插值。

表2 兩種插值方法結果統計

為進一步分析兩種插值方法的效果，對原始降水預報場等間距抽點，形成分辨率為0.5°的格點場，分別采用貝塞爾插值方法及雙線性插值方法將其插值到與原始場分辨率相同的0.25°的格點場。以0.25°原始場為參照，對比分析采用不同插值方法所生成的0.25°格點場。圖10為兩種插值結果與原始降水預報場在MICAPS系統中疊加顯示圖。由圖10（a）及圖10（b）對比可以看出，貝塞爾插值與雙線性插值效果差別比較明顯，這種差別在降雨中心區域表現得更加明顯。

圖10 ECMWF降水預報場（紅線）與兩種插值結果（黑線）疊加顯示圖

表3為降水預報原始場以及經抽點后采用貝塞爾和雙線性插值得到降水預報場的統計分析結果。貝塞爾插值結果相對于降水原始場，方差相對誤差2.225%，平均值相對誤差0.546%；而雙線性插值方差相對誤差4.242%，平均值相對誤差0.664%。從統計結果看，可以同樣得到貝塞爾插值方法的效果優于雙線性插值方法。

表3 兩種插值方法結果統計（抽點試驗）

5 結論

以氣象業務中不同分辨率的氣象衛星（FY2E和FY2D）遙感圖像以及歐洲中期天氣預報中心（ECMWF）降水預報場為例，分別對不同重采樣插值方法進行了分析比較，得出如下結論：

（1）基于間接重采樣的氣象遙感圖像最近鄰點插值法的計算量小于鄰點權重插值方法，而鄰點權重插值方法的效果優于最近鄰點插值方法。

（2）隨遙感圖像時空分辨率的提高，對于高頻次、高分辨率的圖像，采用基于間接重采樣的最近鄰點法在計算速度方面優勢更加明顯，在重采樣插值效果上更加接近鄰點權重法，能快速準確地完成圖像重采樣插值，更實用于客觀自動化預報。

（3）對于數值天氣預報模式輸出的降水預報場，采用貝塞爾插值方法得到的插值結果更接近于原始場，貝塞爾插值方法的插值效果優于雙線性插值方法。

（4）本文采用的氣象遙感圖像與氣象格點場重采樣插值方法具有很好的實用性，便于計算機實現和業務化。

[1]閔錦忠，沈桐立，陳海山，等.衛星云圖資料反演的質量控制及變分同化數值試驗[J].應用氣象學報，2000，11（4）：410-418.

[2]許小峰.中國新一代多普勒雷達網的建設與技術應用[J].中國工程學，2003，5（6）：7-14.

[3]黎光清，董超華.衛星氣象遙感和應用：現狀、問題、前景[J].氣象科技，1990，18（1）：1-7.

[4]王玉斌，萬玉發，羅兵，等.一種等值線填充并行算法[J].計算機工程與應用，2012，48（28）：61-65.

[5]王茂新，沙奕卓，于莉.關于NOAA AVHRR圖像重采樣及投影方法的研究[J].中國圖象圖形學報，1997，2（1）：38-42.

[6]楊學勝，胡江林，陳德輝，等.全球有限區數值預報模式動力框架的試驗[J].科學通報，2008，53（20）：2418-2423.

[7]梁利，林開平，黃海洪.幾種數值預報模式在廣西降水預報效果的比較[J].氣象研究與應用，2012，33（2）：1-4.

[8]樓繡林，黃韋艮，周長寶，等.遙感圖像數據重采樣的一種快速算法[J].遙感學報，2002，6（2）：97-101.

[9]鮑文東，楊春德，邵周岳，等.幾何精校正中三種重采樣內插方法的定量比較[J].測繪通報，2009（3）：71-72.

[10]管志杰，趙政.地圖投影變換及其在GIS中的應用[J].計算機工程與應用，2000，36（6）：50-52.

[11]李得勤，陳力強，周曉珊，等.風電場風速降尺度預報方法對比分析[J].氣象與環境學報，2012，28（6）：25-31.

[12]王岳山.貝塞爾插值及其在數值預報中的應用[J].海洋學報，1995，12（2）：12-17.

[13]許建民，張文健，楊軍，等.風云二號衛星業務產品與衛星數據格式實用手冊[M].北京：氣象出版社，2008.

[14]Coordinate conversions and transformations including formulas[EB/OL].[2012-12-26].http：//www.docin.com/p-25267998. html.

[15]廖宇.一種基于NSCT分解的圖像插值算法[J].計算機工程與應用，2012，48（8）：176-178.

[16]李月安，曹莉，高嵩，等.MICAPS預報業務平臺現狀與發展[J].氣象，2010，36（7）：50-55.

[17]何力，孫涵，黃永璘，等.MODIS1B數據的重采樣方法研究[J].遙感信息，2007（3）：39-44.

YE Jinyin1,QIU Xumin1,HUANG Yong2,ZHANG Chunli1

1.Huaihe River Basin Meteorological Center,Bengbu,Anhui 233040,China
2.Anhui Institute of Meteorology,Hefei 230031,China

Resampling interpolation method is one of the problems in the field of meteorological information processing research. This paper introduces a direct and indirect resampling interpolation methods which are based on map projection coordinate conversion for the meteorological remote sensing images,and introduces bilinear interpolation method and Bessel interpolation method for the meteorological grid point field.Taking meteorological satellite（FY2E and FY2D）remote sensing images with different resolutions and ECMWF（European Centre for Medium-Range Weather Forecasts）precipitation forecast field for example,the paper analyzes different resampling interpolation methods.The results show that the calculation amount of the nearest neighbor algorithm is less than that of the weighted nearest neighbor algorithm based on indirect resampling method,while the weighted nearest neighbor algorithm can get better results than the nearest neighbor algorithm.With the improvement of resolution,the comparative advantage of calculation amount of the nearest neighbor algorithm is more obvious.The weighted nearest neighbor algorithm is more suitable for high-resolution meteorological remote sensing images.The results also show that Bessel interpolation algorithm is better than bilinear interpolation algorithm for the meteorological grid point field.

meteorological remote sensing images;meteorological grid point field;resampling;interpolation method

重采樣插值方法是氣象信息處理領域研究的問題之一。針對氣象遙感圖像，介紹了基于地圖投影坐標轉換的直接重采樣插值方法和間接重采樣插值方法；針對氣象格點場，介紹了雙線性插值方法和貝塞爾插值方法。以氣象業務中不同分辨率的氣象衛星（FY2E和FY2D）遙感圖像以及歐洲中期天氣預報中心（ECMWF）降水預報場為例，分別對不同重采樣插值方法進行了分析比較。結果表明：基于間接重采樣的氣象遙感圖像最近鄰點插值法的計算量小于鄰點權重插值方法，而鄰點權重插值方法的效果優于最近鄰點插值方法；隨著圖像的分辨率提高，最近鄰點插值法與鄰點權重插值方法相比，計算量小的優勢更加明顯；對于高分辨率的氣象遙感圖像建議采用基于間接重采樣的最近鄰點法；對于氣象格點場，貝塞爾插值方法的插值效果優于雙線性插值方法。

氣象遙感圖像；氣象格點場；重采樣；插值方法

A

TP311

10.3778/j.issn.1002-8331.1302-0211

YE Jinyin,QIU Xumin,HUANG Yong,et al.Resampling interpolation methods of meteorological remote sensing image and grid point field.Computer Engineering and Applications,2013,49（18）：237-241.

國家自然科學基金（No.41275030，No.41105098）；淮河流域氣象開放研究基金（No.HRM201103）。

葉金印（1968—），男，高級工程師，主要從事水文氣象學及數值分析的研究。E-mail：yejinyin@sina.com

2013-03-01

2013-05-21

1002-8331（2013）18-0237-05

CNKI出版日期：2013-06-08 http://www.cnki.net/kcms/detail/11.2127.TP.20130608.0953.005.html