張 競(jìng)，杜 東，白耀楠，紀(jì)冬麗

(1.中國(guó)地質(zhì)調(diào)查局 天津地質(zhì)調(diào)查中心，天津 300170； 2.天津城建大學(xué) 市政與環(huán)境工程學(xué)院，天津 300384)

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

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

1 研究區(qū)概況

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

2 研究數(shù)據(jù)與研究方法

2.1 研究數(shù)據(jù)

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

圖1 京津冀地區(qū)DEM及研究樣區(qū)

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

地形起伏度一般基于DEM數(shù)據(jù)采用窗口分析法求取，其表達(dá)式為

Ra=Hmax-Hmin

(1)

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

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

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

1.2.1 對(duì)照組 對(duì)患者每個(gè)月開展1次肺健康知識(shí)宣教，內(nèi)容包括COPD相關(guān)知識(shí)、肺康復(fù)理念，指導(dǎo)戒煙，給予氧療、藥物治療，共治療20周。

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

(2)

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

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

(4)計(jì)算統(tǒng)計(jì)量，公式為

(3)

(4)

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

變點(diǎn)的存在會(huì)使S和Sk的差距增大，S-Sk的最大值對(duì)應(yīng)的分析窗口大小即為最佳分析窗口。

2.3 試驗(yàn)方案

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

表1 地形起伏度分析窗口的設(shè)置

注：圓形窗口計(jì)算面積時(shí)半徑柵格數(shù)需加0.5。

3 結(jié)果與分析

3.1 起伏度與分析窗口面積的統(tǒng)計(jì)特征

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

圖2 平均起伏度-分析窗口面積散點(diǎn)圖

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

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

3.2 最佳分析窗口

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

圖3 S-Sk與分析窗口面積散點(diǎn)圖

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

3.3 地形地貌對(duì)最佳分析窗口的影響

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

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

表4 各樣本的Rm值

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

圖4 最佳分析窗口面積與Rm的關(guān)系

3.4 京津冀地區(qū)的地形起伏度

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

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

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

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

4 結(jié) 論

表5 京津冀地區(qū)各級(jí)地形起伏度面積

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

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