基于DEM的京津冀地區地形起伏度分析

張 競，杜 東，白耀楠，紀冬麗

(1.中國地質調查局 天津地質調查中心，天津 300170； 2.天津城建大學 市政與環境工程學院，天津 300384)

起伏度最佳分析窗口的大小因地形而異：國際地理學會地貌制圖委員會在編制《歐洲1∶250萬地貌圖》時使用了16 km2的分析窗口[6]；陳志朋等[7]在《1∶400萬中國及其鄰區地貌圖》的編制中使用了21 km2的分析窗口；王玲等[8]基于1∶25萬DEM計算的新疆山地最佳分析窗口為2.56 km2；張偉等[9]基于SRTM DEM，認為黃土高原最佳分析窗口為3.88 km2，云貴山區為4.58 km2，青藏高原東緣為5.34 km2；胡最等[10]基于1∶25萬DEM確定的湖南省最佳分析窗口為0.28 km2；王讓虎等[11]基于ASTER GDEM提取的我國東北地區最佳分析窗口為2.62 km2；丁賢法[12]基于SRTM DEM得到的云南省富寧縣最佳分析窗口為0.52 km2；郎玲玲等[13]基于1∶25萬DEM得到福建低山丘陵區的最佳分析窗口為4.41 km2。

既然最佳分析窗口具有明顯的空間差異性，當研究區幅員遼闊、地貌類型復雜時，就可能存在不止一個最佳分析窗口，如涂漢明等[4]指出中國存在5種不同規模的地形起伏度最佳分析窗口——2、10、16、20、22km2，適用于不同的地形地貌。但是，目前的大量研究都采用全局單值最佳分析窗口，其中不乏面積廣大、地形復雜的地區[3,6,10-11,14-15]，此時的最佳分析窗口僅能夠合理地表達一部分地形，對其他地形則失去了“最佳”屬性。本研究選取地貌類型復雜的京津冀地區，針對不同地貌類型提取地形起伏度最佳分析窗口，分析地形對最佳分析窗口的影響，提出京津冀地區地形起伏度的合理表達。

1 研究區概況

京津冀地區位于華北平原北部，戰略地位十分重要。研究區地理范圍介于113°27′31″～119°51′33″E、36°02′47″～42°36′53″N，總面積21.54萬km2，地質構造條件復雜，發育有平原、丘陵、盆地、山區、高原等多種地貌類型，北靠燕山山脈，南部平原展布，西倚太行山，東臨渤海灣，西北部為壩上高原，丘陵主要分布于太行山東側和燕山南側，沿渤海海岸多灘涂、濕地，總體呈現出西北高東南低的地形特點。區內海拔2 000 m以上的山峰有10座，第一高峰位于張家口小五臺山東臺，海拔2 882 m[3]。

2 研究數據與研究方法

2.1 研究數據

本研究DEM數據使用90 m分辨率SRTM3 DEM，絕對高程精度為±16 m，絕對平面精度為±20 m[16]，數據獲取自地理空間數據云網站。京津冀三地行政區劃矢量數據獲取自國家基礎地理信息中心，比例尺為1∶25萬。所有數據統一采用我國CGCS2000_GK_CM_117E投影，在該投影下DEM實際分辨率為84.8 m，使用雙線性插值法重采樣為90 m。利用矢量數據對DEM數據進行裁剪后得到京津冀地區DEM數據(圖1)。

圖1 京津冀地區DEM及研究樣區

2.2 起伏度和最佳分析窗口提取原理

地形起伏度一般基于DEM數據采用窗口分析法求取，其表達式為

Ra=Hmax-Hmin

(1)

式中：Ra為以第a個柵格為中心的一定面積窗口內的地形起伏度；Hmax和Hmin分別為該窗口內的最大高程和最小高程。

起伏度概念的關鍵是最佳分析窗口，在該窗口下提取的地形起伏度可以滿足山體完整性和區域普適性原則。山體完整性原則是指分析窗口起初只包含一部分山體，隨著面積增加，窗口內高差迅速增大，當窗口恰好覆蓋整個山體時，高差達到一較大值，之后隨窗口面積增大高差增大的速度明顯減小，那么這個分析窗口之內的高差就恰到好處地反映了該山體的地形起伏特征。區域普適性原則是指最佳分析窗口需要滿足研究區域內高差最大山體的完整性，此時的分析窗口必然滿足較小山體的完整性原則[1-2]。以往學者在提取最佳分析窗口時使用了不同的方法，如人工作圖法[6,13,17]、最大高差法[13]、最大高差-面積比法[3,18]、模糊數學法[1]、均值變點法[8-9,11,14,19-26]等。趙斌濱等[27]對各種方法的精度進行了驗證，認為人工作圖法和均值變點法得到的最佳分析窗口對應的起伏度更符合實際情況，最大高差法和模糊數學法已不適用。由于人工作圖法的主觀性較強，近年來均值變點法成為最佳分析窗口計算中的主流方法，這種方法對恰有一個變點的檢驗最為有效[28]。本研究采用均值變點法計算地形起伏度的最佳分析窗口，計算原理如下。

(1)依次計算N個遞增分析窗口下單位面積上的平均起伏度，即單位地勢度T。

1.2.1 對照組 對患者每個月開展1次肺健康知識宣教，內容包括COPD相關知識、肺康復理念，指導戒煙，給予氧療、藥物治療，共治療20周。

Ti=ti/si(i=1，2，…，N)

(2)

式中：Ti為第i個分析窗口下的單位地勢度；ti為平均起伏度；si為分析窗口面積。

(2)對數列T取對數lnT得到非線性數列樣本Xk(k=1，2，…，N)。

(4)計算統計量，公式為

(3)

(4)

式中：S為總樣本的離差平方和；Sk為前后兩段樣本的離差平方和之差。

變點的存在會使S和Sk的差距增大，S-Sk的最大值對應的分析窗口大小即為最佳分析窗口。

2.3 試驗方案

根據DEM特征，選取7個包含各類地貌的樣本:樣本①～③為矩形樣本，分別為山區小區、山區平原混合小區、平原小區；樣本④～⑥為京津冀地區依地形劃分的3個地貌大區，分別為高原大區、山區大區和平原大區；樣本⑦為京津冀全區(圖1)。對每個樣本進行分析窗口遞增的地形起伏度計算，窗口使用圓形鄰域。以往研究中最佳分析窗口一般不超過25 km2，本研究將最大窗口控制在35 km2以內，分析窗口設置見表1。起伏度提取過程在ArcGIS中進行，利用Modelbuilder將Focal statistics工具(用于統計最大值、最小值)、Raster Calculator工具(用于計算高差)和Band Collection Statistics工具(用于將起伏度的統計數據輸出)串聯為工作流實現起伏度快速提取，在Excel中利用VBA編程實現均值變點法中S-Sk值的計算。

表1 地形起伏度分析窗口的設置

注：圓形窗口計算面積時半徑柵格數需加0.5。

3 結果與分析

3.1 起伏度與分析窗口面積的統計特征

以每個分析窗口下計算的平均起伏度作為因變量，遞增的分析窗口作為自變量繪制散點圖(圖2)。

圖2 平均起伏度-分析窗口面積散點圖

如圖2所示，平均起伏度隨分析窗口面積的增大而增加，變化曲線呈對數函數或冪函數特征。基于這兩種函數進行回歸分析得到擬合方程(表2)。由表2知，對數函數和冪函數均有較好的擬合效果，R2均大于0.96。相比之下，使用冪函數的擬合效果優于對數函數，且受地貌的影響較小，例如當使用對數函數擬合時，不同樣本的R2在0.964 8～0.985 6之間變化，而使用冪函數擬合時，R2值更高且更集中，變化范圍為0.983 8～0.996 9。當分析窗口面積一定時，不同樣本的平均起伏度存在差異，這種差異在分析窗口面積較小時并不明顯，隨著分析窗口面積遞增，差異逐漸放大，例如分析窗口從0.16 km2增加到33.88 km2時，7個樣本平均起伏度的最大差值由85.54 m增加到578.43 m。

表2 平均起伏度與分析窗口面積的對數和冪函數擬合方程

3.2 最佳分析窗口

圖2中的每條曲線存在一個增大速度由快變慢的點(非數學拐點)，該點對應的窗口即為最佳分析窗口。根據公式(3)和(4)用均值變點法計算得到各樣本的系列S-Sk值，其與分析窗口面積的對應關系見圖3。圖3顯示，各樣本S-Sk值均呈單波峰曲線形態，峰值出現的位置較接近，說明各樣本的最佳分析窗口面積差別不大。但是，統計分析顯示7個樣本仍表現出分區性：樣本③和⑥屬一個區，其他5個樣本屬另一個區，兩區在S-Sk的均值、標準差、偏度、峰度等多個指標上均有明顯差異，這種分區性同樣體現在最佳分析窗口面積上，即樣本③和⑥的最佳分析窗口面積為4.64km2，其他5個樣本為5.35km2(表3)。

圖3 S-Sk與分析窗口面積散點圖

樣本均值標準差偏度峰度最佳分析窗口面積(km2)①山區小區16.476.47-0.59-0.945.35②山區平原混合小區15.426.08-0.58-0.955.35④高原小區14.425.65-0.59-0.945.35⑤山區大區15.796.18-0.60-0.925.35⑦京津冀全區15.576.10-0.59-0.935.35③平原小區21.238.46-0.56-0.984.64⑥平原大區18.517.41-0.54-1.004.64

3.3 地形地貌對最佳分析窗口的影響

顯然，地形差異是出現2種最佳分析窗口的原因。張偉等[9]曾基于SRTMDEM提取我國平原、高原、丘陵、山地4種地貌的最佳分析窗口，結果分別為4.83、4.97、4.79、4.83km2，即山地的分析窗口小于等于高原和平原。此結論與常理恰好相反，也就是說，最佳分析窗口的大小并不能直觀地衡量一個地區的宏觀地貌特征。事實上，根據前人對最佳分析窗口需滿足的山體完整性和區域普適性原則[1]的闡述，最佳分析窗口的大小取決于一個區域內的最大山體，而與小起伏度地區無關，即便后者占據更大面積。因此，張偉等的研究中平原和高原樣本很可能存在局部大起伏山體。

為了量化最佳分析窗口對起伏度的表達，以每個樣本在其最佳分析窗口下的最大起伏度Rm作為該分析窗口能夠表達的最大高差。7個樣本的Rm值見表4，其中：樣本③平原小區最小，為69m；樣本⑥雖同為平原區，但檢查DEM數據發現該樣本在北京順義區大孫各莊鎮北側發育小面積孤山，因而其Rm值遠大于樣本③，達到392m；高原樣本④Rm值為445m；樣本①、②、⑤、⑦均覆蓋一定山區，因而其Rm為高值；樣本⑦京津冀全區包含了區內所有的山體，其Rm值必然是最大的。

表4 各樣本的Rm值

由于DEM和分析窗口增大的步長均為分辨率較粗的離散數據，最佳分析窗口面積并沒有隨Rm單調增大，而是呈階梯狀增大(圖4)，這種階梯特征曾在涂漢明等[1]的研究中被證實。考慮到7個樣本中，樣本③平原小區起伏度非常小，不可能存在Rm值大于樣本⑦的地貌區，因此在SRTMDEM數據下京津冀地區僅存在4.64和5.35km2兩個最佳分析窗口(這2個窗口半徑柵格數為連續值13和14)。此外，樣本⑥和樣本④位于階梯抬升處且兩者Rm值較接近，因此可以大致將Rm值400m作為選擇2種最佳分析窗口的臨界值，即4.64km2窗口可以表達400m以內的高差，5.35km2窗口更適合于表達400m以上的高差。

圖4 最佳分析窗口面積與Rm的關系

3.4 京津冀地區的地形起伏度

在上文分析的基礎上，采用4.64和5.35km2雙窗口方案提取京津冀地區的地形起伏度。首先將4.64km2窗口下提取的起伏度中高差小于400m的部分保留，而將高差在400m以上的區域采用5.35km2分析窗口重新提取起伏度，將兩部分合并后，依據《中國1∶1 000 000地貌圖制圖規范》得到京津冀地區地形起伏度分級圖，并與5.35km2單值窗口提取方案進行對比(圖5)。圖5顯示，雙窗口方案和單窗口方案整體結果較接近，但從細節看，單窗口方案會將平原區起伏度夸大，例如在河北燕郊鎮周邊的平坦地區，單窗口方案識別出7個起伏度大于30m的微起伏柵格，在現有的多個地貌形態劃分方案下，這些柵格均被劃分為丘陵地貌[4]，這是不符合實際情況的，相比之下，雙窗口方案僅識別出2個起伏度大于30m的柵格(圖6)。在整個研究區內，兩種方案得到的平坦區域(起伏度小于30m)面積總計相差1 126km2。可見，雙窗口方案兼顧了京津冀地區半山半平原的地勢特點，因而更具有優勢。

(a)使用單窗口提取 (b)使用雙窗口提取

(a)使用單窗口提取 (b)使用雙窗口提取

京津冀地區的地形起伏度在0～1 145m之間，以平坦區、中起伏區和小起伏區為主，其中平坦區占比為43.99%，中起伏區和小起伏區占比分別為30.46%和14.12%，大起伏區占比為5.02%，極大起伏區極少(表5)。宏觀上各種地貌大致以東北—西南向展布，而在西北—東南向則表現出明顯的分帶性，即中間為中、大起伏山區，東南部為平坦的華北平原，西北部為微、小起伏和平坦的壩上高原。

4 結 論

表5 京津冀地區各級地形起伏度面積

京津冀地區幅員遼闊、 地貌類型復雜， 全局單值最佳分析窗口在該地區地形起伏度的提取中具有局限性。本研究計算了研究區內不同地貌樣本的最佳分析窗口，分析了地形地貌對最佳分析窗口的影響，在此基礎上提取了京津冀地區的地形起伏度。

京津冀地區地形起伏度隨分析窗口面積增大而增大，對兩者對應關系的擬合效果冪函數優于對數函數，前者各樣本R2均大于0.98。區內最佳分析窗口隨地形高差的增大而呈階梯狀增大，存在4.64和5.35km2兩個最佳分析窗口，前者可以表達400m以內的高差，后者更適合于表達400m以上的高差。利用4.64和5.35km2雙窗口方案提取的地形起伏度優于單窗口方案，兩種方案在中、大起伏區結果一致，但在平坦地區前者改善了后者對起伏度的夸大，更好地兼顧了研究區半山半平原的地形特征。分析結果顯示，京津冀地區的地形起伏度在0～1 145m之間，以平坦、中起伏和小起伏為主。