淹沒柔性植被河道阻力系數特性及計算方法研究

趙汗青，王偉杰，趙彥芳，4，馮達騫，4，李今今，徐宇軒，5

（1. 中國長江三峽集團有限公司，科學技術研究院，北京 100038； 2. 水資源高效利用與工程安全國家工程研究中心，江蘇 南京 210024； 3. 中國水利水電科學研究院，流域水循環模擬與調控國家重點實驗室，北京 100038；4. 河北工程大學水利學院，河北 邯鄲 056000； 5. 華北水利水電大學水利學院，河南 鄭州 450045）

0 引 言

植被是河道生態系統中不可或缺的組成部分，它具有許多生態功能。例如通過根系固定來保持河床穩定性，通過表皮的吸收能力來改善水質，并通過為生物提供附著基質和棲息地來豐富生物多樣性特征，因此植被被廣泛應用于河道生態系統的修復和重建。但植被的存在會阻礙水流運動，降低水流流速，抬升河道水位，影響河道的行洪能力。柔性植被的阻水機理更為復雜，在水流影響下表現出一定的運動和彎曲，從而導致植被頂部的水流更加紊亂，影響了河道局部水流的變化。因此，開展柔性植被的阻力特性研究對改善河道生態環境，計算河道防洪能力具有重要意義［1］。已有研究更多側重物理過程描述和基于基理的分析，如Wang［2］利用水槽實驗推斷出局部均勻流中孤立圓柱體的阻力系數（Cd）呈現出接近拋物線的形狀，在最密集的樹冠情況下達到峰值，并將這種Cd（x）變化的結果總結為體積阻力公式；Okamoto 和Nezu［3］使用粒子圖像測速儀（PIV）在有淹沒植被明渠水流中進行湍流測量；Yang 和Choi［4］提出了一種速度關系，用于預測植被中兩層流的速度分布；Wang［5］建立了描述沉水植物歸一化阻力與基于植物有效高度的雷諾數之間關系的方程，并用于計算不同生長階段沉水植物的植被阻力參數。但是大部分研究者沒有考慮植被本身形態的影響，尚未形成含植被形態變化的阻力系數計算公式，因此，將柔性植被的變形加入其中，提出了淹沒柔性植被河道阻力系數的計算方法。該方法可應用于不同水深、不同植被密度、不同植被高度等條件下的淹沒柔性植被河道阻力系數計算，其結果可對河流生態保護及防洪措施提供理論依據。

1 研究方法

1.1 理論分析

經典Darcy-Weisbach 公式由Darcy 和Weisbach 提出，現在被廣泛應用于明渠中摩擦損失的計算［6］。其中的達西-魏斯巴赫阻力系數f可表示為:

在有植被存在的情況下，水流與邊界產生的湍流剪應力遠小于由植被形狀所產生的剪應力，因此可忽略不計，故有:

式中:Cd為植被阻力系數；m為植被密度；D為植被桿徑；hv為植被高度；Uv為植被層水流運動速度。

聯立公式（1）、（2）可得:

其中:Lc=(CdmD)-1為調整長度尺度；

對于淹沒植被的情況，計算植被阻力系數的核心是確定植被層中的水流速度與整體斷面平均流速的比值。定義Us為淹沒植被上方的自由水層的斷面平均流速，定義ΔU=Us-Uv為自由水層流速與植被層流速的差。則有:

式中:hw為水深。

結合式（4）及ΔU的定義，有:

將式（5）代入式（3），得:

考慮到淹沒柔性植被會在水流作用下發生彎曲，故將柔性植被彎曲前的高度引入公式（6）。

式中:hc為柔性植被彎曲前的高度。

由式（7）可以看出，hc/Lc、hv/hc及hc/hw反映了植被阻力系數、植被分布密度、植被桿徑、植被高度、水深、彎曲前后形態變化的綜合作用；而且ΔU/Uv也隨著hc/Lc、hv/hc及hc/hw變化而變化。

由此，可得到三個無量綱因子hc/Lc、hv/hc及hc/hw。

定義:

Wang［7］等 根 據 一 維 閉 合 數 值 模 型 建 立ΔU和Uv的 相 關關系:

式中:c1、c2、c3均為參數。

由式（7）、（8）、（9）可建立f與α、β、γ之間的函數關系，其中，α、β、γ為自變量，f為因變量，得到f的一般公式為:

由式（10）可得f與α、β、γ三個參數相關關系為:

1.2 最大差異性算法

利用最大差異性算法（MDA）將原始數據進行選擇分類。MDA 的目標是從大小為N的數據庫中選擇大小為M的代表性子集。將N維向量組成的數據樣本X={x1，x2，…，xn}，通過該算法得到代表數據多樣性的向量子集{v1，…，vm}，通過從數據樣本{v1}傳輸一個向量來初始化子集，迭代選擇剩余的M-1 個元素，計算數據庫中每個剩余數據與子集元素之間的不相似性，并將最不相似的元素轉移到子集。當算法達到M次迭代時，該過程結束［8］。在該方法中，子集的初始數據作為數據樣本中相對于其他數據樣本的差異和最大的向量。在選擇數據的過程中，計算數據庫中每個剩余向量與子集中每個向量的相異度，并建立數據庫中每個向量與子集中每個向量的唯一相異度來定義最相異的一個。本文即從實驗數據中選擇出訓練組、測試組組和驗證組。

1.3 遺傳算法

遺傳算法是將達爾文遺傳規律和自然淘汰的生物進化過程作為理論依據，將原始數據作為初代種群，通過對交叉和變異過程的模擬，篩選并儲存適應度較好的個體，并對篩選過程重復進行，最終得到個體的最優解。遺傳算法的主要步驟為生成原始種群、選擇、交叉和變異［7］。利用集成遺傳算法的程序Eureqa來尋找公式，具體步驟為［9］:①輸入數據。將由最大差異性算法分類得到的訓練組和測試組數據依次輸入到程序中。②初始化。程序內部的符號函數生成器會隨機的將算子和算符結合起來。本文用到的算符有加、減、乘、除和冪。初始解不用于進化，當子表達式的變化不會導致解的精度超出程序給定的范圍時，該子表達式將會被廢除。③交叉及變異。運用選定的精度標準對表達式的預測值和驗證組的實測值進行比較，舍棄掉不符合標準的表達式，剩下的表達式利用程序內置的概率函數進行雜交，并生成新的表達式。④程序會生成一系列不同精度和復雜度的表達式，需要根據表達式表現的物理意義、精度和復雜度來選擇合適的表達式，隨即終止程序運行。

在尋找目標表達式的過程中需要遵循兩個原則:①為防止公式的過度擬合，需及時終止程序；②選擇表達式同時兼顧解的誤差和公式的復雜度［10］。

2 數據分析與公式構建

2.1 實驗數據整理與分析

利用文獻檢索，收集到不同條件下的試驗數據。試驗數據的相關信息如表1 所示。對收集到的數據進行篩選處理，最終得到89組實驗數據。使用最大差異性算法對數據進行選擇，將訓練組、測試組和驗證組的數據按照4∶4∶2的比例進行分配。3組實驗數據統計見表2。將訓練組和驗證組數據導入到Eureqa軟件中，進行公式擬合。利用Matlab 程序將實驗數據進行分類，分類結果如圖1所示。

圖1 不同數據組分布圖Fig.1 Distribution of different data groups

表1 實驗數據的相關信息列表Tab.1 List of relevant information of experimental data

表2 實驗數據統計Tab.2 Experimental data statistics

2.2 柔性植被阻力系數公式構建

根據最大差異性算法選擇的數據，再利用Eureqa 程序尋找計算公式。Eureqa 程序可以尋找原始數據里變量的關聯性，然后提出來一系列公式來描述出來。在公式搜索過程中，利用適應度函數評估每個候選解，并且還考慮了求解公式的復雜度。當公式中變量的個數、系數和公式中包含的運算符號和類型增加時，公式的復雜度也會隨之增加。軟件會在復雜度相同的情況下保存誤差較小的預測公式。本此公式的搜尋，有14組公式被保留，其中復雜度最小為1，最大為40。求解結果如表3所示。

表3 公式求解Tab.3 Formula solution

圖2 描述了Pareto 前沿，總體看來，誤差與復雜度之間大致呈負相關關系，當公式復雜度增加時，誤差會隨之變小。當公式較為簡單時，會導致計算值與實測值過度擬合，公式的物理意義不大；但當復雜度提升到一定程度時，模擬的精度提高的并不明顯。

圖2 Pareto前沿Fig.2 Pareto frontier

當公式的復雜度分別為16、18、20 時，復雜度由16 變到18時誤差的變化較大，但當復雜度由18 變為20 時誤差變化不明顯，且公式復雜度為18時，復雜度適中，擬合程度良好。綜上所述，選擇復雜度為18的公式作為最終公式，即:

式（12）公式擬合程度較好，表達形式相對簡單，可以準確、清晰的表達出淹沒柔性植被形態特征與阻力系數之間的關系。

在實際工程應用中，曼寧系數的使用相對較多，其與Dacy-Weisbach系數f存在定量關系:

式中:R為水力半徑；g為重力加速度，取9.8 m/s2。

聯立式（12）、（13）得到曼寧系數n與各個無量綱因子之間的定量關系。

2.3 公式計算值與實測值的對比分析

為了評價本文公式的擬合程度，將測試組的數據代入本文公式中來對阻力系數進行計算。采用相關系數（r）及均方根誤差（RMSE）評價f、n與各個無量綱因子之間的擬合程度。

相關系數（r）表示如下:

式中:Xi、Yi為兩個不同的變量（在本發明中分別表示為計算值和實測值）分別為變量Xi、Yi的均值；N為數據長度。

r值的絕對值介于0～1 之間。通常來說，r越接近1，表示X、Y兩個量之間的相關程度就越強，反之，r越接近于0，X、Y兩個變量之間的相關程度就越弱。

均方根誤差（RMSE）是用于表征計算值與實測值曲線的擬合程度，均方根誤差越小，擬合程度越高。公式表示為:

式中:Xi、Yi分別為計算值和實測值；N為數據長度。

將河道阻力系數f和曼寧系數n的公式計算值與試驗測量值進行統計分析，結果見圖3。

圖3 f、n實測值與計算值對比Fig3 Comparison of measured and calculated values of f and n

由圖3 可看出，達西-魏斯巴赫阻力系數f的公式計算值和實驗測量值的相關系數為0.96，均方根誤差為0.07；曼寧系數n的公式計算值和實驗測量值的相關系數為0.95，均方根誤差為0.005 9。由此可看出本文提出的淹沒柔性植被河道阻力系數f與曼寧系數n計算值與實測值擬合較好，即該公式可以很好的用來計算淹沒柔性植被的河道阻力系數。

3 結 論

本文將柔性植被的變形考慮其中，首先分析了達西-魏斯巴赫阻力系數與淹沒度及植被屬性之間的相關關系，然后將數據進行篩選處理，利用最大差異性算法將實驗數據進行選擇分類，運用遺傳算法尋找因變量和自變量之間的關系，再通過對表達式復雜度、誤差的分析，以及不同阻力系數之間的轉換關系，提出并驗證了淹沒柔性植被河道的達西-魏斯巴赫阻力系數公式及曼寧系數公式，為植被化生態河道設計提供理論基礎。