基于衛星影像數據快速構建水動力水質模型研究

譚 鵬，張詩豪，徐宏亮，劉宏寬，彭 虹，周文婷

(1.中電建生態環境集團有限公司，深圳 518102；2.武漢大學水利水電學院，武漢 430072)

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

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

1 研究方法

1.1 Google Earth獲取地形數據

對于沒有河道邊界及缺乏重點斷面實測資料的河段，本次研究采用的是通過利用Google Earth等衛星地圖工具，找到目的河道并勾取河道邊界和重點斷面，得到一系列邊界散點和斷面散點的經緯度坐標(WGS84坐標系)及對應的高程信息，并將各散點的經緯度坐標轉換為模型計算所需的平面坐標。

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

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

X=(N+H)cosBcosL

(1)

Y=(N+H)cosBsinL

(2)

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

(3)

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

(4)

(5)

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

1.2 普通克里金插值

普通克里金插值法(Ordinary Kriging)是克里金插值法的一種，是一種依據協方差函數對隨機過程進行插值的回歸算法，該方法認為可依據區域上的一系列點X1,X2,…,Xn的觀測值，采用一個線性組合來估計區域化變量在X0處的值，即：

(6)

滿足插值的無偏性和最優性，即要求：

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

(7)

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

(8)

(9)

轉化為λi求解的矩陣形式為：

(10)

其中Cij可通過變差函數γ(h)來確定，二維情況下，空間點P在X軸和Y軸上均有變化，可分別對x和y方向進行一維情況下的變差函數計算，即：

(11)

(12)

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

(13)

式中：a為變程，即變量之間具有相關性的區域范圍；h為滯后距；c為基臺值，即變差函數在h大于變程時的值，為塊金值c0和拱高cc之和。

塊金值c0為變量純隨機性的部分，即在很短的距離內有較大的空間變異性，無論h多小，兩個隨機變量都不相關；拱高cc為觀測值的變異幅度范圍。變差函數、變程、滯后距、基臺值、拱高、塊金值之間的關系如圖4所示。

1.3 二維水動力水質模型

本研究采用二維圣維南方程和二維污染物遷移轉化方程模擬寬深比較大、斷面變化較大的大型河道水動力學特征和水質遷移轉化規律，即：

(1)二維圣維南方程包含連續方程和動量方程。

連續方程：

(14)

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

動量方程：

(15)

(16)

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

邊界條件：

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

Q=Qin(t)

(17)

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

z=zout(t)

(18)

初值條件：

φ=φin，t=0

(19)

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

考慮對流、擴散、降解等項，二維水質模型基本方程如下：

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

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

(20)

式中：ci為第i種物質的斷面平均濃度，mg/L；Ex和Ey為河道的縱向和橫向擴散系數，m2/s；SLi、SBi和SKi為第i種物質側流入匯負Z荷、河道中源漏項，以及因水體水生態系統動力轉換項，m·mg/(L·s)。

邊界條件：

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

C=Cin(t)

(21)

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

(22)

初值條件：

Ci=Cin，t=0

(23)

1.4 計算方法

首先，以河道研究可能涉及的最大水位，可以是歷史高 洪水位或最高堤岸線，以此作為計算的外邊界，通過網格生成程序生成模型計算的最大范圍的計算網格。然后，依據圣維南方程，通過初始和邊界條件，利用來流和下游水位及側流入匯過程條件，根據橫斷面的水位變化，通過模型迭代計算，自動識別，判斷出網格的干濕狀態后，濕網格參與計算，干網格不參與計算，依次即可完成水動力水質模型中相應每個濕網格上的水位、流速和水質濃度的計算。

其中谷歌地球提取出來的高程，在校正插值后，構建出研究區域的水下地形，作為模型計算的邊界條件。

2 研究區域

2.1 區域概況

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

2.2 河道地形生成

依據上文所講到的地形獲取方法，根據九江河段的實際位置，使用Google Earth選取河道邊界并沿河道均勻選取30個斷面，得到河道邊界和斷面散點坐標，如圖6所示。

利用獲取的河道邊界，繪制出平面二維正交網格，網格對應實際地形尺寸約為30 m×30 m，如圖7所示。

依據河道邊界散點以及選取的斷面散點，結合網格結點坐標，利用普通克里金法插值計算每個網格結點處的高程信息，生成模型計算所需的地形標準文件。

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

2.3 模型計算參數

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

表1 水動力水質模型參數率定范圍Tab.1 Calibration range of hydrodynamic water quality model parameters

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

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

2.4 結果對比

2.4.1 水動力對比

本次研究區域內，共有CS1、九江站、CS2 3個水文測驗斷面，如圖9所示。其驗證數據為九江水文局在2008年9月6日的水文測驗實測水文資料。

本研究選取這3個斷面的實測水位、流速進行模型驗證。其中，CS1、九江站、CS2 3個水文測驗斷面平均水位驗證結果及誤差如表3所示，流速驗證結果及誤差如表4所示。

表3 兩種地形條件下水位模擬結果誤差統計表Tab.3 Statistics of water level simulation results under two terrain conditions

由表3和表4可以看出，基于插值地形的模型在絕大多實測點的相對誤差值都比基于實測地形的模型要大，其中基于實測地形模型的相對誤差都在10%以下，而本次研究所構建的基于插值地形的模型相對誤差均在20%以下，兩種地形所構建的模型都能夠較好的模擬河道水動力規律。

表4 兩種地形條件下模型流速模擬結果誤差統計表Tab.4 Error statistics of simulation results of model velocity under two terrain conditions

2.4.2 水質對比

模型水質驗證根據九江水文局2015年9月2日實測水質數據進行驗證，選取研究河段內的抗洪廣場、河西水廠、第三水廠(河東水廠)、烏石磯糧庫為水質測驗斷面，如圖10所示。

水質驗證輸入條件如表5所示。

表5 水質驗證輸入條件Tab.5 Input conditions for water quality verification

基于兩種地形的模型在抗洪廣場、河西水廠、第三水廠、烏石磯斷面4個水質監測斷面處的COD及NH3-N斷面平均濃度模擬結果及誤差如表6、7所示。

表6 COD濃度模擬結果誤差統計表Tab.6 Error statistics of COD concentration simulation results

表7 NH3-N濃度模擬結果誤差統計表Tab.7 Error statistics of NH3-N concentration simulation results

從表6和表7可知，本次研究基于插值地形模型在水質模擬方面，相對誤差普遍大于基于實測地形的模型，插值地形的模型結果最大相對誤差出現在河西水廠的NH3-N模擬上，達到了23.71%，小于30%；實測地形模型的最大相對誤差同樣出現在河西水廠的NH3-N模擬上，達到了16.88%，小于20%。

綜合對比水動力模擬精度以及水質模擬精度，本次研究基于插值地形的模型模擬精度不如實測地形，但是其在驗證斷面處，水位模擬誤差小于0.2 m，流速模擬相對誤差小于20%，水質模擬相對誤差小于30%，整體的模擬精度仍然能夠滿足實際需要。所以在一定精度范圍內，本文基于插值地形的模型優于實測地形模型的地方在于，基于地形插值的模型能夠用于缺少地形資料的地區。

3 結 論

本次研究通過二維圣維南方程、二維污染物遷移轉化方程，利用遙感數據獲取河道邊界及斷面地形，采用克里金插值法生成模型計算所需的地形標準文件，建立了水動力水質模型。經過驗證，水動力模擬相對誤差不超過20%，水質模擬相對誤差不超過30%。

將本次研究所構建的基于地形插值模型與基于實測地形模型進行對比，兩者使用相同的水文輸入條件，區別在于二者的地形數據不同。結果表明，本次研究所構建的基于插值地形的模型水質模擬最大相對誤差不超過30%，而基于實測地形模型的最大相對誤差不超過20%，實測地形模型的模擬精度更大，但是在一定的精度范圍內，都可滿足不同實際的需要。下一步希望能夠在更多的研究區域使用本方法構建水動力水質模型，從而多方面驗證基于地形插值模型的精度。

由于我國河流數字化工作不足，大部分河流的地形資料還處于缺乏狀態。本文提出的基于遙感影像數據快速構建水下地形進行水動力水質模擬，雖然在模擬精度上不能有好的效果，但是在為無資料地區的突發水污染預警上還是有一定應用實際意義，可以為污染物運移趨勢提供預測，為污染事故防治提供指導。

特別值得說明的是，本次研究所使用的遙感影像數據來源于Google Earth，所以模擬精度和Google Earth提供的數據精度密切相關。實踐者在使用時，請一定首先校核下數據精度是否在米級變化。我們相信，隨著技術的發展，會有更多其他遙感影像能夠提供精度更高的數據服務于水動力水質模型的構建，水質模型專業人士可以密切關注該領域的科技發展動態，以便更好地服務于水環境保護領域。

□