孫曉卓，曾獻(xiàn)奎，吳吉春，孫媛媛

（南京大學(xué)地球科學(xué)與工程學(xué)院/表生地球化學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023）

地下水?dāng)?shù)值模型能夠定量刻畫地下水流和污染物的時(shí)空分布特征，已被廣泛應(yīng)用于地下水資源的科學(xué)管理和保護(hù)[1?5]。由于實(shí)際地下水系統(tǒng)十分復(fù)雜，限于認(rèn)知不足，通常需對(duì)其進(jìn)行簡化，這將導(dǎo)致建立的水文地質(zhì)概念模型與實(shí)際水文地質(zhì)條件之間存在差異，即模型存在結(jié)構(gòu)誤差，如含水層空間分布的刻畫偏差、化學(xué)反應(yīng)動(dòng)力學(xué)模型存在偏差等。使用帶有結(jié)構(gòu)誤差的模型進(jìn)行預(yù)測時(shí)，即使進(jìn)行參數(shù)識(shí)別，預(yù)測結(jié)果也可能存在偏差[6?8]。因此，為提高模擬結(jié)果的可靠性，需要定量刻畫模型的結(jié)構(gòu)誤差，從而可為地下水資源管理和決策提供更可靠的依據(jù)。

近十年來，眾多研究者開始注重對(duì)模型結(jié)構(gòu)不確定性的研究[9?12]。貝葉斯模型平均（BMA）和基于數(shù)據(jù)驅(qū)動(dòng)的結(jié)構(gòu)誤差統(tǒng)計(jì)學(xué)習(xí)方法是目前常用的處理方法。BMA 在貝葉斯框架下建立一組可行的概念模型，根據(jù)每個(gè)模型的表現(xiàn)賦予相應(yīng)的權(quán)重，從而對(duì)這些模型的預(yù)測結(jié)果進(jìn)行加權(quán)平均[13?18]。研究結(jié)果表明，對(duì)比于其他多模型方法或單一模型法，采用BMA方法能得到更準(zhǔn)確和可靠的預(yù)測結(jié)果[19?20]。盡管如此，BMA 在實(shí)際應(yīng)用過程中存在局限性，如模型空間（概念模型集合）的定義及先驗(yàn)權(quán)重具有主觀性、模型間相關(guān)性未能準(zhǔn)確量化等，從而影響B(tài)MA 的預(yù)測效果[20?21]。

數(shù)據(jù)驅(qū)動(dòng)法（DDM）通過某種統(tǒng)計(jì)模型刻畫模型結(jié)構(gòu)誤差，如支持向量機(jī)回歸[22?23]、人工神經(jīng)網(wǎng)絡(luò)[24]、隨機(jī)森林[25]、高斯過程回歸[26?27]等機(jī)器學(xué)習(xí)方法。其中，高斯過程回歸（GPR）是一種基于貝葉斯理論的監(jiān)督學(xué)習(xí)算法，被廣泛用于模型結(jié)構(gòu)誤差的學(xué)習(xí)。Xu等[26?27]利用GPR 刻畫多個(gè)地下水流模型的結(jié)構(gòu)誤差，通過馬爾科夫鏈蒙特卡洛模擬（MCMC）識(shí)別地下水模型參數(shù)（物理參數(shù)）和GPR 參數(shù)（超參）。結(jié)果表明，GPR 可以刻畫模型結(jié)構(gòu)誤差的時(shí)空相關(guān)性，顯著提高地下水模型的預(yù)測能力。Pan 等[28]通過GPR 刻畫地下水中重質(zhì)非水相液體（DNAPL）運(yùn)移模型的結(jié)構(gòu)誤差，并進(jìn)行針對(duì)地下水DNAPL 污染的人類健康風(fēng)險(xiǎn)評(píng)價(jià)。結(jié)果表明與不考慮模型結(jié)構(gòu)誤差相比，通過GPR 刻畫DNAPL 運(yùn)移模型結(jié)構(gòu)誤差，能夠在關(guān)鍵區(qū)域獲得更準(zhǔn)確的人體健康風(fēng)險(xiǎn)評(píng)價(jià)。

盡管如此，當(dāng)前基于GPR 進(jìn)行模型結(jié)構(gòu)誤差刻畫的研究中，通常假設(shè)物理參數(shù)和超參相互獨(dú)立，且對(duì)這些參數(shù)進(jìn)行聯(lián)立識(shí)別。然而，在實(shí)際情況下物理參數(shù)和超參的獨(dú)立性假設(shè)不夠合理。同時(shí)，超參沒有物理意義，進(jìn)行物理參數(shù)和超參聯(lián)立識(shí)別時(shí)，可能導(dǎo)致參數(shù)的過度識(shí)別，影響模型預(yù)測效果[29?30]。

本文提出兩步識(shí)別DDM 方法進(jìn)行參數(shù)識(shí)別和模型結(jié)構(gòu)誤差定量刻畫，克服了傳統(tǒng)聯(lián)立識(shí)別DDM 方法中的過度擬合問題，并通過一個(gè)地下水溶質(zhì)運(yùn)移案例和一個(gè)地下水水流案例，對(duì)提出方法的可靠性和優(yōu)越性進(jìn)行了驗(yàn)證和分析。

1 研究方法

1.1 高斯過程回歸（GPR）

地下水模型中觀測值f可以表示為：

式中：M—地下水模型；

θ—地下水模型參數(shù)（物理參數(shù)）；

b—模型結(jié)構(gòu)誤差；

x—模型輸入；

φ—超參；

η—隨機(jī)測量誤差。

在GPR 中，假定結(jié)構(gòu)誤差b（x,φ）在空間各點(diǎn)服從均值為μ，協(xié)方差矩陣為C的聯(lián)合多元高斯分布[31?32]。本文均值函數(shù)μ為0，協(xié)方差矩陣C采用平方指數(shù)型：

式中：λ—特征長度；

σ2—控制b（x,φ）的邊緣方差；

n—觀測數(shù)據(jù)的個(gè)數(shù)；

I—n×n階的單位矩陣。

測量誤差η的先驗(yàn)為獨(dú)立同分布。物理參數(shù)θ和超參被認(rèn)為是相互獨(dú)立并服從先驗(yàn)分布p（θ,φ）=p（θ）p（φ）。根據(jù)貝葉斯原理，推導(dǎo)得到θ與φ的后驗(yàn)分布為：

式中：p（θ）—θ的先驗(yàn)函數(shù)；

p（φ）—φ的先驗(yàn)函數(shù)。

1.2 基于聯(lián)立識(shí)別的數(shù)據(jù)驅(qū)動(dòng)方法

傳統(tǒng)的數(shù)據(jù)驅(qū)動(dòng)方法（DDM）中，通常假設(shè)物理參數(shù)θ和超參φ獨(dú)立分布，通過MCMC 模擬對(duì)其進(jìn)行聯(lián)立識(shí)別[26,28]。本次研究中采用改進(jìn)的差分進(jìn)化自適應(yīng)算法DREAMzs（DiffeRential Evolution Adaptive Metropolis）[33?34]進(jìn)行MCMC 模擬，基于聯(lián)立識(shí)別的DDM 基本步驟如下：

（1）定義物理參數(shù)θ和超參φ的先驗(yàn)分布；

（2）通過MCMC 模擬獲得物理參數(shù)θ和超參φ的后驗(yàn)樣本；

（3）對(duì)于新的模型輸入x?，輸入模型獲得輸出M?，根據(jù)式（2）計(jì)算協(xié)方差矩陣C、C?、C??；

（4）根據(jù)式（5）—（7），使用多元正態(tài)抽樣生成b?、η?；

（5）通過式（8）計(jì)算預(yù)測值f*。

1.3 基于兩步識(shí)別的數(shù)據(jù)驅(qū)動(dòng)方法

假設(shè)物理參數(shù)θ和超參φ不獨(dú)立，將其分開進(jìn)行兩步識(shí)別。首先，在φ的全部先驗(yàn)分布空間內(nèi)識(shí)別θ。根據(jù)貝葉斯理論，物理參數(shù)θ后驗(yàn)概率p（θ|f,M）可以表示為：

其中，p（f|θ,M）為模型M（θ）的邊緣似然值，即M（θ）在φ的全部先驗(yàn)分布空間內(nèi)似然函數(shù)的平均值。其表達(dá)式為：

其中，p（f|φ,θ,M）為物理參數(shù)θ和超參φ的聯(lián)合似然函數(shù)，p（φ）為φ的先驗(yàn)分布。

獲得θ的后驗(yàn)概率p（θ|f,M）之后，超參后驗(yàn)概率p（φ|f,M）可以表示為：

其中，p（f|φ,M）是模型M的邊緣似然值，即M（φ）在θ的全部后驗(yàn)分布空間內(nèi)似然函數(shù)的平均值。其表達(dá)式為：

邊緣似然值即式（10）（12）表示似然函數(shù)在高維參數(shù)概率分布空間內(nèi)的積分，由于非線性似然函數(shù)的復(fù)雜性，邊緣似然值的解析表達(dá)式難以獲得。本文采用算術(shù)平均估計(jì)方法AME[35?36]計(jì)算式（10）（12），計(jì)算公式為：

式中：χi—模型M參數(shù)概率分布空間內(nèi)的一個(gè)樣本；

m—χi的樣本數(shù)量。

通過AME 計(jì)算邊緣似然值通常需要大量的樣本，m數(shù)值一般為幾十萬至幾千萬，一個(gè)樣本代表運(yùn)行物理模型一次。此外，通過MCMC 方法獲得θ和φ后驗(yàn)分布（即式（9）（11））時(shí)需多次使用AME。因此，模型邊緣似然值的計(jì)算具有嚴(yán)重的計(jì)算負(fù)荷。本文使用當(dāng)前應(yīng)用廣泛的自適應(yīng)稀疏網(wǎng)格（adaptive SG）技術(shù)[37?41]構(gòu)建邊緣似然值的替代模型，克服計(jì)算耗時(shí)困難。

通過式（9）—（13）獲得θ和φ的后驗(yàn)分布之后，進(jìn)行模型預(yù)測。兩步識(shí)別DDM 法計(jì)算步驟為（圖1）：

圖1 兩步識(shí)別DDM 法計(jì)算步驟圖Fig.1 Flow chart of the two-stage based DDM method

（1）從θ的后驗(yàn)樣本中，選擇一組在整個(gè)超參空間內(nèi)表現(xiàn)最優(yōu)的參數(shù)θmax（即邊緣似然值p（f|θ,M）最大），并在物理模型中將其固定為常量；

（2）重新使用MCMC 識(shí)別φ得到其后驗(yàn)分布；

（3）對(duì)于新的模型輸入x?,輸入模型獲得輸出M?，根據(jù)式（2）計(jì)算協(xié)方差矩陣C，C?，C??；

（4）根據(jù)式（5）—（7），使用多元正態(tài)抽樣生成b?，η?；

（5）通過式（8）計(jì)算預(yù)測值f?。

2 地下水溶質(zhì)運(yùn)移案例分析

假設(shè)真實(shí)情況下的溶質(zhì)運(yùn)移過程為非平衡等溫吸附模型，實(shí)際使用線性平衡吸附模型，因此存在模型結(jié)構(gòu)誤差。分別在不考慮模型結(jié)構(gòu)誤差、考慮模型結(jié)構(gòu)誤差（聯(lián)立識(shí)別DDM、兩步識(shí)別DDM）進(jìn)行參數(shù)識(shí)別及模型預(yù)測。

2.1 溶質(zhì)運(yùn)移案例概況

該案例描述的是在定濃度通量邊界條件下，一維穩(wěn)定流溶質(zhì)運(yùn)移模型。真實(shí)情況下的模型為非平衡等溫吸附模型，解析表達(dá)式見式（14）—（17）。其中，輸入濃度C0為50 mol/L，水流速度v為50 cm/d，彌散系數(shù)D為40 cm2/d，注入污染物的時(shí)間t0為5 d，系數(shù)γ為0，衰變常數(shù)μ為0.1，阻滯因數(shù)R為1.3，初始濃度Ci為0 mol/L。

將真實(shí)模型運(yùn)行后得到濃度觀測值，同時(shí)考慮時(shí)間尺度（t:0～10 h）和空間尺度（x:0～250 cm），如圖2所示，分別選取21 個(gè)和36 個(gè)觀測數(shù)據(jù)用于模型識(shí)別和驗(yàn)證。

圖2 模型識(shí)別和驗(yàn)證數(shù)據(jù)點(diǎn)的位置Fig.2 Location of model observation wells during calibration and validation periods

實(shí)際使用的模型為線性平衡吸附模型，解析表達(dá)式為：

基于該模型描述溶質(zhì)運(yùn)移過程時(shí)，污染物吸附模型的偏差會(huì)導(dǎo)致模擬結(jié)果受到模型結(jié)構(gòu)不確定性的影響。

2.2 案例分析

基于線性平衡吸附模型，分別不考慮和考慮模型結(jié)構(gòu)誤差進(jìn)行參數(shù)識(shí)別。假設(shè)輸入濃度C0和水流速度v為隨機(jī)變量，其余參數(shù)與真實(shí)模型相同。表1所示為物理參數(shù)與超參的先驗(yàn)分布。

表1 模型參數(shù)的先驗(yàn)分布Table 1 Prior distributions of model parameters

采用DREAMzs 算法進(jìn)行參數(shù)識(shí)別，設(shè)置3 條平

行的馬爾科夫鏈，每條鏈的預(yù)熱長度為6 000，預(yù)熱后迭代長度為6 000。對(duì)于不考慮和考慮模型結(jié)構(gòu)誤差情況下的參數(shù)識(shí)別，均選取收斂后的18 000 次樣本統(tǒng)計(jì)邊緣后驗(yàn)分布，結(jié)果見圖3。

從圖3 可知，3 種方法得到的物理參數(shù)后驗(yàn)分布差異明顯，考慮模型結(jié)構(gòu)誤差的情況下，參數(shù)后驗(yàn)范圍明顯變小，呈正態(tài)分布的形式。對(duì)于參數(shù)C0和v，聯(lián)立識(shí)別DDM 得到的參數(shù)后驗(yàn)分布范圍均最小，聯(lián)立識(shí)別DDM 和兩步識(shí)別DDM 具有相近的后驗(yàn)均值。兩種DDM 識(shí)別得到的超參邊緣后驗(yàn)分布也存在顯著差異。對(duì)比聯(lián)立識(shí)別DDM，兩步識(shí)別DDM 得到的超參后驗(yàn)分布范圍均較大。對(duì)于λ，兩步識(shí)別法得到的參數(shù)后驗(yàn)分布集中在0.4，聯(lián)立識(shí)別法的后驗(yàn)分布集中在0.3。對(duì)于σ，兩步識(shí)別法的后驗(yàn)分布趨勢平緩，聯(lián)立識(shí)別法后驗(yàn)分布集中于3.0。對(duì)于σδ，兩種DDM得到的后驗(yàn)分布均集中在0.1，聯(lián)立識(shí)別DDM 具有較大的峰值概率密度。因此，忽略模型結(jié)構(gòu)偏差直接進(jìn)行參數(shù)識(shí)別會(huì)導(dǎo)致參數(shù)補(bǔ)償?；贒DM 考慮模型結(jié)構(gòu)偏差時(shí)，物理參數(shù)和超參的獨(dú)立性假設(shè)會(huì)影響參數(shù)識(shí)別結(jié)果。

圖3 識(shí)別得到的模型參數(shù)和超參邊緣后驗(yàn)分布Fig.3 Identified marginal posterior distributions of model parameters and hyperparameters

基于識(shí)別得到的模型參數(shù)后驗(yàn)分布，獲得模型在識(shí)別期和驗(yàn)證期的模擬值，通過頻率統(tǒng)計(jì)獲得模擬值的置信區(qū)間和平均值等信息。本次研究采用3 個(gè)指標(biāo)評(píng)價(jià)模型在識(shí)別期和驗(yàn)證期的表現(xiàn)：均方根誤差（RMSE）、預(yù)測誤差均值（MAE）、相對(duì)誤差的均值（MRE）。這3 個(gè)指標(biāo)值越小，表明模擬值越接近觀測值，模型預(yù)測性能越好，各指標(biāo)統(tǒng)計(jì)結(jié)果見表2。

表2 模型預(yù)測性能指標(biāo)統(tǒng)計(jì)結(jié)果Table 2 Statistics of model prediction performance

從表2 可以看出，在模型識(shí)別期，不考慮模型結(jié)構(gòu)誤差得到的RMSE 值、MAE 值、MRE 值均最低，識(shí)別效果最好；在模型驗(yàn)證期，不考慮模型結(jié)構(gòu)誤差得到的RMSE 值、MAE 值、MRE 值均最高，預(yù)測效果最差。結(jié)果證明，不考慮結(jié)構(gòu)誤差直接進(jìn)行參數(shù)識(shí)別會(huì)產(chǎn)生參數(shù)補(bǔ)償問題，即參數(shù)過度擬合補(bǔ)償結(jié)構(gòu)誤差，從而影響模型預(yù)測效果。對(duì)于兩種DDM，兩步識(shí)別DDM 法在模型識(shí)別期和驗(yàn)證期均得到較低的RMSE、MAE、MRE。因此，通過兩步識(shí)別DDM 法量化模型結(jié)構(gòu)誤差，能夠減少物理參數(shù)和超參獨(dú)立性假設(shè)導(dǎo)致的參數(shù)過度擬合，從而提高模型的預(yù)測效果。

3 地下水水流案例分析

假設(shè)真實(shí)情況下地下水流模型為有越流過程的定流量非完整井流模型，實(shí)際使用的模型沒有考慮越流過程，因此存在模型結(jié)構(gòu)誤差。分別在不考慮和考慮模型結(jié)構(gòu)誤差（聯(lián)立識(shí)別DDM、兩步識(shí)別DDM）條件下進(jìn)行參數(shù)識(shí)別及模型預(yù)測。

3.1 水流案例概況

該案例刻畫了在第一類越流系統(tǒng)中，定流量非完整井流模型。真實(shí)情況下的模型為有越流過程的定流量非完整井流模型（圖4）。抽水井以定流量Q=8 000 m3/d 向外抽水，過濾管位于深度z=d和z=l之間（d=4 m，l=8 m）。觀測井水位表示含水層中（r,z）處的水頭值，降深s的解析表達(dá)式見式（20）。其中，觀測時(shí)刻t=8 d，承壓含水層儲(chǔ)水系數(shù)μ為0.003，垂向滲透系數(shù)Kzz為5 m/d。

圖4 定流量非完整井流模型示意圖Fig.4 Diagrammatic representation of partially penetrating well pumping at a constant rate

將真實(shí)模型運(yùn)行后得到降深觀測值（圖5），考慮不同方向的空間尺度（r:0～20 m；z:0～20 m），時(shí)間固定在t=8 d，分別選取21，58 個(gè)觀測數(shù)據(jù)用于模型識(shí)別和驗(yàn)證。

圖5 模型識(shí)別和驗(yàn)證數(shù)據(jù)點(diǎn)的位置Fig.5 Location of model observation wells during calibration and validation periods

實(shí)際使用的模型沒有考慮越流過程，即弱透水層的厚度M′為0 m，滲透系數(shù)K′不存在，其余參數(shù)設(shè)置和真實(shí)模型相同，解析表達(dá)式見式（20）?；谠撃Ｐ兔枋龆髁糠峭暾樗鸬乃蛔兓瘯r(shí)，由于對(duì)越流過程刻畫的差異，導(dǎo)致模型存在結(jié)構(gòu)誤差，預(yù)測結(jié)果受模型結(jié)構(gòu)不確定性的影響。

3.2 案例分析

基于沒有考慮越流過程的井流模型，分別不考慮和考慮模型結(jié)構(gòu)誤差進(jìn)行參數(shù)識(shí)別。假設(shè)承壓含水層水平滲透系數(shù)Krr和厚度M為隨機(jī)變量，其余參數(shù)與真實(shí)模型相同。表3 為物理參數(shù)與超參的先驗(yàn)分布。

表3 模型參數(shù)的先驗(yàn)分布Table 3 Prior distributions of model parameters

采用DREAMzs 算法進(jìn)行參數(shù)識(shí)別，設(shè)置3 條平行的馬爾科夫鏈，每條鏈預(yù)熱6 000 次，預(yù)熱后迭代6 000 次。不考慮和考慮模型結(jié)構(gòu)誤差（聯(lián)立識(shí)別DDM、兩步識(shí)別DDM）的參數(shù)邊緣后驗(yàn)分布均取收斂后的18 000 次樣本，結(jié)果見圖6。

圖6 識(shí)別得到的模型參數(shù)和超參邊緣后驗(yàn)分布Fig.6 Identified marginal posterior distributions of model parameters and hyperparameters

從圖6 可以看出，3 種方法得到的物理參數(shù)后驗(yàn)范圍存在差異，且后驗(yàn)分布形態(tài)明顯不同。考慮模型結(jié)構(gòu)誤差時(shí)，物理參數(shù)后驗(yàn)范圍均較小，分布更集中。對(duì)于Krr，兩種DDM 的參數(shù)后驗(yàn)分布均集中在9.2 附近，聯(lián)立識(shí)別DDM 具有較大的峰值概率密度。對(duì)于M，兩種DDM 的參數(shù)后驗(yàn)分布均集中在84 附近，兩步識(shí)別DDM 具有較大的概率密度。此外，兩種DDM 得到的超參邊緣后驗(yàn)分布差異明顯，兩步識(shí)別DDM 得到的超參后驗(yàn)分布范圍均較大。對(duì)于λ，兩步識(shí)別法的后驗(yàn)分布集中在0.36，聯(lián)立識(shí)別法的后驗(yàn)分布集中在0.33。對(duì)于σ，兩步識(shí)別法的后驗(yàn)分布趨勢平緩，聯(lián)立識(shí)別法的后驗(yàn)分布集中于5.0。對(duì)于σδ，兩種DDM 得到的后驗(yàn)分布均集中在0.05，聯(lián)立識(shí)別DDM 具有較大的峰值概率密度。因此，基于DDM刻畫模型結(jié)構(gòu)偏差時(shí)，物理參數(shù)和超參的獨(dú)立性假設(shè)會(huì)影響參數(shù)識(shí)別結(jié)果。

基于DREAMzs 算法得到模型在識(shí)別期和驗(yàn)證期的模擬值，通過頻率統(tǒng)計(jì)獲得模擬值的置信區(qū)間和平均值。模型預(yù)測性能指標(biāo)統(tǒng)計(jì)見表4。

表4 模型預(yù)測性能指標(biāo)統(tǒng)計(jì)結(jié)果Table 4 Statistics of model prediction performance

從表4 可以發(fā)現(xiàn)，與地下水溶質(zhì)運(yùn)移案例結(jié)果類似（表2），在模型識(shí)別期，不考慮模型結(jié)構(gòu)誤差得到的RMSE 值、MAE 值、MRE 值均最低，在模型驗(yàn)證期各指標(biāo)值均最高。表明不考慮模型結(jié)構(gòu)誤差直接進(jìn)行參數(shù)識(shí)別會(huì)影響模型預(yù)測效果。對(duì)于兩種DDM，兩步識(shí)別DDM 法在模型識(shí)別期和驗(yàn)證期均得到較低的RMSE、MAE、MRE。因此，兩步識(shí)別DDM 能提高模型的預(yù)測精度。需要說明的是，本次研究提出的兩步識(shí)別DDM 方法適用于任意形式的數(shù)學(xué)模型，包括解析模型和數(shù)值模型，因?yàn)閷?duì)于實(shí)際復(fù)雜的水文地質(zhì)條件，建立的地下水模型均存在結(jié)構(gòu)偏差。

4 結(jié)論

（1）不考慮模型結(jié)構(gòu)誤差直接進(jìn)行參數(shù)識(shí)別會(huì)產(chǎn)生參數(shù)補(bǔ)償問題，即通過參數(shù)過度擬合補(bǔ)償模型結(jié)構(gòu)誤差，從而影響模型的預(yù)測效果。

（2）基于DDM 刻畫模型結(jié)構(gòu)偏差時(shí)，物理參數(shù)和超參的獨(dú)立性假設(shè)會(huì)影響參數(shù)識(shí)別結(jié)果。相對(duì)于聯(lián)立識(shí)別DDM，兩步識(shí)別DDM 得到模型預(yù)測的均方根誤差（RMSE）、均值（MAE）、相對(duì)誤差均值（MRE）均明顯降低。通過兩步識(shí)別DDM 法量化模型結(jié)構(gòu)誤差，能夠減少物理參數(shù)和超參獨(dú)立性假設(shè)而導(dǎo)致的參數(shù)過度擬合，從而能夠更加準(zhǔn)確地刻畫結(jié)構(gòu)誤差，預(yù)測結(jié)果更加可靠。