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

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

0 引 言

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

1 研究方法

1.1 理論分析

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

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

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

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

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

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

式中:hw為水深。

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

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

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

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

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

由此，可得到三個無量綱因子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與α、β、γ三個參數(shù)相關(guān)關(guān)系為:

1.2 最大差異性算法

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

1.3 遺傳算法

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

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

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

2.1 實驗數(shù)據(jù)整理與分析

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

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

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

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

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

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

表3 公式求解Tab.3 Formula solution

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

圖2 Pareto前沿Fig.2 Pareto frontier

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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