譚 鵬，張?jiān)姾?，徐宏亮，劉宏寬，?虹，周文婷

(1.中電建生態(tài)環(huán)境集團(tuán)有限公司，深圳 518102；2.武漢大學(xué)水利水電學(xué)院，武漢 430072)

21世紀(jì)以來(lái)，隨著中國(guó)經(jīng)濟(jì)的飛速發(fā)展，我國(guó)各大河流的水環(huán)境問(wèn)題日益嚴(yán)峻[1-4]，其中河道水動(dòng)力水質(zhì)模擬作為預(yù)測(cè)以及水污染防治的重要手段之一其作用越來(lái)越凸顯重要性[5]。河道水下地形數(shù)據(jù)是進(jìn)行水質(zhì)模擬的基本條件，地形數(shù)據(jù)的質(zhì)量直接影響水質(zhì)模型模擬精度，地形處理不當(dāng)甚至直接導(dǎo)致模擬計(jì)算的失敗。當(dāng)前，我國(guó)河流實(shí)時(shí)的數(shù)字化工作還相當(dāng)不足，基本處于碎段化狀態(tài)，難以適應(yīng)高速經(jīng)濟(jì)發(fā)展帶來(lái)的環(huán)境問(wèn)題的需求，更不能滿足當(dāng)今我國(guó)對(duì)流域整體實(shí)施水環(huán)境精細(xì)化的管理目標(biāo)。所以，在缺乏實(shí)測(cè)地形資料的地區(qū)如何獲取滿足模型運(yùn)行所需要的河道水下地形數(shù)據(jù)，成為推廣模型的重要限制條件。河道水下地形的測(cè)量，早期采用的方法主要有：光學(xué)定位方法、無(wú)線電定位方法以及深度測(cè)量方法[6]。目前采用比較多的方法有：遙感影像技術(shù)、水聲定位、側(cè)掃聲吶、機(jī)載激光測(cè)量和無(wú)人船測(cè)量[7-12]。受測(cè)量條件和工作時(shí)間限制，水動(dòng)力水質(zhì)模型在缺乏實(shí)測(cè)水下地形資料時(shí)，一般會(huì)采用豐水期和枯水期水面線做概化處理。隨著遙感影像技術(shù)的飛速發(fā)展，遙感影像分辨率越來(lái)越高，由10 m提高到3 m。這為沒(méi)有條件進(jìn)行水下地形實(shí)測(cè)的地區(qū)使用水動(dòng)力水質(zhì)模型提供了一種快速便捷的方法[13]。

本研究提出了一種基于Google Earth等遙感資料，利用普通克里金插值法[14]插值出河道地形[15]，通過(guò)該插值地形構(gòu)建水動(dòng)力水質(zhì)模型[16]。模型可以運(yùn)用于地形相對(duì)較為缺乏的地區(qū)，且模擬結(jié)果能夠達(dá)到一定的精度要求，可在資料缺乏區(qū)域，尤其在工程前期方案的比選中發(fā)揮重要的作用。

1 研究方法

1.1 Google Earth獲取地形數(shù)據(jù)

對(duì)于沒(méi)有河道邊界及缺乏重點(diǎn)斷面實(shí)測(cè)資料的河段，本次研究采用的是通過(guò)利用Google Earth等衛(wèi)星地圖工具，找到目的河道并勾取河道邊界和重點(diǎn)斷面，得到一系列邊界散點(diǎn)和斷面散點(diǎn)的經(jīng)緯度坐標(biāo)(WGS84坐標(biāo)系)及對(duì)應(yīng)的高程信息，并將各散點(diǎn)的經(jīng)緯度坐標(biāo)轉(zhuǎn)換為模型計(jì)算所需的平面坐標(biāo)。

利用Google Earth操作得到邊界散點(diǎn)和斷面散點(diǎn)經(jīng)緯度的如圖1，首先利用Google Earth找到本次研究的研究區(qū)域。隨后，使用Google Earth中的添加路徑與添加地標(biāo)功能勾勒研究區(qū)域的邊界線以及重點(diǎn)斷面(如圖2)，采集的數(shù)據(jù)格式如圖3所示，數(shù)據(jù)內(nèi)容包括經(jīng)緯度以及高程，最后導(dǎo)出邊界點(diǎn)與斷面散點(diǎn)的經(jīng)緯度及高程。

轉(zhuǎn)化公式如公式(1)～(3)所示：

X=(N+H)cosBcosL

(1)

Y=(N+H)cosBsinL

(2)

Z=[N+(1-e2)+H]sinB

(3)

式中：B為緯度，rad；L為經(jīng)度，rad；H為大地高，m；e為地球的第一偏心率；N為地球面卯酉圈曲率半徑。e和N可分別按如下公式求解：

(4)

(5)

式中：a為地球長(zhǎng)半軸，a=6 378 137.0 m；b為地球短半軸，b=6 356 752.314 2 m。

1.2 普通克里金插值

普通克里金插值法(Ordinary Kriging)是克里金插值法的一種，是一種依據(jù)協(xié)方差函數(shù)對(duì)隨機(jī)過(guò)程進(jìn)行插值的回歸算法，該方法認(rèn)為可依據(jù)區(qū)域上的一系列點(diǎn)X1,X2,…,Xn的觀測(cè)值，采用一個(gè)線性組合來(lái)估計(jì)區(qū)域化變量在X0處的值，即：

(6)

滿足插值的無(wú)偏性和最優(yōu)性，即要求：

E[Z(X0)-Z*(X0)]=0

(7)

D[Z(X0)-Z*(X0)]=min

(8)

(9)

轉(zhuǎn)化為λi求解的矩陣形式為：

(10)

其中Cij可通過(guò)變差函數(shù)γ(h)來(lái)確定，二維情況下，空間點(diǎn)P在X軸和Y軸上均有變化，可分別對(duì)x和y方向進(jìn)行一維情況下的變差函數(shù)計(jì)算，即：

(11)

(12)

通常變差函數(shù)有球狀模型、高斯模型、冪函數(shù)模型等，本研究采用球狀模型：

(13)

式中：a為變程，即變量之間具有相關(guān)性的區(qū)域范圍；h為滯后距；c為基臺(tái)值，即變差函數(shù)在h大于變程時(shí)的值，為塊金值c0和拱高cc之和。

塊金值c0為變量純隨機(jī)性的部分，即在很短的距離內(nèi)有較大的空間變異性，無(wú)論h多小，兩個(gè)隨機(jī)變量都不相關(guān)；拱高cc為觀測(cè)值的變異幅度范圍。變差函數(shù)、變程、滯后距、基臺(tái)值、拱高、塊金值之間的關(guān)系如圖4所示。

1.3 二維水動(dòng)力水質(zhì)模型

本研究采用二維圣維南方程和二維污染物遷移轉(zhuǎn)化方程模擬寬深比較大、斷面變化較大的大型河道水動(dòng)力學(xué)特征和水質(zhì)遷移轉(zhuǎn)化規(guī)律，即：

(1)二維圣維南方程包含連續(xù)方程和動(dòng)量方程。

連續(xù)方程：

(14)

式中：H為微小水體的水深，m；u為x方向的流速，m/s；v為y方向的流速，m/s。

動(dòng)量方程：

(15)

(16)

式中：g為重力加速度，m/s2；ρ為水體密度，kg/m3；C為謝才系數(shù)，m0.5/s；f為柯氏力常數(shù)，f=2Ωsinφ，φ為緯度，Ω為地轉(zhuǎn)角速度；ξx、ξy分別為x、y方向上的渦動(dòng)黏滯系數(shù)，m2/s。

邊界條件：

上游邊界條件或入流邊界：

Q=Qin(t)

(17)

下游邊界條件或出流邊界：

z=zout(t)

(18)

初值條件：

φ=φin，t=0

(19)

式中：Q為流量，m3/s；z為水位，m；φ為Q，z中任何一個(gè)變量。

考慮對(duì)流、擴(kuò)散、降解等項(xiàng)，二維水質(zhì)模型基本方程如下：

(2)二維污染物遷移轉(zhuǎn)化方程。

河流中的水質(zhì)成分滿足如下的守恒方程：

(20)

式中：ci為第i種物質(zhì)的斷面平均濃度，mg/L；Ex和Ey為河道的縱向和橫向擴(kuò)散系數(shù)，m2/s；SLi、SBi和SKi為第i種物質(zhì)側(cè)流入?yún)R負(fù)Z荷、河道中源漏項(xiàng)，以及因水體水生態(tài)系統(tǒng)動(dòng)力轉(zhuǎn)換項(xiàng)，m·mg/(L·s)。

邊界條件：

上游邊界條件或入流邊界：

C=Cin(t)

(21)

下游邊界條件或出流邊界：

(22)

初值條件：

Ci=Cin，t=0

(23)

1.4 計(jì)算方法

首先，以河道研究可能涉及的最大水位，可以是歷史高 洪水位或最高堤岸線，以此作為計(jì)算的外邊界，通過(guò)網(wǎng)格生成程序生成模型計(jì)算的最大范圍的計(jì)算網(wǎng)格。然后，依據(jù)圣維南方程，通過(guò)初始和邊界條件，利用來(lái)流和下游水位及側(cè)流入?yún)R過(guò)程條件，根據(jù)橫斷面的水位變化，通過(guò)模型迭代計(jì)算，自動(dòng)識(shí)別，判斷出網(wǎng)格的干濕狀態(tài)后，濕網(wǎng)格參與計(jì)算，干網(wǎng)格不參與計(jì)算，依次即可完成水動(dòng)力水質(zhì)模型中相應(yīng)每個(gè)濕網(wǎng)格上的水位、流速和水質(zhì)濃度的計(jì)算。

其中谷歌地球提取出來(lái)的高程，在校正插值后，構(gòu)建出研究區(qū)域的水下地形，作為模型計(jì)算的邊界條件。

2 研究區(qū)域

2.1 區(qū)域概況

研究河段為長(zhǎng)江干流九江段，如圖5所示，其上邊界斷面位于賽城閘上游4 km處，下邊界斷面為九江烏石磯糧庫(kù)，長(zhǎng)約20 km，寬約1.5～2.5 km。研究河段內(nèi)有鶴問(wèn)湖污水處理廠和老鸛塘污水處理廠兩大污染源。

2.2 河道地形生成

依據(jù)上文所講到的地形獲取方法，根據(jù)九江河段的實(shí)際位置，使用Google Earth選取河道邊界并沿河道均勻選取30個(gè)斷面，得到河道邊界和斷面散點(diǎn)坐標(biāo)，如圖6所示。

利用獲取的河道邊界，繪制出平面二維正交網(wǎng)格，網(wǎng)格對(duì)應(yīng)實(shí)際地形尺寸約為30 m×30 m，如圖7所示。

依據(jù)河道邊界散點(diǎn)以及選取的斷面散點(diǎn)，結(jié)合網(wǎng)格結(jié)點(diǎn)坐標(biāo)，利用普通克里金法插值計(jì)算每個(gè)網(wǎng)格結(jié)點(diǎn)處的高程信息，生成模型計(jì)算所需的地形標(biāo)準(zhǔn)文件。

如圖8為根據(jù)遙感影像生成的河道地形與實(shí)測(cè)地形的對(duì)比圖，圖8(a)為根據(jù)遙感影像生成的地形，圖8(b)為九江水文局2017年實(shí)測(cè)水下地形。從圖8中可以看出，根據(jù)遙感影像生成的水下地形，盡管在趨勢(shì)上與實(shí)測(cè)地形基本一致，但是在河道深泓線的表征方面還是不如實(shí)測(cè)地形連貫合理。

2.3 模型計(jì)算參數(shù)

本次研究以COD和NH3-N兩種污染物為水質(zhì)指標(biāo)進(jìn)行模擬計(jì)算，依據(jù)本河段水力條件，參考該河段相關(guān)研究成果，模型各參數(shù)取值范圍如表1所示。

表1 水動(dòng)力水質(zhì)模型參數(shù)率定范圍Tab.1 Calibration range of hydrodynamic water quality model parameters

鶴問(wèn)湖污水處理廠和老鸛塘污水處理廠，排水規(guī)模分別為10和8 萬(wàn)t/d，兩污水處理廠排污口處排放濃度及江段背景濃度如表2所示。

表2 污染物排放濃度及背景濃度 mg/L

2.4 結(jié)果對(duì)比

2.4.1 水動(dòng)力對(duì)比

本次研究區(qū)域內(nèi)，共有CS1、九江站、CS2 3個(gè)水文測(cè)驗(yàn)斷面，如圖9所示。其驗(yàn)證數(shù)據(jù)為九江水文局在2008年9月6日的水文測(cè)驗(yàn)實(shí)測(cè)水文資料。

本研究選取這3個(gè)斷面的實(shí)測(cè)水位、流速進(jìn)行模型驗(yàn)證。其中，CS1、九江站、CS2 3個(gè)水文測(cè)驗(yàn)斷面平均水位驗(yàn)證結(jié)果及誤差如表3所示，流速驗(yàn)證結(jié)果及誤差如表4所示。

表3 兩種地形條件下水位模擬結(jié)果誤差統(tǒng)計(jì)表Tab.3 Statistics of water level simulation results under two terrain conditions

由表3和表4可以看出，基于插值地形的模型在絕大多實(shí)測(cè)點(diǎn)的相對(duì)誤差值都比基于實(shí)測(cè)地形的模型要大，其中基于實(shí)測(cè)地形模型的相對(duì)誤差都在10%以下，而本次研究所構(gòu)建的基于插值地形的模型相對(duì)誤差均在20%以下，兩種地形所構(gòu)建的模型都能夠較好的模擬河道水動(dòng)力規(guī)律。

表4 兩種地形條件下模型流速模擬結(jié)果誤差統(tǒng)計(jì)表Tab.4 Error statistics of simulation results of model velocity under two terrain conditions

2.4.2 水質(zhì)對(duì)比

模型水質(zhì)驗(yàn)證根據(jù)九江水文局2015年9月2日實(shí)測(cè)水質(zhì)數(shù)據(jù)進(jìn)行驗(yàn)證，選取研究河段內(nèi)的抗洪廣場(chǎng)、河西水廠、第三水廠(河?xùn)|水廠)、烏石磯糧庫(kù)為水質(zhì)測(cè)驗(yàn)斷面，如圖10所示。

水質(zhì)驗(yàn)證輸入條件如表5所示。

表5 水質(zhì)驗(yàn)證輸入條件Tab.5 Input conditions for water quality verification

基于兩種地形的模型在抗洪廣場(chǎng)、河西水廠、第三水廠、烏石磯斷面4個(gè)水質(zhì)監(jiān)測(cè)斷面處的COD及NH3-N斷面平均濃度模擬結(jié)果及誤差如表6、7所示。

表6 COD濃度模擬結(jié)果誤差統(tǒng)計(jì)表Tab.6 Error statistics of COD concentration simulation results

表7 NH3-N濃度模擬結(jié)果誤差統(tǒng)計(jì)表Tab.7 Error statistics of NH3-N concentration simulation results

從表6和表7可知，本次研究基于插值地形模型在水質(zhì)模擬方面，相對(duì)誤差普遍大于基于實(shí)測(cè)地形的模型，插值地形的模型結(jié)果最大相對(duì)誤差出現(xiàn)在河西水廠的NH3-N模擬上，達(dá)到了23.71%，小于30%；實(shí)測(cè)地形模型的最大相對(duì)誤差同樣出現(xiàn)在河西水廠的NH3-N模擬上，達(dá)到了16.88%，小于20%。

綜合對(duì)比水動(dòng)力模擬精度以及水質(zhì)模擬精度，本次研究基于插值地形的模型模擬精度不如實(shí)測(cè)地形，但是其在驗(yàn)證斷面處，水位模擬誤差小于0.2 m，流速模擬相對(duì)誤差小于20%，水質(zhì)模擬相對(duì)誤差小于30%，整體的模擬精度仍然能夠滿足實(shí)際需要。所以在一定精度范圍內(nèi)，本文基于插值地形的模型優(yōu)于實(shí)測(cè)地形模型的地方在于，基于地形插值的模型能夠用于缺少地形資料的地區(qū)。

3 結(jié) 論

本次研究通過(guò)二維圣維南方程、二維污染物遷移轉(zhuǎn)化方程，利用遙感數(shù)據(jù)獲取河道邊界及斷面地形，采用克里金插值法生成模型計(jì)算所需的地形標(biāo)準(zhǔn)文件，建立了水動(dòng)力水質(zhì)模型。經(jīng)過(guò)驗(yàn)證，水動(dòng)力模擬相對(duì)誤差不超過(guò)20%，水質(zhì)模擬相對(duì)誤差不超過(guò)30%。

將本次研究所構(gòu)建的基于地形插值模型與基于實(shí)測(cè)地形模型進(jìn)行對(duì)比，兩者使用相同的水文輸入條件，區(qū)別在于二者的地形數(shù)據(jù)不同。結(jié)果表明，本次研究所構(gòu)建的基于插值地形的模型水質(zhì)模擬最大相對(duì)誤差不超過(guò)30%，而基于實(shí)測(cè)地形模型的最大相對(duì)誤差不超過(guò)20%，實(shí)測(cè)地形模型的模擬精度更大，但是在一定的精度范圍內(nèi)，都可滿足不同實(shí)際的需要。下一步希望能夠在更多的研究區(qū)域使用本方法構(gòu)建水動(dòng)力水質(zhì)模型，從而多方面驗(yàn)證基于地形插值模型的精度。

由于我國(guó)河流數(shù)字化工作不足，大部分河流的地形資料還處于缺乏狀態(tài)。本文提出的基于遙感影像數(shù)據(jù)快速構(gòu)建水下地形進(jìn)行水動(dòng)力水質(zhì)模擬，雖然在模擬精度上不能有好的效果，但是在為無(wú)資料地區(qū)的突發(fā)水污染預(yù)警上還是有一定應(yīng)用實(shí)際意義，可以為污染物運(yùn)移趨勢(shì)提供預(yù)測(cè)，為污染事故防治提供指導(dǎo)。

特別值得說(shuō)明的是，本次研究所使用的遙感影像數(shù)據(jù)來(lái)源于Google Earth，所以模擬精度和Google Earth提供的數(shù)據(jù)精度密切相關(guān)。實(shí)踐者在使用時(shí)，請(qǐng)一定首先校核下數(shù)據(jù)精度是否在米級(jí)變化。我們相信，隨著技術(shù)的發(fā)展，會(huì)有更多其他遙感影像能夠提供精度更高的數(shù)據(jù)服務(wù)于水動(dòng)力水質(zhì)模型的構(gòu)建，水質(zhì)模型專業(yè)人士可以密切關(guān)注該領(lǐng)域的科技發(fā)展動(dòng)態(tài)，以便更好地服務(wù)于水環(huán)境保護(hù)領(lǐng)域。

□