不同尺度DEM對面雨量估算模型的影響分析
——以武夷山為例

葉棟水 林彬彬

(1.福建省氣象信息中心，福建 福州 350001；2.南京信息工程大學地理科學學院，江蘇 南京 210044)

1 概述

高時空分辨率的山地降水空間分布對于山體滑坡等自然災害管理、山區水資源的開發與利用都具有重要意義[1-2]。在宏觀的大氣環流背景下，山地降水會受地形的影響，存在著氣流的動力抬升、降水的再分配以及局地環流等情況，空間分布十分復雜[3-4]。目前，降水空間分布的模擬方法主要有基于氣象站點降水數據的空間插值[5]、遙感觀測[6-7]、數值模擬[8-9]、再分析[10]等。其中，空間插值法對氣象站空間密度依賴性較高，數值模擬存在參數不確定性的問題，遙感觀測降水以其較精準的捕捉大尺度降水空間分布能力，得到了廣泛應用。但由于遙感降水產品空間分辨率較大，無法滿足小區域地形降水對模擬精度的要求，因此國內外學者[11-17]在遙感降水產品的基礎上，引入地形因子，提出了各類型的山地面雨量模型，獲得了更高分辨率、更精準的山地降水空間分布。

面雨量模型所依賴的坡度、坡向、海拔等地形因子，在不同精細程度的地形下存在差異[18-20]，進而會影響模型模擬精度。這樣的地形尺度效應一直以來都是降水模擬過程中的所關心的科學問題。吳昊旻[21]、劉金波[22]、陳躍浩[23]、高學杰[8]等討論了不同水平分辨率氣候模式的精度差異；衣娜娜[24]、黃海波[25]、杜娟[26]等探討了不同空間分辨率對WRF模擬結果的影響。這些研究均表明，降水的模擬過程對地形分辨率的敏感性較為復雜，在精細的地形下，并不一定獲得最理想的降水模擬效果。

不論是氣候模式還是WRF模擬，都是對宏觀地區降水的模擬方法，所分析的地形分辨率大多為30km、50km等，最精細的也只是討論了3km、5km以及9km分辨率的地形對降水模擬的影響[24]。若要精細化模擬小區域的山地降水，需要分析更高分辨率的地形對山地面雨量模型的影響。目前發布的1km及1km以內的地形數據主要有1km分辨率的SRTM30_PLUS、500m分辨率的SRTM15_PLUS，90m分辨率的SRTM3，30m分辨率的NGA-SRTM1、ASTER-GDEM-V2、NASA-SRTM1，15m分辨率的SRTM-X-DLR和12m分辨率的ALOS PALSAR。但這些地形數據的空間分辨率不連續，無法使用這些數據來連貫地討論地形尺度效應對山地面雨量模型的影響。

對此，本文使用二維經驗模態分解算法，獲取相同分辨率但不同精細化程度的地形，即為不同尺度地形；并以武夷山作為研究區域，分別建立各尺度地形下的山地面雨量模型，初步探索地形尺度效應對面雨量模型的影響，旨在挑選最合適的地形尺度，提高面雨量模型的分辨率和模擬精度，為山地自然災害管理和山地水資源開發提供更可靠的數據基礎。

2 研究方法及數據來源

2.1 研究區域與數據

武夷山脈位于江西和福建邊界，是兩省水系的分水嶺，其東部以丘陵為主，西部和中部以山地為主，地勢西北較高，東南較低；山脈呈東北西南走向，長約530km，海拔介于-36m到2153m之間(如圖1中的紅色邊界研究區所示)。該地區屬于中亞熱帶濕潤氣候，年平均氣溫在13～14℃，年降水量在2000mm以上，是福建省降水量最多地區。復雜的武夷山脈地形對福建降水具有重要影響，不僅可以阻滯冷空氣南移，同時還對暖冷氣流強迫抬升，形成顯著的區域內地形降水。本文選取24°N～29°N，115°E～120°E中的一塊區域作為研究區域，初步探討地形尺度效應對面雨量模型的影響。

圖1 武夷山及周邊地區高程及氣象站點分布

(1)氣象站觀測數據。本文使用的2000年至2017年7月份逐日降水數據下載自國家氣象科學數據中心(http://data.cma.cn/)，已對數據進行了質量控制，并將日數據累加至月數據。本文共使用175個氣象站點的降水數據，其中118個外圍氣象站點的降水數據只用來輔助計算研究區內63個氣象站點的模型因子(主導降水方位)，不參與模型構建和驗證統計。另外，研究區內的63個氣象站點中，57個氣象站作為模擬站點，6個氣象站作為驗證站點如圖1所示。

(2)DEM數據。本文使用的DEM數據為2000年的SRTM(Shuttle Radar Topography Mission)數據，下載自地理空間數據云(http://www.gscloud.cn/)，空間分辨率為90m×90m，采用的投影坐標系為Albers投影(中央子午線為110°，標準緯線分別為25°和47°)。

(3)衛星降水產品。本文使用的2000年至2017年TRMM 3B43V7的7月衛星降水數據，下載自美國NASA(http://mirador.gsfcnasa.gov)，空間分辨率為0.25°×0.25°。

2.2 研究方法

2.2.1 DEM尺度轉換算法

經驗模態分解(Empirical Mode Decomposition，EMD)是N.E.Huang[27]等人提出的一種數據驅動的自適應非線性時變一維信號分解方法，隨后J.C.Nunes等人[28]將其發展為處理二維信號的方法(two-dimensional empirical mode decomposition,BEMD)，在圖像降噪、壓縮和融合方面有較廣泛的應用。BEMD是將原始信號分解成有限個頻率從高到低的二維BIMF分量以及每個BIMF所對應的r余量，其中BIMF分量通常為原始信號中不同頻率的信號噪聲，而余量r為一單調函數，一般為原始信號中不同頻率的信號趨勢。

從信號學的角度來講，數字高程模型(Digital Elevation Model，DEM)也可以視為二維的隨機信號，高程值即為信號數值，一系列頻率從高到低的二維BIMF分量即為地形中的精細部分，所對應的一系列r余量即為由微觀到宏觀、精細到粗略的地形趨勢，即為多尺度地形。設x和y分別是橫縱坐標，原始DEM數據為q(x,y)，j為分量BIMF和余量r的標號，N為分解次數，則二維經驗模態分解如式(1)所示：

(1)

分量BIMF和余量r的具體計算過程如下：

①初始化：令j=1，g0(x,y)=q(x,y)；

②剝離第j個BIMF分量，并獲得第j個r余量：

1)令f0(x,y)=gj-1(x,y)，動態變量i=1；

2)在fi-1(x,y)上進行3×3窗口滑動，計算fi-1(x,y)的局部極大、極小值，組成極大極小值點集；

3)根據極大極小值點集，雙線性內插成極大值包絡面ti-1(x,y)和極小值包絡面bi-1(x,y)，得到平均包絡面mi-1(x,y)：

mi-1(x,y)=(ti-1(x,y)+bi-1(x,y))/2

(2)

4)根據公式(3)計算fi(x,y)：

fi(x,y)=fi-1(x,y)-mi-1(x,y)

(3)

5)根據公式(4)計算SD值[23]：

(4)

其中，X為DEM的行數，Y為DEM的列數。

6)SD經驗值為0.2～0.3，本文取0.3。若SD<0.3，則代表了精細地形的分量BIMFj=fi(x,y)，代表了宏觀地形的余量rj=gj-1(x,y)-BIMFj。 若SD≥0.3，則將fi(x,y)看作一個新的信號，并令i=i+1，轉至步驟2)，重復執行步驟2)～6)，直至滿足終止條件。

本文分解次數N為5，即得到了5個分量r1至r5，分別對應5級尺度地形。包含原始DEM在內，共為6級尺度地形。其地形分辨率均為90m×90m，每個柵格數值的含義均為高程，單位均為m。

2.2.2 面雨量模型

地面降水由空氣中的水汽含量、大氣運動等大氣方面的主導因素和海拔、迎背風坡等地形方面的影響因素共同決定。因此，月面雨量的理論估算模型為P′=P0+P1，其中P0為面雨量模型的趨勢項，和大氣緊密相關；P1為模型地形項。Zhu等[13]采用海拔高度、坡向、坡度和主導降水方位的函數表達模型地形項，最終得到分站分月面雨量模型，如式(5)所示：

P′=k0P0+k1h+k2cos(PPD-β)·h+k3sin(PPD-β)·h+k4cos(PPD-β)·sin2α+Δa

(5)

其中P′在建立模型過程中，為地面氣象站點觀測得到的降水量，在根據模型進行面雨量模擬的過程中，P′為所估算的面雨量，空間分辨率為90m×90m；P0為面雨量模型的趨勢項，可用衛星降水數據代替，本文采用TRMM 3B43；k1h+k2cos(PPD-β)·h+k3sin(PPD-β)·h+k4cos(PPD-β)·sin2α即為P1模型地形項，其中h為海拔高度；PPD為主導降水方位；β為坡向；α為坡度；k0-4為各項系數；Δa為余項。

PPD是Zhu等[13]提出的用來指示水汽來源方向的參數。其具體計算方法如下：

1)設定A′氣象站為目標氣象站，其所在位置的坡向為β′(-180°為北，-90°為東，0°為南，90°為西)；設定方位δl為-180°(其方位坐標系與坡向的方位坐標系一致)；

2)搜索A氣象站鄰近的n個氣象站A0～An，計算這些氣象站所在位置的坡向β0～βn。計算β0～βn與δl的夾角，將夾角≤90°所對應的氣象站歸為迎風坡站，反之為背風坡站；

3)統計迎背風坡氣象站的平均降水量差DPl；

4)δl按間距Δδ順時針旋轉，重復2)、3)，直至回到-180°。

5)在所獲得的DPl集中選取最大值，PPD即為最大值所對應的方位δl。

本文n從3～60分別進行嘗試，最終選取40時，最符合水汽來源的自然規律；Δδ選取7.5°。

2.3 技術路線

不同尺度地形的面雨量模型由面雨量模型的地形項和趨勢項構成，其中地形項的計算步驟為：①首先，將原始90m×90m空間分辨率的DEM數據按照本文2.2.1所提到的BEMD經驗模態分解算法進行多尺度地形轉換，得到5級尺度地形r1～r5，并計算得到包含原始DEM在內的6級尺度地形的高程、坡度以及坡向；②融合研究區及周邊區域175個氣象站點的月平均降水數據，獲得6級地形尺度下建模站和擬合站的主導降水方位；③最后根據建模站月平均降水數據以及主導降水方位數據，計算得到面雨量模型的地形項。面雨量模型的趨勢項計算步驟為：將空間分辨率為0.25°×0.25°的TRMM 3B43V7數據進行雙線性內插，重采樣至90m×90m空間分辨率的TRMM數據，該數據即為模型的趨勢項。將建模站月平均降水數據、模型地形項、模型趨勢項代入模型中，進行建模，每個建模站每月獲得一套模型參數。

將建模站的模型參數插值至整個研究區，并根據研究區的地形項和趨勢項，模擬得到面雨量。最后，根據擬合站的面雨量模擬結果以及氣象觀測得到的月平均降水數據，進行模型的驗證與精度評價。技術路線如圖2所示。

圖2 多尺度地形的面雨量模型技術路線

3 結果分析

由于7月水汽充足，降水豐沛，選取7月作為研究月份，并將氣象站2000年至2017年7月的降水量取多年平均作為參照數據，反距離加權插值到研究區上，得到氣象站觀測降水量空間分布圖，如圖3(a)所示。同樣對2000年至2017年7月的TRMM 3B43V7數據取多年平均作為模型的衛星降水數據。由于TRMM 3B43V7數據分辨率較低，無法參與到90m分辨率的面雨量模型中，本文將其雙線性內插至90m×90m分辨率，以參與面雨量模型構建，如圖3(b)所示。

圖3 武夷山及周邊地區的氣象站與TRMM 3B43V7降水量90m×90m空間分布

3.1 各尺度地形

將90m×90m分辨率的DEM作為輸入項，根據BEMD算法，分解得到微觀到宏觀不同精細程度，分別為r1，r2，r3，r4，r5共5個連續變化的地形尺度，如圖4所示。圖4(a)是最為精細的原始DEM數據，細小的山脈邊界都較為明顯。r1剔除了原始DEM中最破碎精細的地形，圖4(b)代表高程較高的白色區域面積增大，意味著相比于DEM，r1所刻畫的山峰被削弱，更多區域被歸于山峰。圖4(c)中的r2進一步剔除在r1中相對較細小的部分，山脈邊界更加柔和。圖4(d)中的r3更加粗略，已經只能辨認清山脈主體的形狀和走勢，山脈邊界更加光滑，而圖4(e)和圖4(f)的r4和r5進一步綜合地形信息，r4西南角的兩座山峰(白色區域)在r5中已經相連，概化為一座山峰，已經無法呈現山脈的原始形狀。從DEM到r1，再逐級到r5，是地形的平滑過程，越往后分解，地形精細度丟失越大，宏觀分區趨勢越來越明顯。另外，從該區域高程的最大最小值分析，隨著地形尺度的變大，高程最大值在逐漸減小，最小值逐漸增大，表明地形精細部分在逐漸去除，符合地形尺度轉換的信息綜合理論，從而證明了BEMD應用于地形尺度轉換的可行性。

圖4 武夷山及周邊地區的各尺度地形

為進一步對比各尺度地形下，氣象站點所在位置海拔、坡向、坡度三個地形因子的區別，本文還進行了統計，其中坡向的坐標系為：北180°，東90°，南0°，西90°，平地181°。部分氣象站的數據如表1所示。

隨著地形尺度的增大，精細化程度的降低，每個站點的海拔、坡向、坡度都發生了變化。同一站點在地形尺度的增大過程中，并不表現為單調遞增或者單調遞減，不同站點的變化規律也不一致。同一站點的坡向在地形宏觀化的過程中，會存在較大的變化，甚至在58734站點的DEM向r1轉變的過程中，發生了215°左右的變化。這可能是由于所在位置的小山坡在地形轉化過程中被概化，其坡向數值從小山坡的坡向變為所在位置的更大尺度的山坡坡向。

3.2 各尺度地形的面雨量模型結果

將BEMD算法所得到的r1到r5地形以及原始的DEM數據作為6層尺度地形數據，分別在各尺度地形下，計算目標氣象站周圍10個氣象站的7月平均降水數據以及所在位置的坡度、坡向和高程，得到該氣象站的PPD，從而進一步計算該站的面雨量模型地形項。接著融合TRMM 3B43V7重采樣的90m分辨率衛星降水數據，構建目標氣象站的面雨量模型。最后，將所有建模氣象站的模型參數反距離加權插值至整個研究區，獲得研究區各尺度地形的面雨量模型模擬結果，如圖5所示。

表1 研究區各尺度地形下氣象站點的地形因子

圖5 武夷山及周邊地區的各尺度地形面雨量模型模擬結果

由DEM以及r1到r5共6個尺度地形構建的面雨量模型模擬結果和氣象站插值降水空間分布較一致，并且DEM以及r1到r3相較于氣象站插值降水空間分布以及TRMM 3B43V7重采樣的結果而言(圖3)，模型模擬的降水結果更能體現實際降水存在迎背風坡分布的規律。6級多尺度地形下模擬得到的降水空間分布遵從了多尺度地形的規律，也表現為從微觀到宏觀、精細到粗略的轉變。其中，圖5(a)中的DEM模型模擬結果細節最為豐富，可以體現細小山脈的迎背風坡降水量差異。圖5(b)中的r1模型模擬結果相較于圖5(a)稍微更加粗略，但在東南部，已將細碎的山脈降水分布情況進行了綜合，體現為從西南到東北走向的條狀分布。r2模型的模擬結果相比于r1，更加綜合，東南區域以及北部區域均出現更粗的條帶狀，山脈的迎背風坡降水差異體現的較為明顯。相比于r2模型，r3模型模擬結果對于北部區域的山脈迎背風坡刻畫更加綜合，但南部山脈的降水空間分布已經趨于均衡，從模擬結果已經無法辨認山脈的走勢。從r4、r5模型的模擬結果空間分布上看，已經無法辨認山脈的迎背風坡，只能呈現西北降水量小、東南降水量大的階梯狀降水空間分布。而且，模擬結果均出現由地形信息被過多去除而引起的模型模擬結果偏差，在圖5(d)和(f)中表現為幾何塊狀。

除此之外，隨著地形尺度的增大，研究區模型模擬得到的降水量最大最小值均不存在單調遞增或者單調遞減的規律。但DEM模型模擬結果的最大降水量值為584.90，最小降水量值為33.35，相比于氣象站插值所得到的最大最小降水量值(圖3(a))最大值更大，最小值更少，而且相差較多。而r2和r5模型結果的最大最小值更接近圖3(a)，且相比于r5，r2模型結果的空間分布更接近圖3(a)。

3.3 面雨量模型最優尺度地形選擇

對未參與建模的6個驗證氣象站點進行誤差分析，統計得到各尺度地形面雨量模型驗證站點的誤差統計量[12]：相關系數r2，平均相對誤差MRE和均方根誤差RMSE，誤差情況如表2所示。

表2 不同尺度地形面雨量模型驗證誤差評價指標

由表2可知，r1和r2面雨量模型對武夷山及其周邊地區7月份降水量模擬效果較好，相關系數均在0.90以上，r2面雨量模型精度最高，相關系數達0.96，r1次之。DEM、r3、r4和r5面雨量模型的相關系數均在0.9以下，相比于r2，下降了至少10.4%。DEM、r3、r4和r5面雨量模型的MRE和RMSE明顯增大，MRE均高于8%， RMSE均超過12.5mm，其中r5模型誤差最大，MRE、RMSE分別達到12.24%、18.32mm。結果表明，r2模型的模擬效果最優，r1模型略次，DEM、r3、r4和r5模型效果依次下降。

結合統計結果和圖5可知，DEM模型模擬結果的降水量最大最小值較極端(圖5(a))，這與模型所依賴的地形尺度直接相關。DEM為最精細的地形尺度，對細小的山脈均刻畫淋漓盡致。但對于面雨量模型而言，過于精細的地形不但不會提高模型精度，反而會給模型帶入噪聲，從而降低了模型精度。同樣，相比于尺度更小的r1模型，r2模型雖然基于更大尺度的地形，但精度更高。相比于r2模型，r3面雨量模型的模擬結果精度降低，這表明在r2至r3的變化過程中，對于降水空間分布具有重要意義的地形信息在綜合過程中被去除，從而影響了模型精度。在重要地形信息的不斷去除過程中，r4、r5模型模擬結果精確度逐漸下降。

為了更直觀展示以上驗證結果，將各地形尺度面雨量模型估算值和氣象站點觀測值繪制成散點圖，如圖6所示。

圖6 驗證站的各地形尺度面雨量模型估算值和氣象站觀測值散點圖

由圖6可知，6個尺度地形的面雨量模型估算的月降水普遍比氣象站觀測的月降水高，但其中有一個站點降水估算均過低，這表明不同尺度地形面雨量模型在同一站點的預測上，存在一定的規律性。另外，DEM、r1模型在150mm以下的降水量估測中劣于r2模型，這是由于150mm以下的降水量主要集中于西北山區，而這塊區域的地形相較于東南部分更加起伏，地形對降水的影響更大。對于這塊區域而言，地形尺度的改變對其影響較大，精細的DEM不僅不會提高模型精度，反而會加入噪聲，影響模型的準確性。綜上所述，各尺度面雨量模型中，r2尺度地形不僅最大程度規避了會產生噪聲的地形細節，而且保留了對面雨量存在影響的地形信息，對武夷山及其周邊地區7月份月降水量模擬的準確性表現最優。

4 結論

本文采用BEMD算法逐步去除90m×90m分辨率的DEM中的精細部分，得到5個尺度的地形，分別為r1、r2、r3、r4和r5。包含原始DEM數據在內，一共6個尺度地形分別參與面雨量模型構建，對武夷山及其周邊地區7月份降水進行模擬和驗證，初步探究地形尺度對面雨量模型的影響。結論如下：

①相比于r2尺度地形，DEM、r1精度較高，微觀地形刻畫更細致，但過于精細的地形信息在面雨量模型中引入噪聲，導致其精度劣于r2尺度模型。

②地形尺度增大到r3、r4和r5后，對降水空間分布具有影響的重要地形信息在綜合過程中逐步被去除，模擬精度逐漸降低。其中r4、r5面雨量模型的模擬結果出現與事實不符的幾何塊狀，無法適用于面雨量模型的模擬。

③研究區內，r2尺度面雨量模型最大程度規避了模型噪聲，保留了與降水相關的地形信息，模型精度最高，為最優尺度。

本文以地形復雜的武夷山及其周邊地區作為研究區，降水量充沛的7月作為研究月份，初步探索了不同尺度地形對面雨量模型的影響，在此基礎上可進一步討論不同月份、地貌地形(山地、盆地等)以及不同降水強度等因素條件下，尺度效應對面雨量模型的影響。