淹沒(méi)柔性植被河道阻力系數(shù)特性及計(jì)算方法研究

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

（1. 中國(guó)長(zhǎng)江三峽集團(tuán)有限公司，科學(xué)技術(shù)研究院，北京 100038； 2. 水資源高效利用與工程安全國(guó)家工程研究中心，江蘇 南京 210024； 3. 中國(guó)水利水電科學(xué)研究院，流域水循環(huán)模擬與調(diào)控國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100038；4. 河北工程大學(xué)水利學(xué)院，河北 邯鄲 056000； 5. 華北水利水電大學(xué)水利學(xué)院，河南 鄭州 450045）

0 引 言

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

1 研究方法

1.1 理論分析

經(jīng)典Darcy-Weisbach 公式由Darcy 和Weisbach 提出，現(xiàn)在被廣泛應(yīng)用于明渠中摩擦損失的計(jì)算［6］。其中的達(dá)西-魏斯巴赫阻力系數(shù)f可表示為:

在有植被存在的情況下，水流與邊界產(chǎn)生的湍流剪應(yīng)力遠(yuǎn)小于由植被形狀所產(chǎn)生的剪應(yīng)力，因此可忽略不計(jì)，故有:

式中:Cd為植被阻力系數(shù)；m為植被密度；D為植被桿徑；hv為植被高度；Uv為植被層水流運(yùn)動(dòng)速度。

聯(lián)立公式（1）、（2）可得:

其中:Lc=(CdmD)-1為調(diào)整長(zhǎng)度尺度；

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

式中:hw為水深。

結(jié)合式（4）及ΔU的定義，有:

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

考慮到淹沒(méi)柔性植被會(huì)在水流作用下發(fā)生彎曲，故將柔性植被彎曲前的高度引入公式（6）。

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

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

由此，可得到三個(gè)無(wú)量綱因子hc/Lc、hv/hc及hc/hw。

定義:

Wang［7］等 根 據(jù) 一 維 閉 合 數(shù) 值 模 型 建 立ΔU和Uv的 相 關(guān)關(guān)系:

式中:c1、c2、c3均為參數(shù)。

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

由式（10）可得f與α、β、γ三個(gè)參數(shù)相關(guān)關(guān)系為:

1.2 最大差異性算法

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

1.3 遺傳算法

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

在尋找目標(biāo)表達(dá)式的過(guò)程中需要遵循兩個(gè)原則:①為防止公式的過(guò)度擬合，需及時(shí)終止程序；②選擇表達(dá)式同時(shí)兼顧解的誤差和公式的復(fù)雜度［10］。

2 數(shù)據(jù)分析與公式構(gòu)建

2.1 實(shí)驗(yàn)數(shù)據(jù)整理與分析

利用文獻(xiàn)檢索，收集到不同條件下的試驗(yàn)數(shù)據(jù)。試驗(yàn)數(shù)據(jù)的相關(guān)信息如表1 所示。對(duì)收集到的數(shù)據(jù)進(jìn)行篩選處理，最終得到89組實(shí)驗(yàn)數(shù)據(jù)。使用最大差異性算法對(duì)數(shù)據(jù)進(jìn)行選擇，將訓(xùn)練組、測(cè)試組和驗(yàn)證組的數(shù)據(jù)按照4∶4∶2的比例進(jìn)行分配。3組實(shí)驗(yàn)數(shù)據(jù)統(tǒng)計(jì)見(jiàn)表2。將訓(xùn)練組和驗(yàn)證組數(shù)據(jù)導(dǎo)入到Eureqa軟件中，進(jìn)行公式擬合。利用Matlab 程序?qū)?shí)驗(yàn)數(shù)據(jù)進(jìn)行分類(lèi)，分類(lèi)結(jié)果如圖1所示。

圖1 不同數(shù)據(jù)組分布圖Fig.1 Distribution of different data groups

表1 實(shí)驗(yàn)數(shù)據(jù)的相關(guān)信息列表Tab.1 List of relevant information of experimental data

表2 實(shí)驗(yàn)數(shù)據(jù)統(tǒng)計(jì)Tab.2 Experimental data statistics

2.2 柔性植被阻力系數(shù)公式構(gòu)建

根據(jù)最大差異性算法選擇的數(shù)據(jù)，再利用Eureqa 程序?qū)ふ矣?jì)算公式。Eureqa 程序可以尋找原始數(shù)據(jù)里變量的關(guān)聯(lián)性，然后提出來(lái)一系列公式來(lái)描述出來(lái)。在公式搜索過(guò)程中，利用適應(yīng)度函數(shù)評(píng)估每個(gè)候選解，并且還考慮了求解公式的復(fù)雜度。當(dāng)公式中變量的個(gè)數(shù)、系數(shù)和公式中包含的運(yùn)算符號(hào)和類(lèi)型增加時(shí)，公式的復(fù)雜度也會(huì)隨之增加。軟件會(huì)在復(fù)雜度相同的情況下保存誤差較小的預(yù)測(cè)公式。本此公式的搜尋，有14組公式被保留，其中復(fù)雜度最小為1，最大為40。求解結(jié)果如表3所示。

表3 公式求解Tab.3 Formula solution

圖2 描述了Pareto 前沿，總體看來(lái)，誤差與復(fù)雜度之間大致呈負(fù)相關(guān)關(guān)系，當(dāng)公式復(fù)雜度增加時(shí)，誤差會(huì)隨之變小。當(dāng)公式較為簡(jiǎn)單時(shí)，會(huì)導(dǎo)致計(jì)算值與實(shí)測(cè)值過(guò)度擬合，公式的物理意義不大；但當(dāng)復(fù)雜度提升到一定程度時(shí)，模擬的精度提高的并不明顯。

圖2 Pareto前沿Fig.2 Pareto frontier

當(dāng)公式的復(fù)雜度分別為16、18、20 時(shí)，復(fù)雜度由16 變到18時(shí)誤差的變化較大，但當(dāng)復(fù)雜度由18 變?yōu)?0 時(shí)誤差變化不明顯，且公式復(fù)雜度為18時(shí)，復(fù)雜度適中，擬合程度良好。綜上所述，選擇復(fù)雜度為18的公式作為最終公式，即:

式（12）公式擬合程度較好，表達(dá)形式相對(duì)簡(jiǎn)單，可以準(zhǔn)確、清晰的表達(dá)出淹沒(méi)柔性植被形態(tài)特征與阻力系數(shù)之間的關(guān)系。

在實(shí)際工程應(yīng)用中，曼寧系數(shù)的使用相對(duì)較多，其與Dacy-Weisbach系數(shù)f存在定量關(guān)系:

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

聯(lián)立式（12）、（13）得到曼寧系數(shù)n與各個(gè)無(wú)量綱因子之間的定量關(guān)系。

2.3 公式計(jì)算值與實(shí)測(cè)值的對(duì)比分析

為了評(píng)價(jià)本文公式的擬合程度，將測(cè)試組的數(shù)據(jù)代入本文公式中來(lái)對(duì)阻力系數(shù)進(jìn)行計(jì)算。采用相關(guān)系數(shù)（r）及均方根誤差（RMSE）評(píng)價(jià)f、n與各個(gè)無(wú)量綱因子之間的擬合程度。

相關(guān)系數(shù)（r）表示如下:

式中:Xi、Yi為兩個(gè)不同的變量（在本發(fā)明中分別表示為計(jì)算值和實(shí)測(cè)值）分別為變量Xi、Yi的均值；N為數(shù)據(jù)長(zhǎng)度。

r值的絕對(duì)值介于0～1 之間。通常來(lái)說(shuō)，r越接近1，表示X、Y兩個(gè)量之間的相關(guān)程度就越強(qiáng)，反之，r越接近于0，X、Y兩個(gè)變量之間的相關(guān)程度就越弱。

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

式中:Xi、Yi分別為計(jì)算值和實(shí)測(cè)值；N為數(shù)據(jù)長(zhǎng)度。

將河道阻力系數(shù)f和曼寧系數(shù)n的公式計(jì)算值與試驗(yàn)測(cè)量值進(jìn)行統(tǒng)計(jì)分析，結(jié)果見(jiàn)圖3。

圖3 f、n實(shí)測(cè)值與計(jì)算值對(duì)比Fig3 Comparison of measured and calculated values of f and n

由圖3 可看出，達(dá)西-魏斯巴赫阻力系數(shù)f的公式計(jì)算值和實(shí)驗(yàn)測(cè)量值的相關(guān)系數(shù)為0.96，均方根誤差為0.07；曼寧系數(shù)n的公式計(jì)算值和實(shí)驗(yàn)測(cè)量值的相關(guān)系數(shù)為0.95，均方根誤差為0.005 9。由此可看出本文提出的淹沒(méi)柔性植被河道阻力系數(shù)f與曼寧系數(shù)n計(jì)算值與實(shí)測(cè)值擬合較好，即該公式可以很好的用來(lái)計(jì)算淹沒(méi)柔性植被的河道阻力系數(shù)。

3 結(jié) 論

本文將柔性植被的變形考慮其中，首先分析了達(dá)西-魏斯巴赫阻力系數(shù)與淹沒(méi)度及植被屬性之間的相關(guān)關(guān)系，然后將數(shù)據(jù)進(jìn)行篩選處理，利用最大差異性算法將實(shí)驗(yàn)數(shù)據(jù)進(jìn)行選擇分類(lèi)，運(yùn)用遺傳算法尋找因變量和自變量之間的關(guān)系，再通過(guò)對(duì)表達(dá)式復(fù)雜度、誤差的分析，以及不同阻力系數(shù)之間的轉(zhuǎn)換關(guān)系，提出并驗(yàn)證了淹沒(méi)柔性植被河道的達(dá)西-魏斯巴赫阻力系數(shù)公式及曼寧系數(shù)公式，為植被化生態(tài)河道設(shè)計(jì)提供理論基礎(chǔ)。