無(wú)資料或少資料區(qū)河流流量監(jiān)測(cè)與定量反演

李甲振，郭新蕾，鞏同梁，王 靜，王 濤，李 慧

（1.中國(guó)水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100038；2.西藏自治區(qū)水利廳，西藏 拉薩 850000；3.西藏自治區(qū)水文水資源勘測(cè)局，西藏 拉薩 850000；4.水利部水利水電規(guī)劃設(shè)計(jì)總院，北京 100120）

1 研究背景

世界范圍內(nèi)80%的人口和65%的河流生態(tài)系統(tǒng)受困于不穩(wěn)定的水資源供應(yīng)，但人們對(duì)河流流量這一關(guān)鍵影響因素的了解甚少［1］，只有不足60%的徑流是在入?？谶M(jìn)行了觀測(cè)，內(nèi)陸地區(qū)的監(jiān)測(cè)點(diǎn)更為稀疏［2］；受經(jīng)濟(jì)、政治因素影響，局部地區(qū)的水文測(cè)站數(shù)量自1980年代起就不斷減少［3］，一些水文站點(diǎn)也因缺乏維護(hù)管理而導(dǎo)致監(jiān)測(cè)數(shù)據(jù)中斷。對(duì)于無(wú)資料或少資料區(qū)，受限于惡劣的氣候、不便的交通或復(fù)雜的地形條件，水文站點(diǎn)的建設(shè)、維護(hù)以及徑流的監(jiān)測(cè)就更加困難了。這些地區(qū)一般是在眾多大江、大河源頭和上游區(qū)域的西部地區(qū)，比如青藏高原-亞洲水塔，對(duì)整個(gè)流域水循環(huán)具有重要意義，因此，如何高效準(zhǔn)確地對(duì)這些地區(qū)的河流流量進(jìn)行監(jiān)測(cè)成為水文研究需要解決的首要問(wèn)題之一。

遙感以其高效、實(shí)時(shí)、觀測(cè)范圍廣、數(shù)據(jù)量大等優(yōu)點(diǎn)，在區(qū)域徑流監(jiān)測(cè)中發(fā)揮了重要的作用?；谶b感技術(shù)獲取河道典型水力參數(shù)，反演河流流量，為上述問(wèn)題提供了一種解決思路。國(guó)內(nèi)外出現(xiàn)了地基高低頻雷達(dá)遙感監(jiān)測(cè)、航空攝影遙感監(jiān)測(cè)、多光譜衛(wèi)星遙感監(jiān)測(cè)、全數(shù)字超高頻雷達(dá)系統(tǒng)監(jiān)測(cè)等手段，可對(duì)河道的長(zhǎng)度和形狀、水面寬、表面流速、水力坡度、地表植被等信息進(jìn)行大數(shù)據(jù)量實(shí)時(shí)采集［4-5］，推動(dòng)了基于遙感的水文徑流監(jiān)測(cè)的發(fā)展；由美國(guó)、法國(guó)、加拿大以及英國(guó)等國(guó)家共同研制，并將于2021年發(fā)射的SWOT（Surface Water and Ocean Topography）衛(wèi)星，其最重要的一個(gè)目標(biāo)就是通過(guò)遙感獲取所需數(shù)據(jù)，反演河流徑流量［6］。根據(jù)流量反演所需水力特征參變量的個(gè)數(shù)，反演模型分為：?jiǎn)巫兞磕Ｐ汀㈦p變量模型和三變量模型［4，7-9］。雙變量和三變量模型中，引入了斷面平均流速、斷面平均水深和水力坡度這些參變量。斷面平均流速需根據(jù)遙測(cè)的表面流速進(jìn)行轉(zhuǎn)換，其率定方法尚未成熟，且表面流速的觀測(cè)受地面風(fēng)向、風(fēng)速等因素影響；另一方面，由于遙感技術(shù)對(duì)平面尺度的信息捕捉和分辨率優(yōu)于垂直方向，因此，基于水位信息的單變量模型和多變量模型應(yīng)用相對(duì)較少，基于水面寬反演河流流量的方法已成為該領(lǐng)域的研究熱點(diǎn)。

Leopold和Maddock［8］通過(guò)對(duì)美國(guó)中西部和西南地區(qū)河流的水面寬和流量關(guān)系進(jìn)行分析，首次建立了水面寬和流量之間關(guān)系的冪函數(shù)關(guān)系式。Park［10］研究了單站水力幾何模型指數(shù)的變化規(guī)律，指數(shù)b的變化范圍為0.00～0.59。在對(duì)單站水力幾何模型研究的過(guò)程中，Gleason和Smith［11］發(fā)現(xiàn)了單站水力幾何模型的系數(shù)與指數(shù)之間存在一定的函數(shù)關(guān)系，建立了多站水力幾何模型，對(duì)密西西比河、亞大巴斯卡河、長(zhǎng)江等河流典型河段進(jìn)行流量反演的相對(duì)誤差均值為20%～30%。之后，Gleason和Wang［12］給出了多站水力幾何模型的數(shù)學(xué)理論依據(jù)。由于上述模型是基于水力學(xué)及地理信息學(xué)的統(tǒng)計(jì)學(xué)方法，源數(shù)據(jù)的處理以及誤差的來(lái)源分析和控制也是研究的重點(diǎn)。Harman等［13］對(duì)澳大利亞?wèn)|南地區(qū)河流進(jìn)行了研究，水面寬、水深、流速和水力坡度的不確定度分別為4%、3%、7%和23%。Afshari等［14］提出了一種辨識(shí)現(xiàn)場(chǎng)調(diào)研數(shù)據(jù)隨機(jī)誤差和系統(tǒng)誤差的遞歸過(guò)濾程序。針對(duì)河流整體連續(xù)信息獲取較為困難這一情況，Pavelsky［15］利用空間不連續(xù)的衛(wèi)星影像構(gòu)建了塔納納河多個(gè)水面寬——流量關(guān)系曲線。國(guó)內(nèi)相關(guān)研究大多是針對(duì)某一河流水力幾何模型的系數(shù)和指數(shù)取值及其變化規(guī)律，例如，吳保生和李凌云［16］、馬元旭和許炯心［17］以及冉立山等人［18］分別研究了黃河下游不同河段、無(wú)定河及其各支流以及黃河上游不同斷面的水力幾何模型的系數(shù)和指數(shù)取值及變化特征。趙春紅和高建恩［19］給出了不同侵蝕溝（細(xì)溝、淺溝、切溝和沖溝）的水力幾何模型的指數(shù)。

針對(duì)無(wú)資料或少資料區(qū)的徑流監(jiān)測(cè)難題，通過(guò)水力分析建立水面寬與流量函數(shù)關(guān)系的一般表達(dá)式，結(jié)合遙測(cè)遙感技術(shù)，提出一種水力學(xué)與遙感相耦合的徑流反演方法，為無(wú)資料或少資料區(qū)的徑流監(jiān)測(cè)提供解決方案。

2 流量反演理論

河流流量計(jì)算公式為：

式中：Q為流量，m3/s；w為水面寬，m；d為斷面平均水深，m；v為斷面平均流速，m/s。

Leopold和Maddock［8］首次建立的水面寬、斷面平均水深和平均流速與流量之間的函數(shù)關(guān)系式為：

式中：a、c和k為系數(shù)；b、f和m為指數(shù)。

將式（2）—式（4）代入式（1），a、c、k和b、 f、m滿足如下關(guān)系式：

典型的河道過(guò)流斷面如圖1所示。圖1（a）是常見(jiàn)于山區(qū)河流的一種過(guò)流斷面形態(tài)。河道的縱比降大、水流流速快、沖刷力較強(qiáng)，水深、水面寬和過(guò)流面積隨流量的波動(dòng)變化大，具有V形河槽的特征。圖1（b）是常見(jiàn)于山區(qū)和平原河流的一種斷面形態(tài)，由于河道走勢(shì)或地質(zhì)的原因，其中一側(cè)邊坡陡，另一側(cè)邊坡緩，兩側(cè)邊坡不同。對(duì)于平原區(qū)、丘陵區(qū)的河流，河道縱比降小，水流流速緩，過(guò)流斷面的寬深比大，常見(jiàn)的斷面形態(tài)如圖1（c）和（d）所示。當(dāng)河道中存在不能忽略的突起時(shí)，過(guò)流斷面如圖1（e）所示。

圖1的（a）（b）可近似概化為三角形，其不同在于，圖1（a）是左右對(duì)稱的，是圖1（b）的一種特殊情況。因此，本文將圖1（a）（b）概化為圖3所示的三角形斷面進(jìn)行分析。圖1（c）（d）（e）所示的寬淺型河道，水面寬與水深的關(guān)系可近似概化為冪函數(shù)型，函數(shù)形式為 y=αxβ。假定過(guò)流斷面的水面寬為100 m，水深為10 m，α和 β取表1給出的數(shù)值時(shí)，斷面如圖2所示?？梢钥闯?，根據(jù)冪函數(shù)y=αxβ與給定的α和 β值所繪制的斷面形狀，與實(shí)際工程中的一些寬淺型河道斷面是相近的。因此，工程中的寬淺型斷面是可以利用冪函數(shù)進(jìn)行擬合的，只不過(guò)系數(shù)α和指數(shù) β會(huì)隨斷面尺寸發(fā)生改變。

圖1 河道過(guò)流斷面的不同形態(tài)

表1 α和β取值

圖2 河道概化斷面示意圖

2.1 三角形斷面 典型的三角形斷面如圖3所示，其過(guò)水?dāng)嗝婷娣eA的計(jì)算公式為：

圖3 三角形斷面

式中：h為水深，m。

過(guò)水?dāng)嗝鏉裰?χ的計(jì)算公式為：

式中 θ1和 θ2為角度，°。

因此，水力半徑R的計(jì)算公式為：

根據(jù)曼寧公式：

式中：J為水力坡度；n為曼寧糙率。

斷面流速為：

對(duì)于某一河道過(guò)流斷面，水力坡度J、糙率n和頂角θ1、θ2均為常數(shù)，因此，流速v與正比例相關(guān)，即：

平均水深d的計(jì)算公式為：

因此，平均水深d與w正比例相關(guān)，即：

聯(lián)立式（12）、式（14）和式（2）—式（4）可得：

式（15）和式（16）代入式（5）可得，三角形斷面的指數(shù)分別為b=0.375，f=0.375，m=0.25。

2.2 冪函數(shù)型斷面 典型的冪函數(shù)型斷面如圖4所示，斷面方程為：

圖4 冪函數(shù)型斷面

過(guò)水?dāng)嗝婷娣eA的計(jì)算公式為：

由于寬淺河道的濕周 χ近似等于河寬w［20］，因此，水力半徑的計(jì)算公式為：

斷面平均水深d的計(jì)算公式為：

因此，斷面平均水深d與wβ正比例相關(guān)，即：

聯(lián)立式（20）、式（22）和式（2）—式（4）可得：

式（23）和式（24）代入式（5）可得冪函數(shù)型斷面的3個(gè)指數(shù)分別為：

當(dāng) β=2時(shí)，河道斷面為拋物線型，根據(jù)式（25）計(jì)算的3個(gè)指數(shù)b、 f和m分別為0.23、0.46和0.31 ，這與文獻(xiàn)［21］給出的結(jié)果是一致的。當(dāng) β=1時(shí)，斷面方程為 y=αx，也就是說(shuō)，過(guò)流斷面為三角形斷面。由式（25）計(jì)算所得的3個(gè)指數(shù)b、 f和m分別為0.375、0.375和0.25，這與2.1的推導(dǎo)結(jié)果是一致的。

2.3 流量反演方法 由于現(xiàn)階段遙測(cè)技術(shù)可較為準(zhǔn)確捕捉的過(guò)流斷面水力特征參數(shù)是水面寬，因此，本文構(gòu)建的流量反演方法是基于遙測(cè)的水面寬。某一時(shí)刻的水面寬w0和流量Q0與t時(shí)刻的水面寬wt和流量Qt代入式（2），可得如下關(guān)系式

聯(lián)立式（26）和式（27）得：

式（28）即為無(wú)資料或少資料區(qū)流量反演的理論基礎(chǔ)。它反映了任一時(shí)刻的流量Qt與幾個(gè)參數(shù)（水面寬wt、指數(shù)b、計(jì)算基數(shù)w0和Q0）之間的函數(shù)關(guān)系式，在確定指數(shù)b和計(jì)算基數(shù)w0、Q0后，即可直接依據(jù)遙測(cè)的水面寬進(jìn)行流量反演。這就解決了傳統(tǒng)設(shè)立水文測(cè)站進(jìn)行測(cè)流所面臨的站點(diǎn)建設(shè)和運(yùn)營(yíng)維護(hù)困難的問(wèn)題，僅需進(jìn)行一次測(cè)量就可以實(shí)現(xiàn)未來(lái)河流流量的遠(yuǎn)程實(shí)時(shí)監(jiān)測(cè)。因此，本文提出了一種基于遙測(cè)、遙感與水力學(xué)模型相耦合的河流流量監(jiān)測(cè)與定量反演方法。對(duì)于無(wú)資料區(qū)，也就是沒(méi)有任何可參閱水文數(shù)據(jù)的地區(qū)，具體步驟為：

（1）選取河流斷面，建設(shè)地面增強(qiáng)標(biāo)識(shí)，作為遙感影像識(shí)別特征點(diǎn)。

（2）利用遙感或其他技術(shù)獲得河流典型位置的斷面特征，包括斷面形態(tài)、水位H、水深D和水面寬w0，進(jìn)而計(jì)算獲得過(guò)水?dāng)嗝婷娣e A0、濕周 χ0和水力坡降J0。

（3）根據(jù)曼寧公式，計(jì)算當(dāng)前水面寬w0對(duì)應(yīng)的河道流量Q0：

式中n為糙率。糙率系數(shù)n根據(jù)工程界積累的不同邊界材料及使用情況下的取值表格確定［22-23］。

（4）采用冪函數(shù)方程 y=αxβ對(duì)實(shí)測(cè)斷面進(jìn)行擬合，確定系數(shù)α和冪指數(shù) β后，即建立了水面寬與水深、過(guò)流面積等參數(shù)間一一對(duì)應(yīng)的函數(shù)關(guān)系；根據(jù)冪指數(shù) β，由式（25）計(jì)算指數(shù)b。

（5）利用無(wú)人機(jī)或衛(wèi)星得到的高精度遙感影像并結(jié)合地面增強(qiáng)標(biāo)識(shí)點(diǎn)的尺寸辨識(shí)該斷面任一時(shí)刻的水面寬wt。

（6）根據(jù)式（28）反演未來(lái)時(shí)刻的河道流量Qt。

針對(duì)少資料區(qū)，也就是以前有過(guò)水文測(cè)站或數(shù)據(jù)記錄的地區(qū)，過(guò)流斷面形狀和計(jì)算基數(shù)可查閱相關(guān)資料獲取，通過(guò)步驟（1）、（4）、（5）和（6）進(jìn)行流量反演。

3 反演模型應(yīng)用實(shí)例

利用模型實(shí)驗(yàn)和野外實(shí)測(cè)資料，對(duì)提出的無(wú)資料或少資料區(qū)河流流量監(jiān)測(cè)與定量反演方法的應(yīng)用進(jìn)行展示和校驗(yàn)，分析反演方法的精度。

3.1 模型實(shí)驗(yàn) 模型實(shí)驗(yàn)在中國(guó)水利水電科學(xué)研究院水力學(xué)試驗(yàn)室展開(kāi)，實(shí)驗(yàn)平臺(tái)為自循環(huán)系統(tǒng)，主要包括變速泵供水系統(tǒng)、有機(jī)玻璃水槽和水位控制尾門等，如圖5所示。變速泵供水系統(tǒng)由水泵和變頻器組成，可調(diào)節(jié)水槽來(lái)流量；有機(jī)玻璃水槽長(zhǎng)17.3 m、寬1.2 m、高0.25 m，其過(guò)流斷面為復(fù)式斷面；水位控制尾門采用推拉式結(jié)構(gòu)，通過(guò)調(diào)節(jié)推拉門的開(kāi)度控制下游水位。水槽來(lái)流量通過(guò)電磁流量計(jì)監(jiān)測(cè)，水面寬采用卷尺進(jìn)行測(cè)量。

利用冪函數(shù)形式對(duì)復(fù)式斷面進(jìn)行擬合，復(fù)式斷面及其擬合斷面如圖6所示，函數(shù)關(guān)系式為y=1.65298x3.78366。由式（25）得，擬合斷面的三個(gè)指數(shù)分別為b=0.13687、 f=0.51788和m=0.34525。

圖5 模型實(shí)驗(yàn)平臺(tái)

圖6 實(shí)際斷面與擬合斷面

本文提出的方法是基于遙測(cè)的水面寬進(jìn)行流量反演的，應(yīng)用實(shí)例中，模型實(shí)驗(yàn)及野外實(shí)測(cè)均使用了測(cè)量的水面寬作為遙測(cè)數(shù)據(jù)。模型實(shí)驗(yàn)測(cè)試了9種不同工況的水面寬和流量，以工況5的水面寬w5和流量Q5作為計(jì)算基數(shù)w0和Q0，利用式（28）反演其他8種工況的流量，反演流量和實(shí)測(cè)流量如圖7所示。其中，6種工況的反演流量與實(shí)測(cè)流量的相對(duì)誤差小于20%。分別以工況1～9的水面寬和流量作為計(jì)算基數(shù)，對(duì)其他8種工況的流量進(jìn)行反演，反演流量與實(shí)測(cè)流量的相對(duì)誤差如圖8所示。以工況1的實(shí)測(cè)水面寬w1和流量Q1作為計(jì)算基數(shù)w0和Q0，反演其他工況流量時(shí)，反演流量均小于實(shí)測(cè)流量，且相對(duì)誤差大于20%；同樣地，以工況9的實(shí)測(cè)水面寬w9和流量Q9作為計(jì)算基數(shù)

w0和Q0，反演其他工況流量時(shí)，反演流量均大于實(shí)測(cè)流量，且誤差較大。分別以工況2～8的實(shí)測(cè)水面寬和流量作為計(jì)算基數(shù)反演其他工況流量時(shí)，相對(duì)誤差均值為18.77%；56組反演流量中，40組與實(shí)測(cè)流量的相對(duì)誤差小于20%，所占比重為71.43%。

圖7 實(shí)測(cè)流量與反演流量

圖8 反演流量與實(shí)測(cè)流量的相對(duì)誤差

3.2 野外實(shí)測(cè) 奴各沙站位于雅魯藏布江上游河段，左岸為山體，右岸為河灘地，上下游500 m范圍內(nèi)水流平穩(wěn)，斷面均勻，河道水面與岸邊對(duì)比度清晰。第二次青藏高原綜合科學(xué)考察水利資源與工程科考分隊(duì)于2017年11月在該站搭建了遙感水文增強(qiáng)野外標(biāo)識(shí)設(shè)施，并對(duì)典型斷面進(jìn)行了現(xiàn)場(chǎng)測(cè)量。奴各沙站的過(guò)流斷面與圖1的斷面（b）相近，可將其近似為三角形斷面，過(guò)流斷面及其擬合斷面如圖9所示。根據(jù)2.1的推導(dǎo)，三角形斷面的水面寬與流量?jī)绾瘮?shù)關(guān)系式的指數(shù)b為0.375。

實(shí)測(cè)的74組水面寬與流量數(shù)據(jù)如圖10所示。以水面寬為基準(zhǔn)，對(duì)上述實(shí)測(cè)的74組數(shù)據(jù)進(jìn)行排序，水面寬最小的工況記為1，水面寬最大的工況記為74。分別以工況1、10、20、……、60和70的實(shí)測(cè)水面寬和流量作為計(jì)算基數(shù)，利用式（28）反演其他73種工況的流量，反演流量與實(shí)測(cè)流量如圖11所示。以工況1和工況10的實(shí)測(cè)水面寬和流量作為計(jì)算基數(shù)進(jìn)行流量反演時(shí)，反演流量小于實(shí)測(cè)流量；以工況70的實(shí)測(cè)水面寬和流量作為計(jì)算基數(shù)時(shí)，反演流量大于實(shí)測(cè)流量。以相對(duì)誤差小于20%和30%為基準(zhǔn)，計(jì)算使用不同計(jì)算基數(shù)的反演流量合格率，結(jié)果如圖12所示。以工況20、30、40、50和60的數(shù)據(jù)為計(jì)算基數(shù)時(shí)，相對(duì)誤差均值為20.71%；64.66%的反演流量與實(shí)測(cè)流量的相對(duì)誤差小于20%，86.03%的反演流量與實(shí)測(cè)流量的相對(duì)誤差小于30%。

圖9 實(shí)測(cè)斷面與擬合斷面

圖10 實(shí)測(cè)流量與水面寬

圖11 反演流量與實(shí)測(cè)流量

圖12 反演流量合格率

4 結(jié)論

針對(duì)現(xiàn)階段由于復(fù)雜地形、惡劣氣候或經(jīng)濟(jì)、政治等因素造成的無(wú)資料或少資料區(qū)所面臨的河流流量監(jiān)測(cè)難題，將常見(jiàn)的過(guò)流斷面進(jìn)行概化，通過(guò)水力分析構(gòu)建了一種水力學(xué)與遙感相耦合的徑流反演方法。

（1）對(duì)于三角形斷面，水面寬、斷面平均水深和斷面平均流速與流量之間函數(shù)關(guān)系式的三個(gè)指數(shù)分別為0.375、0.375和0.25；對(duì)于方程為 y=αxβ的冪函數(shù)型斷面，三個(gè)指數(shù)的取值分別為

（2）無(wú)資料區(qū)河流流量監(jiān)測(cè)與定量反演方法包括：①建設(shè)地面增強(qiáng)標(biāo)識(shí)點(diǎn)；②實(shí)測(cè)斷面特征，包括斷面形態(tài)、水位H、水深D和水面寬度w0，計(jì)算過(guò)流面積 A0、濕周 χ0、水力坡度J0；③根據(jù)河床及植被情況確定糙率系數(shù)n，計(jì)算過(guò)流量Q0；④采用冪函數(shù)方程 y=αxβ對(duì)實(shí)測(cè)斷面進(jìn)行擬合，確定系數(shù)α和冪指數(shù) β，即建立了水面寬與水深、過(guò)流面積等參數(shù)一一對(duì)應(yīng)的函數(shù)關(guān)系，根據(jù)冪指數(shù) β計(jì)算指數(shù)b；⑤監(jiān)測(cè)水面寬wt；⑥反演t時(shí)刻的流量Qt。對(duì)于少資料區(qū)，根據(jù)已有資料可確定斷面形態(tài)和計(jì)算基數(shù)w0、Q0時(shí)，按照步驟①、④、⑤和⑥進(jìn)行流量反演。

（3）根據(jù)模型實(shí)驗(yàn)水面寬反演的56組流量，相對(duì)誤差平均值為18.77%；其中，40組反演流量與實(shí)測(cè)流量的相對(duì)誤差小于20%，占比為71.43%。根據(jù)野外實(shí)測(cè)水面寬反演的流量，相對(duì)誤差均值為20.71%；其中，64.66%的反演流量與實(shí)測(cè)流量的相對(duì)誤差小于20%，86.03%的相對(duì)誤差小于30%。

（4）提出的方法主要適用于河道過(guò)流斷面可近似為三角形或冪函數(shù)型的情形，且具有一定的精度；在現(xiàn)階段水文數(shù)據(jù)匱乏的無(wú)資料或少資料區(qū)，該方法可期成為一種有效的流量監(jiān)測(cè)手段。當(dāng)斷面為其他形態(tài)或易沖淤時(shí)，可能會(huì)產(chǎn)生較大的誤差，這也是進(jìn)一步研究需要解決的問(wèn)題。