劉墨陽，蔣四維，林云發(fā)，張生穩(wěn)，晉 濤，黃迎春， 蘭蓓蓓，周曉英，蔣飛卿，童冰星， 周文靜

(1.長江水利委員會水文局漢江水文水資源勘測局，湖北 襄陽 441025； 2.河海大學水文水資源學院，江蘇 南京 210098)

為提高流量自動化監(jiān)測水平，國內(nèi)部分水文站陸續(xù)投入定點式聲學多普勒流速剖面儀(ADCP)進行比測。常見的定點式ADCP為水平式，也稱水平聲學多普勒流速剖面儀(H-ADCP)。該儀器安裝在水中某一高度，能連續(xù)不斷地發(fā)射聲波信號，獲取其安裝高程對應水層沿橫斷面流速分布，具有精度高、速度快、數(shù)據(jù)傳輸頻次高、數(shù)據(jù)量大的特點[1-2]。

H-ADCP測得的斷面某一水層流速分布與斷面平均流速有相關性但不等價，因此數(shù)據(jù)獲取后需經(jīng)過分析處理才能用于推求斷面流量。斷面流量推求方法主要有數(shù)值法和代表流速法，數(shù)值法也稱流速剖面法，該方法基于H-ADCP獲取的流速數(shù)據(jù)，根據(jù)理論垂向流速分布，推算得到全斷面流量[3]；代表流速法也稱指標流速法，通過分析流速網(wǎng)格數(shù)據(jù)關系找到合適區(qū)間得到代表流速，將代表流速與該測站常規(guī)測流方法得到的斷面平均流速建立函數(shù)關系，進而推求流量。代表流速法實際應用時通過代表流速推求斷面平均流速，與相應時刻水位下查算出的斷面面積相乘得到斷面流量。國內(nèi)部分水文站如高壩洲[4]、宜昌[5]、黃陵廟[6]、長沙壩[7]、小河壩[8]等基于代表流速法建立了比較簡單實用的推流公式。這些應用成果表明，實測斷面平均流速與代表流速有較強的線性相關性，一般用一元一次或者二次方程就能達到滿足規(guī)范和生產(chǎn)要求的推流精度。

近年來由于測流斷面受上下游水利工程調(diào)蓄影響，部分水文站點測流斷面水情復雜，代表流速與斷面平均流速的關系無法采用代表流速-斷面平均流速線性方程或一元二次方程描述；此外，其他影響因素如洪水漲落、回水頂托、連續(xù)洪水壅水、分洪潰口等引起水面比降改變，影響斷面流速時空分布的均勻性，導致斷面流速的非線性變化，增大了H-ADCP應用的難度。為解決該問題，國內(nèi)外不少學者做了相關嘗試，主要有兩種途徑：①結合先進測量儀器更好地控制斷面流速監(jiān)測質(zhì)量。例如：韋立新等[9]將H-ADCP和垂向定點式ADCP相結合，分別得到水平平均流速和垂線平均流速，建立多元線性回歸模型，得到南京水文實驗站不同流速級代表流速與斷面平均流速的關系；肖林[10]為解決滇池口流量在線監(jiān)測問題，在新運糧河布設流量在線監(jiān)測系統(tǒng)，通過非接觸式雷達波傳感器與H-ADCP分別監(jiān)測表面流速和水下代表流速，采用監(jiān)測數(shù)據(jù)模擬斷面流速分布。②探索新方法率定基于代表流速的斷面平均流速計算公式。這類方法主要是改進擬合方法，充分利用多種數(shù)據(jù)分析方案得到合理的推流公式。例如：杜興強等[11]分不同水位級率定了高壩洲水文站的代表流速與斷面平均流速的關系，建立了與低枯水水位-流量相結合的推流方案；李輝等[12]考慮斷面水位變化，引入變化水位流速標定系數(shù)，建立了虛擬梯形斷面平均流速與測驗斷面平均流速相關關系，并與直接建立代表流速與斷面平均流速定線關系的傳統(tǒng)方法進行對比，在上海蘇州河溫州路水文站取得了較好的應用效果；徐剛等[13-14]采用小波分析法進行濾波處理，通過BP神經(jīng)網(wǎng)絡建立斷面平均流速與H-ADCP代表流速關系,實現(xiàn)了復雜水流河段(三峽—葛洲壩段)的流量自動測報及資料整編；鄧山等[15]從粗大誤差剔除、線性回歸模型損失函數(shù)、代表單元格選取三方面優(yōu)化代表流速法，并以小河壩站為例，證明了所提出的優(yōu)化方法可顯著提高模型精度；袁德忠等[16]采用支持向量機、BP神經(jīng)網(wǎng)絡、極限學習機等機器學習方法，根據(jù)清泉溝水文站H-ADCP數(shù)據(jù)模擬其斷面流速分布，探究了機器學習方法與傳統(tǒng)水文測驗結合的可能；Hoitink等[17]結合代表流速法和流速剖面法，采用邊界層模型基于點流速計算局部流量，通過回歸模型得到斷面流量并應用于感潮河段；Sassi等[18-19]在Hoitink工作的基礎上，進一步考慮邊壁與水位的影響推求斷面流量，并在受變動回水影響的內(nèi)陸河流進行了應用推廣。

從現(xiàn)有研究可以看出引入機器學習以及深度學習等新方法不失為提升H-ADCP在生產(chǎn)實踐中應用精度的可行途徑，然而已有研究方法需要大量數(shù)據(jù)支撐，若參與訓練樣本量較小，模型容易出現(xiàn)過擬合現(xiàn)象，且構建推流方案時甚至會出現(xiàn)震蕩現(xiàn)象。為解決該問題，本文采用較為簡單實用的公式結構，在不同水情下選擇合適的回歸模型確定模型參數(shù)，探索了復雜水情下受水利工程影響測站的H-ADCP推流方案，可為實現(xiàn)流量在線監(jiān)測提供參考。

1 研究方法

本文通過構建以下兩種模型來研究復雜水情下受水利工程影響測站的H-ADCP推流方案：①除代表流速外，引入落差、H-ADCP儀器入水深等相關因素，建立斷面平均流速與代表流速、落差、H-ADCP儀器入水深等相關因素之間的多元線性回歸模型，根據(jù)最小二乘法得到模型系數(shù)；②考慮所有相關性較強的流速單元格，即每個參與率定的流速單元格在公式中具有不同的權重，應用LASSO回歸模型求解參數(shù)，從而擬合出多因素影響下的斷面平均流速回歸方程。此外，率定模型前，分析H-ADCP流速與斷面平均流速相關性，為挑選流速區(qū)間提供依據(jù)，根據(jù)實際水情選擇相應模型建立推流公式，得到研究時段的推流結果并進行精度評估。

1.1 最小二乘回歸模型

最小二乘回歸模型遵循以殘差平方和最小確定直線位置的原則，除了計算比較方便外，得到的估計量還具有無偏性、最小方差性等優(yōu)良特性，但這種方法對異常值非常敏感。最小二乘法是求解線性回歸模型應用最廣泛的方法，其優(yōu)化的目標函數(shù)為

(1)

式中：y1為輸出，代表斷面平均流速；x1為輸入，代表流速、落差、水深等；β1為線性回歸方程的系數(shù)。該方法試圖找到一條直線，使所有樣本到直線上的歐氏距離之和最小，可求得平方損失函數(shù)的極值點，進而得到待估參數(shù)[20]。在變量數(shù)目小于樣本數(shù)目時，應用效果較好。

H-ADCP代表流速大多只能反映安裝高度對應水層的流速情況。當河流流態(tài)穩(wěn)定且水情單一時，可通過H-ADCP單元格流速和斷面平均流速建立相關關系推流。然而水情復雜時，僅采用單一固定區(qū)間代表流速作為自變量無法滿足推流精度，也很難反映測驗斷面受頂托影響的程度和上下游水利工程狀況等，因此有必要綜合考慮多種水力因素的影響。本文參考落差指數(shù)法原理，除選取代表流速外，加入總落差、上游落差、下游落差等因素反映蓄放水及回水頂托影響，引入H-ADCP入水深反映本斷面流速分布代表性，建立多元線性回歸模型，采用最小二乘法求解模型參數(shù)。

1.2 LASSO回歸模型

LASSO回歸模型最初由Tibshirani[21]于1996年提出，目的是為了消除最小二乘回歸中的冗余變量，適用于參數(shù)數(shù)目縮減和參數(shù)選擇。在水文氣象領域，Hammami等[22-23]最早將其應用于模型變量選擇和模型簡化，而在國內(nèi)的應用較為少見。

LASSO回歸模型損失函數(shù)為L1范數(shù)，其優(yōu)化的目標函數(shù)為

(2)

式中：y2為實測輸出序列矩陣；x2為實測輸入序列矩陣；λ為控制LASSO回歸模型復雜度的正則參數(shù)，一般通過迭代選擇最佳的λ，本文采用網(wǎng)格搜索優(yōu)化該參數(shù)；β2為模型回歸系數(shù)。

LASSO回歸模型的L1正則化項為向量中各個元素絕對值之和，導致?lián)p失函數(shù)在零點位置不可導，其常見求解方法為坐標軸下降法和最小角回歸法，本文采用坐標軸下降法求解回歸系數(shù)β2[24]。坐標軸下降法沿某一坐標方向進行搜索，固定其他坐標方向，找到函數(shù)的局部極小值，不需要對目標函數(shù)求導，計算效率較高，最大迭代次數(shù)設置為2 000。

在流量較小時，引入多個水力因素后仍然無法較好擬合實測流量樣本。小流量時流速較小且流場紊亂，如果將固定區(qū)間流速取均值后可能無法反映實際情況，故將所有相關性較強的單個流速單元格納入考慮，采用LASSO回歸模型求解。

1.3 評價指標

采用相對偏差、系統(tǒng)誤差、隨機不確定度等評價指標評價流量模擬精度，并進行符號檢驗、適線檢驗和偏離數(shù)值檢驗[25]。除評估實測點擬合效果外，還收集上下游水文站點及水利工程資料進行逐日平均流量對比，并對照上下游水文站點的年徑流量，綜合評估推流成果合理性。

2 應用實例

2.1 研究區(qū)概況

白河水文站是漢江流域的重要水文站點，以往水位流量關系呈一條或多條有規(guī)律可循的臨時曲線。近年來，上游受白河電站蓄放水調(diào)蓄影響，且上游200 m處新建白鄖大橋影響測驗斷面流速分布，下游受新建孤山電站試驗性蓄水頂托影響，水情復雜，白河水文站水位流量關系點散亂(圖1)，特別是小流量時期，同水位級對應多個流量，無法應用臨時曲線法推流。

圖1 白河水文站7—12月實測水位流量

白河水文站于2020年7月3日安裝H-ADCP，選用600 kHz換能器收集數(shù)據(jù)，設置128個單元格，單元格寬度為1 m，安裝高程為172.50 m。8月左右白河進入主汛期，出現(xiàn)流量峰值4 000 m3/s左右的洪水過程，伴隨著含沙量增大，600 kHz換能器聲波能量衰減較快，其剖面有效測量范圍縮短至原有的1/3，所收集到的流速數(shù)據(jù)代表性減弱，無法滿足數(shù)據(jù)采集需要，改為300 kHz換能器，設置90個單元格，單元格寬度為2m。隨著下游孤山電站蓄水水位逐漸升高，H-ADCP安裝相對位置不斷下降，經(jīng)分析發(fā)現(xiàn)該水層流速數(shù)據(jù)逐漸失去代表性，10月14日將儀器安裝高程提高至173.77 m，單元格設置為120個。安裝位置及儀器測量范圍如圖2所示。

圖2 白河水文站H-ADCP安裝位置及測量范圍

2.2 數(shù)據(jù)選擇及預處理

選取白河水文站及其上游大王溝水位站、下游將軍河水位站、下游孤山橋水位站2020年8—12月水文資料作為研究數(shù)據(jù)。數(shù)據(jù)預處理步驟如下：

a.比測數(shù)據(jù)選擇。選取走航式ADCP流量測次，數(shù)據(jù)選取遵循以下原則：①回放測驗時的走航式ADCP原始數(shù)據(jù)，采用在H-ADCP正常工作時段內(nèi)的測次；②當走航式ADCP在GGA模式下流量小于BTM(底跟蹤)模式下流量時，選用BTM模式計算流量；③對于在水情平穩(wěn)時測得的數(shù)據(jù)，根據(jù)SL337—2006《聲學多普勒流量測驗規(guī)范》要求，施測 2 個測回，任一半測回BTM或GGA模式下測量值與平均值的相對誤差不應大于5%，否則補測同向的一個測次流量。當?shù)谝粶y回任一次BTM或GGA模式下測量值與平均值的相對誤差小于2%時，可不施測第二測回，流量取施測第一測回的平均值。滿足以上條件的視為閉合流量，否則做不閉合流量處理。將走航式 ADCP 測流流量除以相應水位對應過水斷面面積，得到斷面平均流速作為實測數(shù)據(jù)，為下一步公式率定做準備。

b.流速單元格選取。采用走航式ADCP流量對應起始時間段內(nèi)H-ADCP流速網(wǎng)格數(shù)據(jù)(剔除偽值)計算每個流速單元格在該時間段的流速均值；分析每個單元格流速與實測斷面平均流速的相關關系，計算兩者相關系數(shù)，選取與實測流量斷面平均流速相關性較好的單元格區(qū)間，準備好模型率定時所需H-ADCP流速數(shù)據(jù)。

c.計算相關水力因素。引入上下游相關輔助水位站點同時刻的摘錄水位(若時間不對應，對摘錄水位做線性插值處理)計算落差，用于反映上下游水利工程蓄放水、回水頂托對斷面流速的影響；計算H-ADCP儀器入水深和斷面平均水深，用于反映變化水位下相對水深對斷面流速分布的影響。

2.3 推流方案

推流方案流程如圖3所示。當H-ADCP流速數(shù)據(jù)沿橫斷面具有較好的分布規(guī)律時，采用最小二乘回歸模型推流較為簡便；當H-ADCP流速數(shù)據(jù)沿橫斷面分布正負不定時，采用LASSO回歸模型推流。該方法引入單個流速單元格，待率定參數(shù)個數(shù)較多。

圖3 推流方案流程

a.最小二乘回歸法。先采用互相關系數(shù)確定代表流速預選區(qū)間[20]，選擇相關系數(shù)較高的網(wǎng)格區(qū)間計算區(qū)間代表流速均值；然后分析模型自變量與因變量相關性，為確定公式項提供依據(jù)，并以提高擬合精度為目標，逐步增加自變量個數(shù)，采用最小二乘回歸法率定回歸模型。

b.LASSO回歸法。白河水文站流量較小時測流斷面流速分布較紊亂，單個網(wǎng)格數(shù)值在±0.2 m/s之間波動，H-ADCP流速區(qū)間均值代表性不佳，使用最小二乘回歸法很難得到滿足要求的推流精度，因此改變思路，將斷面相關性較強的單元格流速數(shù)據(jù)放入LASSO回歸模型中進行訓練。該方法在模型計算時并未采用流速區(qū)間均值，而是考慮所有相關性較強的單個流速單元格，即每個參與率定的流速單元格在公式中具有不同的權重。此外，將相關性較好的單個流速網(wǎng)格數(shù)據(jù)與儀器入水深、平均水深、白河孤山橋落差等相關因素組合，進一步提高擬合精度。

需要注意的是，兩種推流方法均屬于多元線性回歸模型求解方法，主要區(qū)別是對代表流速處理方法不同、公式形式不同(LASSO回歸法引入單個流速網(wǎng)格，其公式項數(shù)較多)和模型參數(shù)率定方法不同。

2.4 結果分析

由于儀器型號、儀器安裝位置改變以及水情變化的影響，H-ADCP獲取的流速資料不具有一致性，故根據(jù)儀器型號、水情和安裝位置，將資料劃分為7月3日到8月9日(600 kHz換能器，安裝高程172.50 m)、8月10日至9月24日(300 kHz換能器，安裝高程172.50 m)、9月25日至10月13日(300 kHz換能器，安裝高程172.50 m)和10月14日至12月31日(300 kHz，安裝高程173.77 m)4個時段。其中，8月10日至10月13日儀器型號與安裝高程未改變，但考慮到9月22—24日受孤山電站調(diào)節(jié)影響，白河水文站水位下降，退至H-ADCP安裝高度以下，9月24日后電站恢復蓄水，斷面水位上升至正常蓄水位，H-ADCP開始正常運行，水情有較大改變，故以9月24日為分界點再次進行時段劃分。

600 kHz換能器正常運行時間較短，因此主要對8月10日至12月31日收集的后3個時段300 kHz流速數(shù)據(jù)進行分析推算。

2.4.1時段1(8月10日至9月24日)

8月10日至9月24日共測流49次，測次均閉合，實測斷面平均流速范圍為0.547～2.19 m/s。去除測流期間H-ADCP流速數(shù)據(jù)異常或中斷的測次，共選用34次實測流量進行相關關系擬合。

擬合前對H-ADCP單元格流速數(shù)據(jù)進行相關性分析，計算每個單元格與實測斷面平均流速相關性(圖4)，以此作為單元格挑選依據(jù)。

圖4 時段1單元格流速與實測斷面平均流速相關系數(shù)

由圖4可知，從單元格20開始，相關性逐漸增強，經(jīng)過主泓段后，相關性逐漸減弱，到單元格67以后相關性銳減，流速相關性系數(shù)分布與斷面形狀較為一致。選取相關系數(shù)在0.8以上的單元格42～57區(qū)間流速均值作為代表流速，該區(qū)間基本位于主泓段，代表性較好，同時考慮上下游落差、斷面形狀、換能器入水深等相關因素，通過分析這些因素與斷面平均流速的相關性，以提高擬合精度為原則逐步引入相關變量，與斷面平均流速進行多元回歸分析。擬合公式為

(3)

表1 時段1單元格系數(shù)

使用擬合公式推算時段內(nèi)逐時流量，并進行平滑摘錄，得到該時段流量過程線。以實測流量作為真值，推算流量過程線相應時段流量均值作為推求流量(下同)，進行誤差檢驗。時段1使用擬合關系推算流量的系統(tǒng)誤差為0.41%，隨機不確定度為8.68%，符號檢驗、適線檢驗、偏離數(shù)值檢驗均通過，滿足規(guī)范要求。

2.4.2時段2(9月25日至10月13日)

9月25日至10月13日共測流15次，其中閉合流量14次，實測斷面平均流速范圍為0.011～0.808 m/s。去除測流期間H-ADCP流速數(shù)據(jù)缺失或異常的測次后，共有閉合流量9次，不閉合流量1次。由于該時段樣本數(shù)量較少，無法擬合出精度較好的相關關系，因此采用單測回流量進行擬合。同樣，首先進行單元格流速相關性分析，如圖5所示。

圖5 時段2單元格流速與實測斷面平均流速相關系數(shù)

從圖5可以看出，靠近左岸的主泓段相關性較好，中間部分受上游在建白鄖大橋橋墩影響相關系數(shù)較為散亂，因此選取左岸相關性較好的14～33單元格進行分析計算。因該時段總體流量較小，上下游落差變幅較小，使用時段1的分析方法擬合時精度不佳。為提高擬合精度，將14～33單元格流速值乘以H-ADCP相對水深，增強其代表性，并采用LASSO回歸法進行優(yōu)化求解。擬合公式為

(4)

式中：bj為H-ADCP流速14～33單元格系數(shù)，取值見表2；vind,j為H-ADCP單元格j的流速。由表2可知，參與率定的20個單元格系數(shù)不為0的共有17個。

表2 時段2單元格系數(shù)

對時段2擬合的相關關系進行誤差分析和檢驗，對于不閉合流量測次，檢驗時取測次均值進行誤差統(tǒng)計(下同)。推算流量的系統(tǒng)誤差為-0.55%，隨機不確定度為9.02%，符號檢驗、適線檢驗、偏離數(shù)值檢驗均通過，滿足規(guī)范要求。

2.4.3時段3(10月14日至12月31日)

10月14日調(diào)整H-ADCP安裝高度后共測流39次，其中閉合流量34次，不閉合流量5次，實測斷面平均流速范圍為0.037～0.799 m/s，其中11月13日測流期間H-ADCP異常(數(shù)據(jù)中斷)，不參與后續(xù)分析。各單元格流速序列與實測斷面平均流速的相關系數(shù)如圖6所示。

圖6 時段3單元格流速與實測斷面平均流速相關系數(shù)

由圖6可以看出，提高H-ADCP安裝高度后，斷面中間部分受上游在建白鄖大橋影響，流速分布依然散亂。選取相關系數(shù)0.6以上單元格區(qū)間6～40、46～55、76～93進行分析，分析方法與時段2相同，擬合公式如下：

(5)

式中：bj、bm、bn分別為H-ADCP單元格6～40、46～55、76～93對應的系數(shù)，取值見表3；vind,j、vind,m、vind,n分別為H-ADCP單元格j、m、n的流速。

由表3可知，參與率定的63個單元格系數(shù)不為0的一共為37個，可見LASSO回歸模型在保證計算精度的同時有效地減少了變量數(shù)目，并且主泓段對應的流速單元格參數(shù)權重較大，與實際情況相符。

表3 時段3單元格系數(shù)

對相關關系進行誤差分析和檢驗，推算流量的系統(tǒng)誤差為0.73%，隨機不確定度為9.37%，符號檢驗、適線檢驗、偏離數(shù)值檢驗均通過，滿足規(guī)范要求。

2.4.4推流成果

為消除流速脈動影響[26]，應用H-ADCP流速數(shù)據(jù)推流時，需對H-ADCP原始流速數(shù)據(jù)進行平滑處理后代入率定公式推流。平滑處理公式為

(6)

8月10日12:10之前，采用連時序法繪制水位-流量關系線進行推流，之后(H-ADCP正常運行期間)，由構建的相應時段公式采用連實測流量過程線法推流，H-ADCP中斷期間采用連時序法繪制水位-流量關系線推流。

將白河水文站8月10日后逐日平均流量與上游來水逐日平均流量(白河電站+長沙壩)進行對照(圖7)，兩者趨勢相應性較好。

圖7 上游來水逐日平均流量與白河水文站逐日平均流量對比

白河水文站洪水水文要素摘錄流量和實測流量如圖8所示,實測流量基本位于逐日平均流量過程線之上。從上下游水量來看，2020年白河水文站年徑流量233.7億m3，白河水文站下游黃家港水文站以上各站合計來水238.1億m3，白河水文站下游黃家港(二)站年徑流量271.6億m3。黃家港水文站以上區(qū)間各站集水面積合計82 708 km2，占黃家港(二)站集水面積的87%，上游來水量238.1億m3，占黃家港(二)站年徑流量的88%，上下游水量基本平衡，說明推流成果合理。

圖8 白河水文站洪水水文要素摘錄流量與實測流量對比

推流結果主要受到以下兩個不確定性因素的影響：①下游電站試驗性蓄水的影響。白河水文站水情變化較大，H-ADCP安裝高程并非一直保持在較為理想的位置。儀器使用時需要根據(jù)流速剖面分布[27]，合理選擇安裝高程，以更好地監(jiān)測代表流速。②應用LASSO回歸模型時，需要確定算法本身的相關參數(shù)，參數(shù)優(yōu)選范圍以及優(yōu)化方法對模型計算結果可能有一定的影響。

3 結 語

本文分析了不同水情下斷面平均流速與H-ADCP流速數(shù)據(jù)的相關關系，考慮落差、儀器入水深等因素建立最小二乘回歸模型，解決受水利工程影響測站的推流問題；引入代表性較好的單個流速網(wǎng)格，利用LASSO回歸模型進行參數(shù)估計，解決小流量下推流問題，構建了白河水文站2020年8月10日至12月31日H-ADCP推流方案?；诒疚姆桨竿扑愕牧髁砍晒到y(tǒng)誤差、隨機不確定度均滿足SL247—2012《水文資料整編規(guī)范》精度要求，并通過符號檢驗、適線檢驗和偏離數(shù)值檢驗；從日水量和年徑流量來看，推算流量滿足上下游水量平衡，有較好應用效果。本文所構建的推流方案可為白河水文站及類似受水利工程影響測站的H-ADCP推流方案構建提供參考。

H-ADCP流速代表性以及對應的實測流量成果精度直接影響推流關系率定，決定推流成果質(zhì)量。在水利工程影響下，測流斷面選擇、儀器參數(shù)設置以及安裝位置等均會影響數(shù)據(jù)收集成果，進而影響分析成果。未來需規(guī)范比測流程，依據(jù)定點式ADCP相關測量規(guī)范，確定H-ADCP安裝高程和有效工作范圍，保證數(shù)據(jù)采集質(zhì)量，收集不同水情下數(shù)據(jù)進行分析對照，更好地控制H-ADCP流速代表性。同時嘗試采用新方法應對水利工程等影響下的復雜水情，得到推流精度高且穩(wěn)健的方案。