基于四維變分同化法的巢湖流域南淝河水質(zhì)模擬*

賴錫軍，何國(guó)建

(1：中國(guó)科學(xué)院南京地理與湖泊研究所，流域地理學(xué)重點(diǎn)實(shí)驗(yàn)室，南京 210008) (2：清華大學(xué)水利水電工程系，水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084)

采用數(shù)學(xué)模型模擬河流水質(zhì)過(guò)程是診斷和預(yù)測(cè)河流水環(huán)境問(wèn)題的重要手段[1]. 當(dāng)前，河流水質(zhì)模型已被廣泛應(yīng)用于面向流域水環(huán)境系統(tǒng)的科學(xué)研究、工程規(guī)劃和管理決策等方面，已成為流域水環(huán)境管理不可或缺的手段[2]. 在應(yīng)用傳統(tǒng)的數(shù)學(xué)模型模擬和預(yù)測(cè)河流水環(huán)境過(guò)程時(shí)，需要完整的邊界條件和初始條件信息以及可靠的模型參數(shù). 然而，由于實(shí)際河流的復(fù)雜性，開(kāi)展河流水質(zhì)模擬分析時(shí)常常難以全面收集到河流水質(zhì)初始條件和邊界條件等定解資料. 特別是入河污染負(fù)荷，因其由點(diǎn)面源污染綜合產(chǎn)生，具有過(guò)程復(fù)雜、不確定性大等特點(diǎn)，是影響水質(zhì)模擬精度的重要不確定性源[3]. 當(dāng)前，我國(guó)對(duì)于大型工業(yè)企業(yè)、污水處理廠等大型排口有較好的監(jiān)測(cè)，可較準(zhǔn)確地獲取其排放過(guò)程，但是對(duì)于一些小型排口、地表漫流以及地下水滲入等，其入河規(guī)模和位置等情況都難以準(zhǔn)確估計(jì). 這些缺失的或者具有高度不確定性的污染負(fù)荷數(shù)據(jù)使常規(guī)的河流水質(zhì)模型運(yùn)行缺少了必要的準(zhǔn)確可靠的輸入數(shù)據(jù)，導(dǎo)致常規(guī)方法不再有效. 考慮到當(dāng)前河流關(guān)鍵斷面水質(zhì)濃度有相對(duì)完整的觀測(cè)資料，值得思考如何更好地去利用這些數(shù)據(jù)，探求河流水質(zhì)模擬的新方法，為污染負(fù)荷缺失條件下河流水質(zhì)模擬提供新的可能.

基于偏微分方程最優(yōu)化控制理論的變化同化方法(又稱“四維變分同化”，4D-Var)[4]可以通過(guò)同化不同類型觀測(cè)數(shù)據(jù)，識(shí)別模型模擬所需的時(shí)空變化參數(shù)，改進(jìn)模型動(dòng)力過(guò)程模擬預(yù)測(cè)能力. 其在大氣科學(xué)、海洋科學(xué)、水文學(xué)和水力學(xué)等不同門類學(xué)科已有成功的應(yīng)用[5-9]. 不過(guò)，采用變分同化方法開(kāi)展河流水質(zhì)模擬的報(bào)道仍很少，Cho等[10]在其地表水水質(zhì)模擬數(shù)據(jù)同化的回顧論文中提到了2篇利用三維變分同化方法(3D-Var)模擬一維渠道的水質(zhì)輸運(yùn)過(guò)程[11]和水庫(kù)的水質(zhì)模擬[12]. 鑒于四維變分同化方法既具備同化觀測(cè)數(shù)據(jù)，又可保持物理動(dòng)力規(guī)律約束的能力[13-14]，值得加強(qiáng)其在河流水質(zhì)模擬中的應(yīng)用研究.

南淝河位于巢湖流域西北部合肥市境內(nèi)，為城市河流，水環(huán)境問(wèn)題突出. 南淝河是巢湖的主要入湖河流，其水質(zhì)達(dá)標(biāo)是巢湖流域水污染治理的重點(diǎn)任務(wù)[15]. 掌握其水環(huán)境變化過(guò)程及其入河負(fù)荷對(duì)水質(zhì)的影響是開(kāi)展水環(huán)境治理的基礎(chǔ). 但是，與國(guó)內(nèi)城市河流類似，其存在污染來(lái)源難以全面排查清楚等問(wèn)題[16]. 按照常規(guī)的傳統(tǒng)模型模擬方法，需要提前估算污染入河負(fù)荷，才能實(shí)現(xiàn)水質(zhì)過(guò)程的模擬. 對(duì)于這個(gè)入河污染負(fù)荷存在很高不確定性的河流，若沒(méi)有可靠的污染入河數(shù)據(jù)，將不可能準(zhǔn)確模擬其水質(zhì)動(dòng)態(tài). 本文以南淝河干流董鋪水庫(kù)以下城區(qū)河流水質(zhì)模擬為例，探究基于變分同化法模擬水質(zhì)動(dòng)力過(guò)程的性能，為入河污染負(fù)荷缺失或具高不確定性條件下河流水質(zhì)模擬提供可選方法. 研究顯示，方法可在提升河流水質(zhì)模擬精度的同時(shí)識(shí)別不確知的入河污染負(fù)荷，為巢湖流域等環(huán)境條件復(fù)雜的流域系統(tǒng)中污染定量解析、水質(zhì)預(yù)測(cè)預(yù)警及污染管控等提供重要支持.

1 研究區(qū)域

南淝河位于合肥市境內(nèi)，是巢湖重要的入河支流，也是典型城市河流. 南淝河干流河長(zhǎng)約70 km，發(fā)源于江淮分水嶺大潛山余脈，從西北向東南流入董鋪水庫(kù)、流經(jīng)合肥市主城區(qū)后于入湖口施口處匯入巢湖. 南淝河自董鋪水庫(kù)大壩至入湖口施口全長(zhǎng)42.0 km，其間依次有四里河、板橋河、二十埠河、店埠河和長(zhǎng)樂(lè)河等支流匯入.

南淝河所在流域多年平均降水量為1022 mm. 盡管降水較為充沛，但是由于上游董鋪水庫(kù)為合肥市的主要飲用水水源，壩下游河段仍缺乏充足的清潔水源補(bǔ)給. 南淝河水源主要為污水處理廠處理后的生產(chǎn)和生活用水的尾水，以及區(qū)域內(nèi)的城市和農(nóng)村的地表徑流等. 受區(qū)域密集人口及高速發(fā)展的社會(huì)經(jīng)濟(jì)影響，南淝河水質(zhì)近年多處于Ⅴ類和劣V類水平，水環(huán)境問(wèn)題突出. 目前，干流計(jì)算區(qū)域內(nèi)有望塘、清溪、王小郢、小倉(cāng)房4座污水處理廠尾水入河；店埠河有撮鎮(zhèn)和聯(lián)禧污水處理廠尾水入河(圖1). 其中干流望塘、清溪、王小郢和小倉(cāng)房污水處理規(guī)?？偭繛?10萬(wàn)m3/d. 從2019年運(yùn)行數(shù)據(jù)可知，各個(gè)污水廠排水基本穩(wěn)定，達(dá)到了準(zhǔn)Ⅳ類水質(zhì)要求.

圖1 研究區(qū)南淝河水系及關(guān)鍵斷面分布Fig.1 The Nanfei River and the main gauging points

2 模型方法與構(gòu)建

2.1 河流水質(zhì)變分同化方法

南淝河水質(zhì)沿程動(dòng)態(tài)變化過(guò)程可采用具有反應(yīng)源項(xiàng)的一維物質(zhì)輸運(yùn)方程來(lái)表達(dá). 為了簡(jiǎn)便，僅考慮水體中污染物質(zhì)沿河縱向的隨流輸移運(yùn)動(dòng)以及沿程的生化降解過(guò)程，具體為：

(1)

式中，C為某水質(zhì)指標(biāo)的濃度，mg/L；t為時(shí)間，s；D為縱向離散系數(shù)，m2/s；U為斷面平均流速，m/s；K為某水質(zhì)指標(biāo)的一階衰減系數(shù)，s-1；S為某水質(zhì)指標(biāo)的外部源匯項(xiàng)，g/(m3·s)；x為空間坐標(biāo).

水質(zhì)數(shù)據(jù)同化的目標(biāo)主要是通過(guò)實(shí)際觀測(cè)水質(zhì)濃度等數(shù)據(jù)和水質(zhì)動(dòng)態(tài)模型相互融合以達(dá)到提高模型模擬和預(yù)測(cè)真實(shí)物理過(guò)程精度. 也就是在水質(zhì)控制方程約束下，尋找最優(yōu)的控制變量P(模型初始條件、邊界條件和過(guò)程參數(shù))，使表征模擬結(jié)果與觀測(cè)差異的代價(jià)函數(shù)J為最小，即：

(2)

式中，T為積分(模擬)時(shí)段；H為觀測(cè)算子；?Ω表示計(jì)算域的邊界.

如果能計(jì)算出代價(jià)函數(shù)J關(guān)于控制變量P的梯度. 即可用梯度下降類方法求解優(yōu)化求解代價(jià)函數(shù)極小值. 變分同化法通過(guò)求解伴隨方程來(lái)快速計(jì)算代價(jià)函數(shù)關(guān)于控制變量的梯度[14]. 方程(1)的伴隨方程可表達(dá)為：

(3)

原方程(1)，伴隨方程(3)和相應(yīng)的初始、邊界條件構(gòu)成了雙向的積分系統(tǒng). 其即為一個(gè)關(guān)于控制變量P的最優(yōu)化控制系統(tǒng)[17]. 利用最優(yōu)化算法可實(shí)現(xiàn)該最優(yōu)控制系統(tǒng)的求解，即模擬得到與觀測(cè)資料最為接近的水質(zhì)過(guò)程. 基于該最優(yōu)化系統(tǒng)建立的模型即為水質(zhì)變分同化模型. 本文建模采用了適用于大規(guī)模問(wèn)題計(jì)算的擬牛頓算法L-BFGS[18]作為最優(yōu)化算法，同化模型整體架構(gòu)可參考文獻(xiàn)[17]，這里不再贅述. 模型采用Fortran語(yǔ)言編碼.

2.2 南淝河模型構(gòu)建

根據(jù)南淝河的實(shí)際可用資料情況，模型計(jì)算范圍設(shè)定為董鋪水庫(kù)以下5.3 km至入湖口施口(圖1)，模擬干流段總長(zhǎng)34.2 km. 具體包括南淝河干流(從望塘污水廠下游監(jiān)測(cè)站開(kāi)始至施口斷面止)、四里河入河口段(臨泉路)、板橋河入河口段(北二環(huán)路橋)、二十埠河入河段(合馬路)、店埠河入口段(撮鎮(zhèn)). 河道模擬計(jì)算分為14個(gè)河段，干流內(nèi)部的分段節(jié)點(diǎn)分別為四里河入河口、板橋河入河口、長(zhǎng)江東街、南一環(huán)、銅陵路、當(dāng)涂路、范拐、二十埠河入口以及店埠河入口. 各河段又劃分成計(jì)算斷面，各斷面間距根據(jù)河道形態(tài)取800～1200 m不等.

南淝河的水質(zhì)計(jì)算邊界主要包括干支流入流的水質(zhì)濃度過(guò)程(第一類邊界條件)、南淝河入湖口的自由輸出邊界條件(第二類邊界條件)以及側(cè)向污水處理廠和長(zhǎng)樂(lè)河等支流的入流負(fù)荷邊界條件. 入流濃度過(guò)程包括干流上游、四里河臨泉路、板橋河北二環(huán)路橋、二十埠河合馬路和店埠河撮鎮(zhèn)共5個(gè)邊界水質(zhì)指標(biāo)濃度過(guò)程.

入河負(fù)荷是水質(zhì)計(jì)算的主要影響要素. 南淝河污染包括污水處理廠等排口以及地表漫流. 已知的污水廠排口數(shù)據(jù)作為點(diǎn)源考慮，包括計(jì)算區(qū)域內(nèi)直接匯入干流的清溪、望塘、王小郢、小倉(cāng)房以及匯入店埠河的聯(lián)熹和撮鎮(zhèn)污水處理廠尾水. 對(duì)于一些小型排口及地表漫流，其入河負(fù)荷情況難以掌握. 在本次計(jì)算中將這些未知污染負(fù)荷計(jì)入非點(diǎn)源處理，作為模型的控制變量，通過(guò)變分同化模型進(jìn)行識(shí)別計(jì)算，解決了傳統(tǒng)模型若數(shù)據(jù)缺失則難以可靠模擬的問(wèn)題.

模型的其他計(jì)算參數(shù)，如河流非恒定的水動(dòng)力參數(shù)由水動(dòng)力模型計(jì)算提供；本次計(jì)算氨氮和總磷一階衰減系數(shù)[19]分別取0.06和0.04 d-1.

2.3 數(shù)據(jù)選用

本文以氨氮和總磷作為代表性水質(zhì)指標(biāo)，模擬時(shí)段為2019年10-12月. 南淝河水質(zhì)模擬的邊界條件和初始條件等數(shù)據(jù)根據(jù)自動(dòng)監(jiān)測(cè)站采集的逐日平均氨氮和總磷數(shù)據(jù)給定. 污水處理廠尾水排放數(shù)據(jù)也基于逐日監(jiān)測(cè)的過(guò)程數(shù)據(jù). 模型同化觀測(cè)資料來(lái)自于南淝河干流自動(dòng)監(jiān)測(cè)站，本次選用了方興大道、板橋碼頭、范拐、當(dāng)涂路、銅陵路、南一環(huán)、長(zhǎng)江東街、宿州路橋、亳州路橋、四里河濉溪路橋、板橋河鴛鴦路橋、二十埠河入南淝河口、店埠河二電廠(圖1)等站的自動(dòng)監(jiān)測(cè)數(shù)據(jù)，時(shí)間頻次為d.

3 結(jié)果分析與討論

3.1 模擬結(jié)果

采用模型分別計(jì)算了分散點(diǎn)源和面源入河污染負(fù)荷數(shù)據(jù)缺失條件下南淝河水質(zhì)過(guò)程，以及將這些非點(diǎn)源負(fù)荷作為控制變量模擬的結(jié)果. 由于實(shí)際問(wèn)題非線性強(qiáng)，迭代收斂性能受制約，本次南淝河模擬經(jīng)過(guò)5000次迭代，代價(jià)函數(shù)相對(duì)值從1降低至0.416，降幅為59.4%. 以此結(jié)果作為采用變分法得到的計(jì)算成果進(jìn)行了分析. 圖2和圖3展示了分散點(diǎn)源和面源入河污染負(fù)荷數(shù)據(jù)缺失條件下初始計(jì)算結(jié)果和采用變分法計(jì)算結(jié)果與實(shí)際過(guò)程的比較.

圖2 變分同化前后各站點(diǎn)氨氮過(guò)程線與實(shí)測(cè)值的比較Fig.2 Comparison of measured ammonia nitrogen concentration hydrograph with the modeled ones before and after variational assimilation

圖3 變分同化前后各站點(diǎn)總磷過(guò)程線與實(shí)測(cè)值的比較Fig.3 Comparison of measured total phorsphorus concentration hydrograph with the modeled ones before and after variational assimilation

從各監(jiān)測(cè)斷面水質(zhì)過(guò)程來(lái)看，若分散點(diǎn)源以及面源的入河數(shù)據(jù)缺失，計(jì)算得到各斷面水質(zhì)過(guò)程在長(zhǎng)期無(wú)降雨時(shí)段(11月)，當(dāng)涂路以上站點(diǎn)的水質(zhì)過(guò)程相對(duì)比較一致. 不過(guò)，下游段，特別是板橋碼頭和范拐斷面的計(jì)算值較實(shí)測(cè)值偏低，說(shuō)明在該河段內(nèi)旱季也有污染匯入；而方興大道低估又有減少，這可能和南淝河會(huì)出現(xiàn)水流倒灌入河的情況有關(guān). 在降雨影響的時(shí)段，計(jì)算結(jié)果偏差較大. 統(tǒng)計(jì)表明(表1)，各斷面模擬結(jié)果的均方根誤差(RMSE)明顯偏大，納什效率系數(shù)(NSE)除個(gè)別斷面?zhèn)€別指標(biāo)外，多數(shù)小于0. 這也說(shuō)明，南淝河水質(zhì)模擬僅考慮污水處理廠等尾水排放還不夠. 受降雨影響的時(shí)段，南淝河沿程仍有大量雨污水匯入南淝河，造成南淝河的污染物濃度升高，水質(zhì)變差，需對(duì)這些分散式的污染負(fù)荷加以考慮.

表1 變分同化前后各站點(diǎn)氨氮和總磷模擬計(jì)算精度比較

采用變分法同化觀測(cè)數(shù)據(jù)后，模擬得到的各斷面水質(zhì)過(guò)程計(jì)算結(jié)果得以大幅校正，捕捉到了各個(gè)污染指標(biāo)波動(dòng)的峰值，較好地模擬了南淝河實(shí)際的水質(zhì)變化過(guò)程. 雖然存在計(jì)算偏高和偏低等問(wèn)題，特別是范拐斷面的氨氮模擬結(jié)果偏低較為明顯，可能由于實(shí)際問(wèn)題的復(fù)雜性限制，使模型并沒(méi)有收斂到理想狀態(tài). 但是，從總體來(lái)看，計(jì)算結(jié)果均方根誤差有明顯下降. 氨氮RMSE從各斷面平均的1.03 mg/L降至0.69 mg/L，NSE則從平均-0.66上升到0.16. 總磷R(shí)MSE從各斷面平均的0.07 mg/L降至0.06 mg/L，NSE則從平均-1.50 上升到-0.21. 從各斷面的各自變化來(lái)看，干流長(zhǎng)江東街到當(dāng)涂路河段的水質(zhì)模擬精度改進(jìn)更為明顯. 氨氮和總磷的NSE均從小于0分別提高至0.62和0.50. 下游施口、方興大道等斷面模擬精度較差，可能與巢湖水體出現(xiàn)倒流等因素有關(guān).

3.2 負(fù)荷識(shí)別結(jié)果

本次模擬以入河非點(diǎn)源負(fù)荷作為控制變量, 采用變分同化模型模擬得到南淝河水質(zhì)過(guò)程的同時(shí), 也識(shí)別得到了南淝河的入河非點(diǎn)源負(fù)荷過(guò)程. 各河段在模擬時(shí)段內(nèi)(10-12月)入河負(fù)荷總通量如圖4所示. 南淝河沿程各河段均有污染輸入. 其中，氨氮負(fù)荷在當(dāng)涂路-范拐河段、二十埠河入南淝河段、四里河口-板橋河口河段、板橋河入南淝河河段以及南淝河入湖段均較高；總磷負(fù)荷較高的河段則為南淝河入湖段、二十埠河入南肥河段以及店埠河入南肥河段等. 通過(guò)同化不僅可得到各河段的負(fù)荷總通量數(shù)據(jù)，而且直接得到了入河污染負(fù)荷通量隨時(shí)間的變化過(guò)程. 這里選擇南淝河入湖段以及當(dāng)涂路-范拐河段的入河過(guò)程為例加以說(shuō)明模型的性能. 從圖5中可以看出，污染負(fù)荷隨降雨過(guò)程的變化比較明顯，對(duì)應(yīng)10月5日、11月24-30日以及12月16-25日每場(chǎng)降雨，各河段的氮、磷均可形成不同量級(jí)的負(fù)荷峰值. 說(shuō)明降雨引起的雨污溢流問(wèn)題對(duì)南淝河水質(zhì)有重要影響. 此外，本次計(jì)算還識(shí)別出了當(dāng)涂路-范拐段河段10月末存在一個(gè)很高的氨氮入河峰值，為647.45 g/s. 雖然該負(fù)荷未知來(lái)源，但是根據(jù)范拐及板橋碼頭斷面在10月底出現(xiàn)的氨氮濃度峰值，可以確認(rèn)這是合理的. 這也從側(cè)面說(shuō)明了該模型對(duì)未知的污染來(lái)源有很強(qiáng)的識(shí)別能力，可在今后應(yīng)用于河流污染來(lái)源的定量追溯分析. 對(duì)于巢湖等流域入湖污染控制與河流水質(zhì)預(yù)測(cè)預(yù)警具有重要價(jià)值.

圖4 南淝河各河段識(shí)別的污染負(fù)荷Fig.4 Identified pollutant loads in each reach of the Nanfei River

圖5 南淝河入湖段及當(dāng)涂路-范拐河段的入河負(fù)荷過(guò)程Fig.5 Hydrograph of pollutant loads in the two typical reaches and the corresponding rainfall during this period

4 結(jié)論

以巢湖流域南淝河干流為研究對(duì)象，基于變分同化方法構(gòu)建了南淝河水質(zhì)模型，模擬了2019年10-12月間南淝河干流及支流入河段的水質(zhì)時(shí)空動(dòng)態(tài)變化. 與實(shí)際觀測(cè)數(shù)據(jù)以及基于已知污染負(fù)荷數(shù)據(jù)代入常規(guī)模型模擬結(jié)果的對(duì)比分析表明：

1)基于變分同化方法建立的南淝河水質(zhì)模型可模擬入河污染負(fù)荷數(shù)據(jù)缺失條件下南淝河水質(zhì)過(guò)程，成功捕捉水質(zhì)波動(dòng)變化的不同峰值，水質(zhì)濃度模擬精度明顯改進(jìn). 該方法為缺資料地區(qū)河流模擬提供了可能的手段.

2)模型在計(jì)算水質(zhì)過(guò)程的同時(shí)，也得到了南淝河各河段的入河負(fù)荷動(dòng)態(tài)變化過(guò)程，為河流水質(zhì)污染溯源及成因分析提供了定量數(shù)據(jù)支持. 模擬結(jié)果表明，本文提出的方法對(duì)入河污染敏感，可用于關(guān)鍵河段污染定量溯源，可支撐巢湖等流域復(fù)雜河流的入湖污染控制.

3)基于變分同化法的水質(zhì)模擬通過(guò)同化河道內(nèi)水質(zhì)指標(biāo)濃度觀測(cè)數(shù)據(jù)，代替常規(guī)模型需要以入河污染負(fù)荷數(shù)據(jù)作為輸入條件的應(yīng)用前提. 亦即以河道內(nèi)的濃度觀測(cè)數(shù)據(jù)作為輸入條件代替污染入河負(fù)荷作為輸入條件，避開(kāi)了入河負(fù)荷準(zhǔn)確估算的難點(diǎn)，使得在復(fù)雜環(huán)境條件下河流水質(zhì)模擬更為可靠. 南淝河的水質(zhì)模擬研究展現(xiàn)了該模型在復(fù)雜環(huán)境條件下的應(yīng)用潛力，對(duì)流域河流水質(zhì)預(yù)測(cè)預(yù)警及污染管控具有重要意義.