渭河干流楊凌段水邊線提取及其分維數(shù)計(jì)算方法研究

夏積德， 江仕榮， 周 波

(1.楊凌職業(yè)技術(shù)學(xué)院， 陜西 楊凌 712100； 2.西北農(nóng)林科技大學(xué)， 陜西 楊凌 712100)

1 研究背景

分形理論是對復(fù)雜形體不規(guī)則性的量度，以探索和揭示隱藏在復(fù)雜現(xiàn)象不規(guī)則形態(tài)背后的規(guī)律及本質(zhì)聯(lián)系[1]，并可用分形維數(shù)來反映復(fù)雜形體占有空間的有效性。水邊線是復(fù)雜形體，通過對水邊線盒分形維數(shù)的計(jì)算有助于了解河流的侵蝕發(fā)育階段，為水土流失防治提供參考和依據(jù)。由于多光譜遙感影像具有信息量豐富的特點(diǎn)，目前在大多數(shù)研究中均采用歸一化差值水體指數(shù)(NDWI)方法對多光譜遙感影像數(shù)據(jù)進(jìn)行處理以提取水邊線信息[2-3]。在實(shí)際操作過程中，多光譜遙感影像又具有采集周期長、占用儲存空間大、計(jì)算要求高等缺點(diǎn)。相比之下，RGB遙感影像(如GOOGLE影像、BAIDU影像)更容易獲取，但目前基于RGB遙感影像提取水邊線的研究較少。

在采集水邊線信息時(shí)，多采用影像波段選取法、人工交互閾值分割法、Roberts算子法、Sobel算子法、Canny算子法、Laplace算子法、面向?qū)ο蠓诸惙ǖ萚4]。其中，Sobel算子法是根據(jù)灰度圖邊緣表現(xiàn)不連續(xù)的規(guī)律來提取邊緣信息，多用于小范圍的檢測[5-8]；人工交互閾值分割法是通過觀察圖像的灰度直方圖確定閾值，根據(jù)計(jì)算結(jié)果對閾值進(jìn)行不斷修正的水邊線提取方法；面向?qū)ο蠓诸惙▽⒕哂羞B續(xù)色調(diào)的單色影像按一定密度范圍分割成若干等級，然后根據(jù)等級再將地物分割[9]。這3種方法更為常用。

為探究在缺少多光譜影像時(shí)能否運(yùn)用Sobel算子法、人工交互閾值分割法、面向?qū)ο蠓诸惙▽λ吘€信息進(jìn)行提取，本文以渭河干流楊凌段為例進(jìn)行研究，并與DEM數(shù)據(jù)結(jié)果進(jìn)行了對比分析，以期為研究區(qū)水沙運(yùn)移和河道管理提供依據(jù)。

2 研究區(qū)概況與數(shù)據(jù)來源

2.1 研究區(qū)概況

渭河是黃河最大的一級支流 ，發(fā)源于甘肅省渭源縣鳥鼠山，自西向東流經(jīng)甘肅、寧夏、陜西，于陜西潼關(guān)注入黃河，全長818 km，流域總面積13.5×104km2。渭河流域楊凌段位于渭河中游段，區(qū)間無較大支流匯入，流域面積38 017 km2，占渭河流域總面積的28.2%，渭河自西向東從楊凌示范區(qū)南側(cè)穿過，全長約11.8 km，河床寬500～1 000 m，多年平均流量約為136.5 m3/s，年總徑流量約為43.06×108m3。楊凌位于陜西黃土高原南部、關(guān)中平原西部，是我國唯一的農(nóng)業(yè)高新技術(shù)產(chǎn)業(yè)示范區(qū)，也是國家批準(zhǔn)為向亞太經(jīng)合組織開放的十大工業(yè)園區(qū)和國家重點(diǎn)支持的五大高新區(qū)之一。區(qū)內(nèi)地勢整體趨于平緩，北高南低，西高東低。東臨武功，西倚扶風(fēng)，南面渭河，北部為塬地，總面積135 km2。楊凌全區(qū)屬于溫帶半濕潤半干旱氣候，降水多為夏季和秋季，年平均降水量約630 mm，年均氣溫為13 ℃，夏季炎熱，冬季寒冷，全區(qū)塬、坡、灘地交錯，土壤肥沃，適宜多種農(nóng)作物生長。

2.2 數(shù)據(jù)來源

DEM數(shù)據(jù)分別為地理空間數(shù)據(jù)云網(wǎng)站下載的2009年ASTER GDEM數(shù)據(jù)(30m×30 m分辨率)和楊凌區(qū)2009年1∶10000地形圖數(shù)字化后經(jīng)過普通克里金插值生成(5 m×5 m分辨率)。

RGB影像數(shù)據(jù)來源于高分單波段Google Earth(以下簡稱GE數(shù)據(jù))。其中，2009、2010、2014年GE數(shù)據(jù)空間分辨率分別為0.536 m×0.536 m，2011年GE數(shù)據(jù)空間分辨率為1.072×1.072 m。

3 研究方法

3.1 Sobel算子法

Sobel算子法基本原理[10-12]：假設(shè)連續(xù)圖像函數(shù)為f(x,y),函數(shù)在(x,y)處的梯度是一個(gè)具有方向和大小的矢量，即：

(1)

式中：i、j分別為x、y方向的單位矢量。

梯度算子的幅值和方向分別為：

(2)

(3)

式中： grad[f(x,y)]值即為圖像在(x,y)處的邊緣數(shù)據(jù)，將邊緣數(shù)據(jù)與事先設(shè)置好的閾值對應(yīng)，大于閾值則輸出為邊緣像素，否則輸出為非邊緣像素。

3.2 人工交互閾值分割法

閾值分割法[13]是根據(jù)圖像的灰度直方圖，按照一定的規(guī)則選擇兩個(gè)峰值之間的谷底作為閾值，對圖像進(jìn)行檢測的方法。其表達(dá)式為：

(4)

式中：g(x)為運(yùn)算后的二值圖；T為閾值。

3.3 面向?qū)ο蠓诸惙?/h3>

面向?qū)ο蠓诸惙╗14]以建立統(tǒng)計(jì)識別函數(shù)為理論基礎(chǔ)，依據(jù)典型樣本訓(xùn)練方法進(jìn)行分類的技術(shù)，根據(jù)預(yù)先給定的訓(xùn)練樣本，得到各類柵格的判別函數(shù)g1,g2, …，gc。根據(jù)判別函數(shù)計(jì)算各柵格X的值g1(x),g2(x),…，gc(x)，選擇其中的最大值，對X進(jìn)行分類。

3.4 DEM河網(wǎng)法

根據(jù)流域?qū)嶋H對DEM數(shù)據(jù)進(jìn)行填洼處理，計(jì)算出各中心柵格與周圍柵格之間的距離權(quán)落差，選取差值最大的方向作為水流流向。根據(jù)柵格水流流向計(jì)算出各個(gè)點(diǎn)的流量值，結(jié)合各柵格的水位單元體積，得到區(qū)域內(nèi)的匯流累計(jì)量水系河網(wǎng)。距離權(quán)落差計(jì)算公式為：

S=ΔZ/D

(5)

式中:ΔZ為中心柵格與周邊柵格之間的高程差，m；D為柵格中心之間的距離,m。

3.5 盒分維計(jì)算

通過不同邊長δ的柵格所覆蓋的水系圖，計(jì)算相應(yīng)的覆蓋水系圖柵格數(shù)量N(δ)。根據(jù)分形理論，二者關(guān)系式為：

δ∝N(δ)

(6)

由正比關(guān)系可知，對關(guān)系式(6)的兩邊同時(shí)取以e為底的對數(shù)，也滿足正比關(guān)系。即：

ln (δ)∝ln (N(δ))

(7)

由柵格邊長δ的變化得到數(shù)集(δ1,δ2, …,δx-1,δx)和柵格數(shù)量集(N(δ1),N(δ2),… ,N(δx-1),N(δx))，同時(shí)取兩者的對數(shù)得到ln (δ)和ln (N(δ))數(shù)集，以ln (δ)為橫坐標(biāo)、ln (N(δ))為縱坐標(biāo)繪制散點(diǎn)圖，并擬合得趨勢線方程與判定系數(shù)R2。

ln (N(δ))=-Dln (δ)+B

(8)

式中：由趨勢線方程所得斜率D即為研究區(qū)地貌形態(tài)分形信息維數(shù)。

4 結(jié)果分析與討論

4.1 盒分維值計(jì)算與比較

以2009年GE數(shù)據(jù)和DEM為例，分別選取邊長δ為 0.1、0.2、0.4、0.8、2、4、8、10、16、32、64，得到對應(yīng)的N(δ)以及l(fā)n (δ)-ln (N(δ))關(guān)系圖，見圖1。

根據(jù)圖1(a)，各方法所得曲線最佳的邊長δ在(1，5)的范圍內(nèi)，其中采用RGB-GE數(shù)據(jù)的3種方法，其拐點(diǎn)對應(yīng)邊長均大于原始圖像分辨率，而DEM數(shù)據(jù)拐點(diǎn)對應(yīng)邊長均小于原始圖像分辨率。根據(jù)圖1(b)中的線性關(guān)系，確定各方法的盒分維值D及判定系數(shù)R2，見表1 。

圖1 不同柵格邊長δ條件下δ-N(δ)與ln (δ)-ln (N(δ))關(guān)系圖

表1 5種計(jì)算方法的盒分維表

根據(jù)馮小慶[15]對水系地貌侵蝕發(fā)育階段的劃分方法，當(dāng)水系分維數(shù)小于1.6時(shí)，流域地貌屬于侵蝕發(fā)育階段的幼年期；分維數(shù)值越趨近于1.6，河流地貌的侵蝕階段越趨近于幼年晚期，河流的下蝕作用減弱，旁蝕作用加強(qiáng)，地面分割的程度越來越碎。由表1可知，不同數(shù)據(jù)和方法所計(jì)算的盒分維值D∈[1.0035，1.3090]，D值大小排序?yàn)椋篋EM河網(wǎng)法<人工交互閾值分割法<面向?qū)ο蠓诸惙?Sobel算子法，且這4種方法的判定系數(shù)R2均大于0.99，這說明4種方法用來提取水邊線均具有良好的可行性。其中，DEM河網(wǎng)法的判定系數(shù)最高，接近于1。這說明與其他方法相比，DEM方法更能較好地反映河流的侵蝕發(fā)育階段，這和王林等[16]的研究結(jié)果一致。當(dāng)研究區(qū)數(shù)據(jù)采用RGB-GE影像數(shù)據(jù)時(shí)，人工交互閾值分割法所得分維值的判定系數(shù)較面向?qū)ο蠓诸惙ê蚐obel算子法略高，具有較高的提取精度。不同方法提取的渭河干流楊凌段水邊線的局部圖如圖2所示。

由圖2可以看出，3種方法所提取水邊線長度關(guān)系為：Sobel算子法>人工交互閾值分割法>面向?qū)ο蠓诸惙?。Sobel算子法提取水邊線信息時(shí)受噪聲影響提取了多余信息，使得水邊線長度增加(圖2(a))；面向?qū)ο蠓诸惙o法完整的提取水邊線信息(圖2(c))；人工交互閾值分割法提取水邊線的效果最佳(圖2(b)。

4.2 盒分維年際變化

根據(jù)上述分析結(jié)果，以2009年RGB-GE數(shù)據(jù)為例，基于人工交互閾值分割法的渭河干流楊凌段水邊線提取結(jié)果見圖3。

以不同年份楊凌區(qū)水邊線參數(shù)見表2。

圖2 基于RGB-GE數(shù)據(jù)的渭河干流楊凌段水邊線提取局部比較

圖3 基于人工交互閾值分割法的渭河干流楊凌段水邊線提取結(jié)果(以2009年為例)

表2 基于人工交互閾值分割法的不同年份GE數(shù)據(jù)水邊線參數(shù)

由圖3和表2可以看出，人工交互閾值分割法能夠完整地提取出水邊線信息，各年份水邊線盒分維值小于1.4，研究區(qū)水邊線侵蝕發(fā)育階段均處于幼年期，各判定系數(shù)R2均大于0.99；對比各年降水量數(shù)據(jù)，楊凌區(qū)水邊線長度與河網(wǎng)密度值具有良好的契合度。因此，人工交互閾值分割法具有良好的可行性。

4.3 討 論

目前可獲取的多光譜影像分辨率多為30 m，而0.5、1 m乃至更高分辨率的RGB影像更容易收集，且多光譜影像相較于RGB影像的數(shù)據(jù)量大，存儲要求高，處理難度也更高。本文應(yīng)用兩種數(shù)據(jù)進(jìn)行了實(shí)例比較，結(jié)果表明，除了多光譜影像數(shù)據(jù)和地形圖數(shù)據(jù)外，RGB影像(如GE)數(shù)據(jù)亦可用于水邊線信息的精確提取，本文的研究結(jié)果可為地區(qū)的水沙運(yùn)移和河道管理提供可靠的數(shù)據(jù)和方法依據(jù)。然而，在研究面積較大時(shí)RGB影像是否也能適用于其他參數(shù)值的計(jì)算，如何提高RGB影像數(shù)據(jù)的提取精度與效率等問題，仍有待進(jìn)一步系統(tǒng)研究。

5 結(jié) 論

本文以渭河干流楊凌段為研究區(qū)，采用Sobel算子法、人工交互閾值分割法、DEM河網(wǎng)法和面向?qū)ο蠓诸惙ǎ枚喙庾V影像、RGB-GE影像和地形圖數(shù)據(jù)進(jìn)行了水邊線信息提取和分維數(shù)計(jì)算。

(1)Sobel算子法、人工交互閾值分割法、面向?qū)ο蠓诸惙?、DEM河網(wǎng)法所計(jì)算得到盒分維方程的判定系數(shù)均大于0.99，說明這4種方法均能有效用于提取河流水邊線信息。

(2)Sobel算子法因?yàn)樵肼暤挠绊?，?dǎo)致盒維數(shù)偏大；人工交互閾值分割法需要人工參與，工作量較大，若大面積使用需要投入較多人力；面向?qū)ο蠓诸惙ǖ挠?jì)算量較大，計(jì)算時(shí)間長，對計(jì)算機(jī)要求較高；DEM河網(wǎng)法由于數(shù)據(jù)分辨率較低，對河網(wǎng)細(xì)節(jié)的反應(yīng)能力較弱，因此推測盒維數(shù)可能小于真實(shí)值。

(3)綜合來看，當(dāng)數(shù)據(jù)精度較高、研究面積較小時(shí)，人工交互閾值分割法優(yōu)于Sobel算子法和面向?qū)ο蠓诸惙?。渭河干流楊凌段水邊線分維值在1.0～1.4之間，河網(wǎng)地貌屬于侵蝕發(fā)育階段的幼年期，水邊線長度在[201 173.86 m, 544 210.45 m]，河網(wǎng)密度為0.080 m/m2，符合客觀實(shí)際，具有合理性。