摘要：水文分析實(shí)戰(zhàn)教程河流平均比降計(jì)算本章導(dǎo)讀河流比降指的是河流任意兩端點(diǎn)間的高程差與兩點(diǎn)間的水平距離之比值，簡(jiǎn)單的概括為單位河長(zhǎng)內(nèi)的落差，其單位一般都是。讀者可以參考水文分析實(shí)戰(zhàn)教程河流提取與河網(wǎng)分級(jí)的做法去提取。

本章導(dǎo)讀：河流比降指的是河流任意兩端點(diǎn)間的高程差與兩點(diǎn)間的水平距離之比值，簡(jiǎn)單的概括為單位河長(zhǎng)內(nèi)的落差，其單位一般都是‰ 。 河流（縱）比降對(duì)地質(zhì)災(zāi)害、水運(yùn)、水電、礦砂和新構(gòu)造運(yùn)動(dòng)研究有實(shí)用價(jià)值。比降是水文學(xué)上重要的一項(xiàng)指標(biāo)。在絕大部分的GIS通用平臺(tái)上，都沒(méi)有實(shí)現(xiàn)平均比降的計(jì)算，原因在于計(jì)算比降的條件比較苛刻，不符合通用GIS的要求。本章同樣是一個(gè)實(shí)操章節(jié)，筆者將會(huì)利用ArcGIS軟件，結(jié)合其ArcPy編寫(xiě)一個(gè)平均比降計(jì)算的工具，在通用GIS軟件上實(shí)現(xiàn)河流平均比降計(jì)算。 BY 李遠(yuǎn)祥

河流比降是很多水文研究者研究的內(nèi)容，其水文意義非常的大。這里研究的河流比降指的是河流的縱比降，是河流任意兩端點(diǎn)間的高程差與兩點(diǎn)間的水平距離之比值，即單位河長(zhǎng)的落差。

正是因?yàn)橹苯颖磉_(dá)的是單位河長(zhǎng)的落差，所以在水文學(xué)上其最直接的影響就是流速，帶來(lái)更多的水動(dòng)力引發(fā)的能量轉(zhuǎn)移，如勢(shì)能與動(dòng)能的轉(zhuǎn)化（涉及到水電站的建設(shè)）；河水對(duì)河床和河岸的沖刷，影響泥沙量；大面積持續(xù)降雨時(shí)，洪峰到來(lái)的時(shí)間和強(qiáng)度等。所以，比降在水文學(xué)上從來(lái)就不缺少話(huà)題。

河流的比降分為床面比降和水面比降。床面比降，用以表示河床縱斷面地形的變化；水面比降即河流中任意兩端點(diǎn)間的瞬時(shí)水面高程差與其相應(yīng)水平距離之比，用以表明河流全程或分段的水面坡度，故又稱(chēng)水力坡度，通常說(shuō)的河流比降就是河流水面比降。

在恒定均勻流情況下，水面坡度恰好等于河底坡度。所以，在沒(méi)有大量斷面數(shù)據(jù)支撐的情況下，根據(jù)地形和河流線型去研究的河流比降都是床面比降。

筆者在這里研究的都是床面比降，因?yàn)橹豢紤]地形與河流本身的幾何形狀，沒(méi)有任何的站點(diǎn)數(shù)據(jù)支撐。但在非洪水情況下，這個(gè)床面比降與水面比降基本上是吻合的。

在河流相距較短的兩個(gè)點(diǎn)上計(jì)算河流比降，意義不是十分大。一般需要計(jì)算的是整個(gè)河流或者某一河段的平均比降，這對(duì)于空間跨度較大的河流來(lái)說(shuō)才有研究?jī)r(jià)值。

本章研究的內(nèi)容也就是河流（段）的平均比降。研究河流平均比降的方法有幾種，比較流行的是【約翰斯通-克羅斯】法。其原理是根據(jù)比降大致均勻的原則，將河段劃分為若干個(gè)子河段，然后去計(jì)算整體的平均比降。其公式如下：

公式中個(gè)參數(shù)如下
S---- 河段平均河底比降
Li----第i個(gè)子河段長(zhǎng)度
Si----第i個(gè)子河段河底比降
N---- 子河段數(shù)量

另一種方法是【勞倫斯】法，其需要經(jīng)驗(yàn)系數(shù)和流速，相對(duì)來(lái)說(shuō)比較復(fù)雜，需要結(jié)合測(cè)站的一些數(shù)據(jù)。由于筆者沒(méi)有這些數(shù)據(jù)，就不作論述。

本章考慮的是在ArcGIS中實(shí)現(xiàn)【約翰斯通-克羅斯】法的計(jì)算。

利用GIS手段來(lái)研究水文，其核心方式是使用DEM數(shù)據(jù)去模擬水文過(guò)程。GIS的空間疊加能力非常的管用，非常簡(jiǎn)單的就能獲取到河流線對(duì)應(yīng)的DEM像元值，也就是河流對(duì)應(yīng)的高程值，這樣用于計(jì)算比降是非常有用的。最起碼解決了河流高程點(diǎn)獲取的問(wèn)題。

河流沿線高程的獲取，不能直接使用原始的DEM數(shù)據(jù)。原因有兩點(diǎn)：

如果河流數(shù)據(jù)與DEM數(shù)據(jù)不是同一份數(shù)據(jù)，那么像元點(diǎn)與河流線的重合程度有待商榷。

如果河流數(shù)據(jù)與DEM數(shù)據(jù)吻合（或者河流數(shù)據(jù)直接從DEM數(shù)據(jù)提取出來(lái)的），其高程的獲取所依賴(lài)的DEM必須是填洼后的無(wú)凹陷點(diǎn)DEM。

河流線屬于一維數(shù)據(jù)，DEM在平面上屬于二維數(shù)據(jù)，其相交之后會(huì)產(chǎn)生一個(gè)像元點(diǎn)的連續(xù)像元。如果了解過(guò)D8單流向算法的話(huà)，就可以很容易理解，如果直接從原始DEM中獲取像元值，其值有可能是一個(gè)匯點(diǎn)（凹陷點(diǎn)）。然而使用填洼后的DEM，可以確保該點(diǎn)與周邊一點(diǎn)的高程持平。

下圖是一個(gè)測(cè)試數(shù)據(jù)，黃色部分為河流線，點(diǎn)為河流的節(jié)點(diǎn)。包含兩份DEM，一份是原始的DEM，一份是填洼過(guò)的DEM(Fill_DEM) 。河流是使用了前面章節(jié)介紹到的河流提取方法提取出來(lái)的線狀河流。讀者可以參考《ArcGIS水文分析實(shí)戰(zhàn)教程（6）河流提取與河網(wǎng)分級(jí)》的做法去提取。
測(cè)試的主要內(nèi)容是分兩種情況打印河流各節(jié)點(diǎn)的高程值。一是獲取原始DEM的高程值，而是獲取填洼后的DEM的高程值。通過(guò)高程值進(jìn)行對(duì)比選用不同數(shù)據(jù)引起的差異。

第一張圖是原始DEM的沿河高程值，第二張圖是填洼后的DEM的沿河高程值?？梢钥吹绞褂迷糄EM的時(shí)候，有可能出現(xiàn)下游高程比上游高程還要高的情況，也就是說(shuō)這個(gè)上游地勢(shì)低的點(diǎn)剛好就是一個(gè)匯（凹陷點(diǎn)）。而填洼后的DEM則會(huì)將這個(gè)點(diǎn)給填平，確保不會(huì)斷流，從而使得下游高程最少不會(huì)高于上游高程。對(duì)比其他高程值，是沒(méi)有變化的，變化的部分就只是出現(xiàn)【匯】的點(diǎn)。

如果原始DEM的精度足夠高的話(huà)，這種情況一般很少碰到。但為了穩(wěn)妥起見(jiàn)，筆者還是建議對(duì)原始DEM做填洼，生成無(wú)凹陷點(diǎn)的DEM。讀者可以參考之前章節(jié) 《ArcGIS水文分析實(shí)戰(zhàn)教程（4）地形預(yù)處理》

如果所使用的河流是從DEM中提取出來(lái)的，那么在矢量河流提取的時(shí)候使用【柵格河流矢量化】時(shí)是否需要使用【簡(jiǎn)化折線】參數(shù)，這個(gè)參數(shù)是默認(rèn)選上的，如下圖所示：

不選擇簡(jiǎn)化折線，河流線會(huì)基本上就是按照流量線走向，河流折點(diǎn)密度會(huì)比較高；選擇簡(jiǎn)化的話(huà)，河流線的折點(diǎn)會(huì)被抽希一部分。這個(gè)怎么取舍？如果在地形變化比較大的情況下，不建議簡(jiǎn)化；如果地勢(shì)趨平，那就直接簡(jiǎn)化線。這種方式會(huì)直接影響到最終的比降數(shù)值。

下圖黃線部分是簡(jiǎn)化節(jié)點(diǎn)后的河流線，紅線是沒(méi)有簡(jiǎn)化節(jié)點(diǎn)的河流線，其形狀基本上是重合的，而且在關(guān)鍵的河流節(jié)點(diǎn)上，也是完全重合。

也就是說(shuō)，簡(jiǎn)化與不簡(jiǎn)化，其參與計(jì)算的河段數(shù)量是不一樣的。上圖數(shù)據(jù)中簡(jiǎn)化后的子河段數(shù)為32，不簡(jiǎn)化子河段的數(shù)量是92 。如果人為計(jì)算的話(huà)，那么可定是后者的計(jì)算量比較大了。

ArcGIS 除了前面章節(jié)提及到的水文分析工具及其相關(guān)工具之外，還有非常龐大的工具集對(duì)空間計(jì)算進(jìn)行支撐。由于沒(méi)有斷面數(shù)據(jù)的支撐，這里只考慮通過(guò)地形和河流線型進(jìn)行計(jì)算，也就是說(shuō)實(shí)現(xiàn)【約翰斯通-克羅斯】法。

根據(jù)【約翰斯通-克羅斯】法的計(jì)算公式，可以推導(dǎo)出其分析流程

通過(guò)河流線提取河流折點(diǎn)，利用這些折點(diǎn)與DEM數(shù)據(jù)疊加，獲取高程值

利用河流折點(diǎn)，獲取其坐標(biāo)值，根據(jù)兩點(diǎn)距離公式，計(jì)算每一段子河段的長(zhǎng)度

根據(jù)每個(gè)點(diǎn)的高程值，計(jì)算出相鄰兩點(diǎn)之間的高差

根據(jù)高差與子河段長(zhǎng)度，計(jì)算出各子河段的河底比降

匯總這些數(shù)據(jù)計(jì)算出整個(gè)河段的平均河底比降

在這5各步驟中，只有第1步與ArcGIS工具有關(guān)。河流的折點(diǎn)提取，可以使用【數(shù)據(jù)管理】--【要素】--【要素折點(diǎn)轉(zhuǎn)點(diǎn)】工具來(lái)實(shí)現(xiàn)，如下圖所示：

至于利用折點(diǎn)與DEM數(shù)據(jù)疊加，獲取其高程值，可以使用【數(shù)據(jù)管理】--【柵格】--【柵格屬性】--【獲取像元值】工具來(lái)獲取。這個(gè)工具每次只能實(shí)現(xiàn)一個(gè)坐標(biāo)點(diǎn)與DEM數(shù)據(jù)疊加查詢(xún)。如下圖

所以，為了方便計(jì)算，筆者采用ArcPy去寫(xiě)一個(gè)腳本工具，實(shí)現(xiàn)參數(shù)的靈活輸入。其中heightList就是用來(lái)保存高程值的數(shù)組。

    # 獲取河流的節(jié)點(diǎn)，以用作計(jì)算高差和河長(zhǎng)
    for row in arcpy.da.SearchCursor(outLine1, ["SHAPE@XY"]):
        # 從點(diǎn)圖層數(shù)據(jù)獲取所有的點(diǎn)的XY
        x, y = row[0]
        #從地形數(shù)據(jù)中獲取其高程的象元值
        pointStr = str(x) +" "+str(y)
        result = arcpy.GetCellValue_management(dem,pointStr)
        riverPoints.append(row[0])
        heightList.append(str(result))

第2-5步基本上就是一些數(shù)學(xué)的運(yùn)算，利用python就很好實(shí)現(xiàn)了。為了讓廣大水文研究者用起來(lái)，筆者將這段代碼公開(kāi)，并且制作成為工具。在代碼中筆者已經(jīng)預(yù)留了ArcGIS的工具參數(shù)，方便制作成為腳本工具。

河流平均比降【約翰斯通-克羅斯】法的實(shí)現(xiàn)代碼如下

#############################################################
#                                                           #
#   該工具利用【約翰斯通-克羅斯】方法計(jì)算河底平均比降           #
#   在使用的時(shí)候必須先選擇一條河流，工具會(huì)判斷輸入的河流數(shù)量     #
#   大于一條河流時(shí)不予以計(jì)算                                  #
#   用于計(jì)算的DEM為無(wú)凹陷點(diǎn)DEM                                #
#   BY 李遠(yuǎn)祥                                                #
#                                                            #
##############################################################


# -*- coding: gbk -*-
import arcpy
import math
inputLine =arcpy.GetParameterAsText(0) #河流圖層

dem = arcpy.GetParameterAsText(1) #無(wú)凹陷點(diǎn)DEM
riverLen = 7873.232038

arcpy.AddMessage(u"..................")
arcpy.AddMessage(u"......開(kāi)始執(zhí)行...........")

# 記錄高程值的數(shù)據(jù)
heightList = []

iputCout=0
mmCursor = arcpy.da.SearchCursor(inputLine, ["SHAPE@LENGTH"])
for row1 in mmCursor:
    riverLen = row1[0]      #順便獲取河流長(zhǎng)度
    iputCout=iputCout+1

arcpy.AddMessage(u"......獲取了河流長(zhǎng)度......")
#如果選中的是一個(gè)線要素，則執(zhí)行計(jì)算，否者不計(jì)算多條河流
if iputCout==1:
    ##先將線轉(zhuǎn)為折點(diǎn)
    outLine1 = arcpy.FeatureVerticesToPoints_management(inputLine, "in_memory/outline")
    arcpy.AddMessage(u"......河流線已經(jīng)轉(zhuǎn)換為點(diǎn)集.......")

    #河流節(jié)點(diǎn)
    riverPoints = []

    # 獲取河流的節(jié)點(diǎn)，以用作計(jì)算高差和河長(zhǎng)
    for row in arcpy.da.SearchCursor(outLine1, ["SHAPE@XY"]):
        # 從點(diǎn)圖層數(shù)據(jù)獲取所有的點(diǎn)的XY
        x, y = row[0]
        #從地形數(shù)據(jù)中獲取其高程的象元值
        pointStr = str(x) +" "+str(y)
        result = arcpy.GetCellValue_management(dem,pointStr)
        riverPoints.append(row[0])
        heightList.append(str(result))

    arcpy.AddMessage(u"......已經(jīng)獲取了折點(diǎn)的高程值......")
    #子河段長(zhǎng)度

    subRiverList = []
    hDiffList = []
    # 計(jì)算每一個(gè)子河段的長(zhǎng)度
    for i in range(0,len(riverPoints)-1):
        x,y = riverPoints[i]
        x1,y1 = riverPoints[i+1]
        #計(jì)算河段長(zhǎng)度
        dis =math.sqrt(pow((x1-x),2)+pow((y1-y),2))
        subRiverList.append(dis)

    arcpy.AddMessage(u"......河流子河段長(zhǎng)度計(jì)算完成完成.......")

    # 計(jì)算從高程從低到高部分的差距
    listLen = len(heightList)
    for j in range(0,listLen-1):
        h1 = float(heightList[j])
        h2 = float(heightList[j+1])
        xheight = 0
        xheight = h1-h2
        hDiffList.append(xheight)

    arcpy.AddMessage(u"......各折點(diǎn)的高程差計(jì)算完成......")

    # 計(jì)算每個(gè)河段的比降
    subbijiang = []
    for i in range(len(hDiffList)):
        Si= hDiffList[i]/subRiverList[i]
        subbijiang.append(Si)

    arcpy.AddMessage(u"......子河段比降計(jì)算完成......")
    # 計(jì)算最終的比降
    # 計(jì)算分子
    sumUp = 0
    for i in range(len(hDiffList)-1):
        sumUp = sumUp + subRiverList[i]*pow(subbijiang[i],0.5)

    S = pow((sumUp/riverLen),2)

    arcpy.AddMessage(u"......該河流的總長(zhǎng)為： " + str(riverLen)+u"   平均比降為: " +str(S))
    arcpy.AddMessage(u"......計(jì)算完畢.......")

else:
    arcpy.AddMessage(u"......您選中了 "+str(iputCout) +  u"  條河流......")
    arcpy.AddMessage(u"......工具已經(jīng)停止計(jì)算,沒(méi)有算出您想要的結(jié)果......")
    arcpy.AddMessage(u"......因?yàn)槟壳斑x中了多條河流。只能選中一條河流進(jìn)行計(jì)算!......")

從代碼中可以看到，這個(gè)工具針對(duì)的是一條河流進(jìn)行計(jì)算的。如果想指定河流的某一段進(jìn)行計(jì)算，那么就需要將河流轉(zhuǎn)為折點(diǎn)，選擇對(duì)應(yīng)的折點(diǎn)，生成新的河段才能參與計(jì)算。這個(gè)過(guò)程需要用到【點(diǎn)集轉(zhuǎn)線】工具進(jìn)行重構(gòu)要分析的河段。

可以將上面的python代碼制作成為腳本工具，方便在ArcGIS中進(jìn)行可視化的調(diào)用。工具的制作可以參考筆者另一篇博文《ArcMap自定義腳本工具制作》，里面有詳細(xì)的論述。

該工具用到的ArcGIS工具涉及到ArcInfo級(jí)別 和 Spatial 擴(kuò)展模塊，請(qǐng)先確保ArcGIS軟件具備相應(yīng)的授權(quán)。

制作好工具之后可以雙擊調(diào)用該工具了，【約翰斯通-克羅斯】法需要用到無(wú)凹陷點(diǎn)DEM和河流線型數(shù)據(jù)。

從上面代碼可以看到，工具會(huì)判斷選中的河流數(shù)量，如果數(shù)量大于1，就不作計(jì)算，確保計(jì)算的是某一河段。ArcGIS的交互性非常的優(yōu)秀，只要選中某一條河流，就會(huì)將選中的部分作為參數(shù)輸入。

對(duì)應(yīng)工具的流程，這里第一步需要選中要分析的一條河流。如下圖

接著雙擊調(diào)用工具，兩個(gè)參數(shù)分別是河流圖層和無(wú)凹陷點(diǎn)DEM。一定要注意，這個(gè)DEM是填洼過(guò)的DEM，不是原始DEM數(shù)據(jù)。如下圖

點(diǎn)擊確定就可以執(zhí)行計(jì)算結(jié)果，結(jié)果會(huì)在最終的提示框中找到，如下圖所示：

至此，復(fù)雜的河流平均比降就計(jì)算完畢了！

水文分析是一個(gè)復(fù)雜的過(guò)程，需要有深厚的水文學(xué)基礎(chǔ)作為背景。而GIS是一個(gè)通用工具，提供非常龐大的分析功能。二者要結(jié)合起來(lái)，不但需要學(xué)習(xí)兩個(gè)專(zhuān)業(yè)的知識(shí)，還需要具備知識(shí)整合的能力。

代碼和制作ArcGIS工具的方法已經(jīng)在文中貼出來(lái)了。讀者可以自行參考制作。當(dāng)然想省事的話(huà)，筆者也已經(jīng)上傳制作好的工具，可以點(diǎn)擊跳轉(zhuǎn)到 下載頁(yè)面 。為了加強(qiáng)交流，防止伸手，多動(dòng)手操作，筆者已經(jīng)將下載分值設(shè)置為最高分10分。 積分不夠的讀者請(qǐng)多動(dòng)手制作 ^_^ 
如果讀者對(duì)水文分析感興趣的話(huà)，可以持續(xù)關(guān)注CSDN的GIS制圖樂(lè)園，以及微信公眾號(hào)【GIS制圖樂(lè)園】。BY 李遠(yuǎn)祥