基于多源數據的山前水系提取實驗:以濟南市及鄰區為例

王澤利, 許洪泰, 肖澤后, 倪永進

(1.山東科技大學地球科學與工程學院, 青島 266400; 2.山東省地震工程研究院, 濟南 250000)

作為構造地貌的重要標志,河流已成為新構造運動研究中的一種主要研究對象[1-4]。河流的流量、彎曲度、分級、流域面積、縱剖面形態等特征對區域的構造隆升、斷層活動等控制因素具有敏感的反應[5-6]。因此,快速、準確地提取河流信息,并通過對提取獲取的水系形態特征與斷裂構造之間的關系進行分析,來討論其對構造運動的響應,進而為研究區域構造穩定性提供地貌上的證據至關重要。

隨著空間數據采集方法的進步,基于數字高程模型(digital elevation model,DEM)的河流信息提取已在造山帶地區受到廣泛的研究與應用。但在沉積平原區,由于地貌平緩和人為影響干預較大的原因,該類方法應用受到極大限制。文獻[7-16]分別借助雷達影像和光學衛星影像開展了對水系的提取的研究工作,形成了如閾值法、濾波法、分類器法等提取方法,但以上多關注于小面積范圍的水系信息提取。值得指出的是,上述方法在各自針對的單一地質地貌區內都可以達成可接受的方法,但針對包含山地和平原的多相復雜地貌,均難以獨立實現良好的水系提取效果。

濟南市及鄰區位于魯西隆起區與濟陽凹陷的邊界部位,水系和斷裂構造均較發育。擬針對該區域,選用DEM數據和遙感光學影像數據,嘗試進行多源數據河流提取實驗,以期探索山前平原及山間地帶河流自動提取的方法。并通過分析提取出的水系在斷裂附近的形態特征,探討河流對構造活動的響應。

1 研究區概況

如圖1所示,選取山東省中部的濟南市及鄰區作為研究區,范圍為36°N～37.5°N,116.5°E～118°E。濟南及鄰區地處魯中南低山丘陵與魯西北沖積平原交界地帶,構成了山前地帶河流提取實驗的良好實驗地區。區內南部為泰山山脈,徂徠山山脈及兩者所圍限的萊蕪盆地,有牟汶河、贏汶河等大汶河主要支流流經。北部多為平原地區,區內分布北東向的黃河,徒駭河等河流。如圖2、圖3所示,DEM高程圖與坡度圖表現出相應特征,區內海拔自北西向南東遞增,最低點僅3 m,最高點為泰山山脈主峰玉皇頂,海拔1 545 m,與山前平原相對高差1 500 m以上,高程圖(圖2)顯示的1 517 m應是由于DEM數據分辨率所限而導致的誤差;坡度增加趨勢與海拔相同,研究區北部基本上都是較低的區域,高坡度地區大都集中在研究區的南部。

圖2 高程圖Fig.2 Elevation map

圖3 坡度圖Fig.3 Slope map

2 數據來源與研究方法

2.1 數據來源

結合數字高程數據與遙感影像數據,嘗試提取山前地帶的河流水系。數字高程模型最早于20世紀50年代提出,在地質災害,土地利用,景觀格局等方面具有重大作用[17]。本次實驗所用DEM數據于地理空間數據云中下載,名稱為SRTMDEM(shuttle radar topography mission digital elevation model) 90 m分辨率原始高程數據,中心經度為117.5°,中心維度為37.5°,條帶號60,行編號5,數據標識為srtm_60_05。

Landsat8衛星是美國陸地衛星計劃(Landsat Program)的第8顆衛星,由美國國家航空航天局(National Aeronautics and Space Administration,NASA)與美國地質調查局(United States Geological Survey,USGS)合作開發、軌道科學公司建造。該衛星于2013年2月11號在加利福尼亞范登堡空軍基地由Atlas-V火箭搭載發射成功,攜帶陸地成像儀(operational land imager,OLI)和熱紅外傳感器(thermal infrared sensor,TIRS)。OLI陸地成像儀包括9個波段,空間分辨率為30 m,其中包括一個15 m的全色波段,成像寬幅為185 km×185 km;熱紅外傳感器TIRS包括2個單獨的熱紅外波段,分辨率100 m[18]。本次實驗所用的Landsat8亦是于地理空間數據云中下載,云量低于5%,成像時間為2019年9月27日,條帶號均為122,行編號分別為34、35,數據標識分別為LC81220342019270LGN00、LC81220352019270LGN00。

2.2 提取方法

理想情況下,地表水流向與地形坡向一致,并傾向于沿局部最低海拔和最大坡度形成地表徑流,基于這一認識,發展出利用DEM提取線狀河流的方法。在ArcGIS10.2平臺的水文分析模塊中,對下載的DEM數據進行鑲嵌、填洼等預處理,然后計算流向、流量,確定流量閾值分析河網,對河流進行分級和矢量化工作。

流向的確定是采用的D8算法,8個數字分別代表8個方向,以中心柵格為參照,計算與周圍柵格的距離權落差,將結果最大的設定為流出位置,生成流向柵格數據。

流量閾值的確定工作是影響河網提取結果的關鍵因素,本實驗采用河網密度法[19-22]來確定流量的具體閾值。河網密度指單位面積流域上的河流總長度,其計算公式為

(1)

式(1)中:D為河網密度;L為對應閾值下的河流總長度;S為研究區域面積。

利用ArcGIS生成不同閾值下的河網,統計河流長度并計算出河網密度;再擬合函數關系,繪制河網密度與閾值的圖像;最后根據拐點位置找出河網密度隨閾值變化趨于平緩的坐標,以該點的閾值進行河網提取可得到較好的準確度。依據此方法結果閾值提取河網并將其分級、矢量化后即可形成依據DEM數據提取的最終河流數據。

不同光學性質的地物對不同波段光線的反射率存在差異,構成了基于多波段遙感影像水體提取的理論基礎,當前主要的方法包括歸一化差異水體指數[9](normalized difference water index,NDWI)、改進的歸一化差異水體指數[10](modified normalized difference water index,MNDWI)、Gauss歸一化差異水體指數[23](Gauss normalized difference water index,GNDWI)、增強水體指數[24](enhanced water index,EWI)等,各自的特征和適用性具體如下。

(1)NDWI法。水體反射率從綠光到近紅外波段逐漸遞減,在近紅外波段最弱;而植被的反射率則相反,在近紅外波段反射率最強[9]。因此,本方法采用綠光和近紅外波段構建模型,突出顯示水體的同時抑制了植被等地物,特點是能夠抑制植被信息,突出水體,適用于地勢開闊平坦地區。

(2)MNDWI法。該方法以短波紅外波段(short-wave infrared,SWIR)替換近紅外波段(near infrared,NIR),顯著增大了水體與建筑物的反差。該方法特點是能夠較好地去除居民地和土壤的影響,有利于提取城市水體。

(3)GNDWI法。針對線狀河流水體的精確提取,可較好保留河流完整性,受陰影影響較大,需要用到高程信息,實現過程比較麻煩。

(4)EWI法。能夠抑制居民地、土壤和植被噪聲,易受到陰影及淺灘的影響,適合半干旱地區水體提取。

考慮研究區面積較大,覆蓋了山野、村落和中大型城市,既有陡峭山地,又包括開闊平原,選取NDWI法和MNDWI法。

NDWI的數學模型可表示為

(2)

MNDWI的數學模型可表示為

(3)

式中:Green為綠光波段反射率,對應Landsat8數據的第3波段;NIR為近紅外波段反射率,對應Landsat8數據的第5波段;SWIR為短波紅外波段反射率,對應landsat8數據的第6波段。

遙感影像處理在ENVI5.3平臺上進行,在完成輻射定標、大氣校正、圖像鑲嵌等預處理后,依據式(2)、式(3)分別計算獲得NDWI和MNDWI圖像。

3 提取結果與分析

3.1 DEM提取結果

經過預處理、填洼及流向分析后,分別選取100、200、400、600、…、3 600共19個集水流量閾值對研究區進行水文信息提取,計算出河網密度,如表1所示。繪制出閾值與河網密度的散點圖,利用Origin對其進行曲線擬合,結果如圖4所示,擬合后確定冪函數R2=0.999 7,此數值越接近一說明擬合效果越好,可以看出,河網密度開始隨著閾值的增大而急劇減小,后下降趨勢趨于平緩,對擬合曲線求二階導數可知,閾值在約1 000時,其數值幾乎接近于零,并且隨著閾值增大其變化趨勢也不明顯,因此此法求得的適宜閾值為1 000,并以此結果作為流量閾值進行水文信息提取,然后完成河流分級及河網矢量化工作,即可獲得區內基于DEM提取出的矢量化水系信息圖,如圖5所示。

表1 河網密度計算結果

圖4 擬合曲線Fig.4 Fitting curves

圖5 DEM提取圖Fig.5 DEM extraction diagram

3.2 遙感數據提取結果

利用ENVI對遙感影像數據完成預處理后,按照上述比值模型公式,分別計算出NDWI、MNDWI的結果,如圖6、圖7所示,去除異常值后兩圖像的像元亮度值(digital number,DN)均在-1～1范圍內,經多次試驗調整后,確定NDWI閾值為0,MNDWI閾值為0.4,DN值大于閾值的像素判定為水體。

（2）經產母豬情期受胎率：采用人工授精為86.93%，采用自然交配為78.69%，差異顯著（P<0.05）。

圖7 MNDWI提取圖Fig.7 MNDWI extraction diagram

從NDWI與MNDWI兩者的結果來看,研究區內黃河的提取效果兩者幾乎相同,但對于徒駭河這種次級水系MNDWI提取結果的連續程度、完整程度要好于NDWI,此外對濟南市區和山區的噪音處理方面MNDWI表現較為突出,結合上述兩者的表現情況,認為MNDWI相較于NDWI體現出了更加豐富的水文信息,即MNDWI更加適用于研究區,所以將會用MNDWI方法的提取結果來參與接下來的研究步驟。

3.3 適用地貌判定

比較基于DEM方法與MNDWI方法提取的水文信息,發現無論是何種方法,提取水文信息的準確程度在研究區中的不同區域是不同的,并且提取水文信息準確程度的高低與地貌類別具有極高的相關性,為了判定不同方法對不同地貌的適用程度,在研究區內選取了兩處不同地貌的典型區域平原(M1)、山地(M2),如圖8所示,并通過目視解譯提取出M1、M2中的水文信息,通過兩種方法分別在M1、M2中的提取結果與目視解譯結果相對比來分析適用性,進而判定適用地貌,其準確率結果如表2所示。

表2 提取準確率比較

圖8 不同地貌提取結果對比Fig.8 Comparison of extraction results of different landforms

如圖8、表2所示,在平原(M1)中MNDWI方法效果較好,準確率可達90%以上,雖然由于遙感影像分辨率、人為建筑等因素的影響,有些細小水系沒有提取出,但是沒有像DEM方法這樣出現大量與實際不符的偽河道,山區(M2)中DEM方法效果較好,準確率也可達90%以上,河流形態符合實際,而MNDWI方法只提取出了面狀水系,故此準確率只有27%。這兩種方法存在的局限性,導致它們均難以獨自實現對此地貌的水系自動提取。而兩種提取水文信息的方法在同一地貌表現出的互補特性,為將兩者以地貌為依據進行綜合應用提供了可能。

3.4 山前水系的綜合提取

依據前述結果,在山地丘陵地貌區,利用DEM數據提取河流的準確率顯著高于遙感影像自動提取,而在平原地區結果相反。因此在不同地貌區利用最優方法提取后再行合成可達到自動識別的最佳效果。為達成此目標,如何正確劃分地貌區至關重要。

利用地形起伏度區分山地和平原是近十年來被廣泛探討的方法[25-26]。地形起伏度系指特定范圍內的最大高度差,是描述區域地形特征的宏觀指標,其特點在于具有尺度效應,即通過調節計算窗口的大小,可適用于不同面積的工作范圍,利用均值變點法[27]來確定適合本區的計算窗口大小。

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

(4)

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

(3)對每個k(k≥ 2),將樣本數列分為前后兩段,即X1,X2,…,Xk-1和Xk,Xk+1,…,XN,分別計算前后兩段數列的算術平均值Xk1、Xk2及總樣本的算數平均值X。

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

(5)

式(5)中:Xt為某個具體分析窗口單位地勢度T的對數值。

(6)

式(6)中:Sk為前后兩段樣本的離差平方和之和。

S為總樣本的離差平方和,變點的存在會使S與Sk之間的差距增大,S-Sk的最大值對應的分析窗口,即為最佳分析窗口。

ArcGIS中“焦點分析”模塊下計算并統計2×2～32×32,面積為3.24×104～829.44×104m2共31個遞增矩形窗口,統計數據如表3所示,繪制出S-Sk隨窗口大小變化的散點圖并對其進行非線性曲線擬合,在窗口大小為11×11時,S-Sk具有最大值16.85,即適宜計算地形起伏度的窗口大小為0.98 km2。

表3 適宜計算窗口統計結果

選取大小為11×11網格計算出的區內平均起伏度數據,以30 m為界線將研究區地貌劃分為平原與山地丘陵區,起伏度小于30 m的地區圈劃為平原,否則圈劃為山區[28],在山區應用DEM并保留MNDWI提取的湖泊等面狀水系,在平原地區應用MNDWI水體指數,將兩者基于ArcGIS10.2進行空間綜合修正合并,得到圖9所示的研究區內綜合水系圖。

從合成后的結果(圖9)來看,山地丘陵區水系提取較完整,湖泊與河流空間位置吻合良好;平原地區受遙感影像分辨率、地表植被、建筑物和人工覆蓋的影響,流量較大的河流(如黃河等)可以被近完整的提取,但流量較小的河流斷續存在。兩類地貌接觸地帶展現了類似的規律,大規模河流(如玉符河、北大沙河等)基本可以完整的連接,說明無論基于DEM亦或遙感影像自動提取,均正確識別了該河流。更多流量較小的河流被截斷在山體丘陵區邊界,即依據遙感影像自動提取未能識別該河流,為核實該類河流是否真實存在,對代表性區域進行了人工解譯,結果發現該類河流真實存在,但由于是季節性河流或者被周邊植被遮擋等原因,導致沒有被遙感影像自動提取。全圖未出現終止在山地丘陵區邊界的平原河流,可以判斷在轉換帶附近基于DEM數據提取方法獲得的數據更豐富。值得關注的是,在轉換帶附近人工建筑對提取結果的干擾在遙感影像自動提取方法上更加明顯。如濟南市區,該市地表徑流豐富,甚至形成了大明湖等較大面積的湖泊,但受人工建筑干擾的影響,流經該區的河流并未完整銜接。

3.5 新構造活動響應

當代構造地貌學研究普遍認識到,隨著構造運動的進行,河流會持續記錄地貌演化過程,因此對河流地貌分析可反映出長時間尺度的構造活動信息[29-31]。在均衡狀態下,河流的空間形態和縱剖面表現相對平滑;當河流流經的斷層兩盤發生足以影響地表形態的相對運動時,這一均衡狀態被打破,進入一種瞬時不均衡狀態,即在縱剖面上形成裂點或在地表形成流向變化,這種瞬態地貌會隨著時間發生空間遷移,如裂點的溯源侵蝕[32]和彎曲度的變化[33-34],并逐漸消失已達到新的均衡狀態。

研究區南部為泰山山脈、徂徠山山脈及由二者圍限的萊蕪盆地。泰山基底為古老的太古宇泰山巖群,受太古宙-元古宙多次巖漿活動、褶皺構造的影響,變質變形作用明顯,主要由斜長角閃巖組成,其次為黑云變粒巖、角閃變粒巖等,其原巖為超基性、基性火山巖和火山凝灰巖建造,后因泰山大幅抬升,部分蓋層風化剝蝕,使得基底重新出露;蓋層為古生界寒武系-奧陶系灰巖、頁巖,主要分布于北側張夏崮山一帶,從老到新出露饅頭組、毛莊組、徐莊組、張夏組、崮山組、長山組、鳳山組,與基底呈角度不整合接觸。受泰山山前斷裂等限制,所述山脈普遍表現為南斷北超的單斜斷塊山系[35-36]。

研究區北部為大面積的沖洪積平原。黃河、徒駭河由該部分北東向經過,經黃河等沖、洪積形成了巨厚的第四系沉積,可達60～310 m,自底向頂可分為平原組、黑土湖組、寒亭組、黃河組等巖石地層。該部分地勢平緩,起伏較小,最低海拔約為5 m。

研究區以廣饒-齊河斷裂為界,北部為華北東部盆地斷坳區,區內斷裂以近東西向為主,如臨邑斷裂、興隆斷裂等,第四紀以來不活動或僅在早、中更新世發生過活動;南部魯西斷塊隆起區內北西向萊蕪斷裂、泰山西麓斷裂、長清斷裂等也多為早、中更新世活動,但北西向的銅冶店-孫祖斷裂和近東西向的泰山山前斷裂等被認為是晚更新世活動斷裂[37-39]。

將區內的地層與第四紀斷裂資料完成矢量化后疊加到綜合水系圖中,得到圖10所示的綜合水系地質圖,通過對比研究區的斷層和水系分布,分析河流對活動構造的響應:從水系與斷裂水平空間分布結果來看,黃河在長清斷裂附近、徒駭河在桑梓店斷裂附近流向發生了顯著變化,但這種變化與走滑斷裂引起的河流變形現象存在較明顯差異,后者多表現為短距離的快速曲折,斷層影響帶兩側流向整體保持一致。黃河在廣饒-齊河斷裂、千佛山斷裂附近,徒駭河在夏口斷裂、林南斷裂附近、馬頰河在寧南斷裂附近、玉符河在千佛山斷裂附近等,水系未發生明顯的流向、彎曲程度變化。泰山山前斷裂、萊蕪斷裂本身位于兩種地貌區的邊界,因此有石汶河、贏汶河河流等終止于該斷裂,此種河流變化的起因更多的是基于DEM和遙感影像提取水系效果的差異。但斷裂與地貌分區邊界吻合良好這一現象本身也可說明其對地貌可能存在較大的控制作用。

圖10 綜合水系地質圖Fig.10 Geological map of integrated drainage system

在斷裂與河道交匯處利用GIS軟件提取出斷裂兩側河道的高程數據,發現斷裂兩側河道高程變化不明顯,表明河道的現在垂直形態沒有表現出受斷層影響后的特征。

從總體的河流水平及垂直形態的變化來看,研究區范圍內識別的河流目前未有受斷層影響的顯著表現,至少說明研究區作為新構造運動的弱活動區,斷層走滑及垂直運動造成的地貌影響已充分衰減,河流已經達到了平衡狀態。

4 結論

以濟南及鄰區作為研究區,以SRTMDEM 90 m分辨率原始高程數據和Landsat 8 OLI_TIRS衛星影像作為源數據,通過對比不同方法(DEM、NDWI、MNDWI)的特點及在各地貌區的應用效果,嘗試了一種能夠自動提取上述地區水系的方法。結果顯示,在區內的山地丘陵區和平原區,所對應的最優水系提取方法分別為基于DEM和MNDWI的提取方法,兩分區的分界線可通過地形起伏度進行確定。通過對比研究區的斷層和水系分布,認為研究區內的早中更新世斷裂對河流造成的影響已經衰減,并且河流已達到了新的平衡狀態。