譚 鵬,張?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ā)揮重要的作用。
對(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。
普通克里金插值法(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所示。
本研究采用二維圣維南方程和二維污染物遷移轉(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)
首先,以河道研究可能涉及的最大水位,可以是歷史高 洪水位或最高堤岸線,以此作為計(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ì)算的邊界條件。
研究河段為長(zhǎng)江干流九江段,如圖5所示,其上邊界斷面位于賽城閘上游4 km處,下邊界斷面為九江烏石磯糧庫(kù),長(zhǎng)約20 km,寬約1.5~2.5 km。研究河段內(nèi)有鶴問(wèn)湖污水處理廠和老鸛塘污水處理廠兩大污染源。
依據(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è)地形連貫合理。
本次研究以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.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ū)。
本次研究通過(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)域。
□