一、引言

　　隨著近幾十年計算機軟硬件和數值計算方法的發展，地下水建模成為人類深入理解開放、復雜的地下水系統的有效手段。但任何數學模型或數值模擬必然都只是真實地下水系統的簡化，不確定性不可避免地由各種模型輸入因素引入，導致模型的輸出或預測結果也帶有不確定性。為了更好地利用有限的資金和資源去最優化地降低預測結果的不確定性，需要測算各個模型輸入因素的重要性，找出最重要即預測結果對之變化最敏感的不確定模型輸入因素，需要引入敏感性分析。

　　從不確定性角度來理解，敏感性分析是指比較或計算不同的模型不確定性輸入因素對輸出結果不確定性的貢獻的方法。現有的在地下水領域應用的敏感性分析方法大致分為 2 種：局部敏感性分析和全局敏感性分析。全局敏感性方法相較于局部敏感性方法雖然更復雜、計算量更大，但它能夠得到適用于模型參數全部取值范圍的結論，并考慮輸入因素之間的相互作用，因此在最近數十年得到了重視和應用。

　　盡管全局敏感性分析已經被廣泛應用于地下水建模，但多數的敏感性分析研究僅針對地下水模型參數的不確定性，而忽略了另外 2 類地下水建模中所遇到的主要模型輸入不確定性：地下水模型本身的不確定性和情景的不確定性。在地下水建模中，情景不確定性被定義為一種影響地下水系統，并存在于未來的多種可能的自然或者人為(比如降雨、抽水等)外部狀態導致的不確定性。模型的不確定性是指由于缺乏模擬區域或系統的數據或者對其的理解，多種可能的概念或數學模型并存導致的不確定性。

　　DAI 等提出了先情景、后模型，最后參數的不確定性多層級關系，并將方差分析法與多層級的不確定性分析框架相結合，推導出新的全局敏感度系數來評估地下水模型中參數的敏感性。但該研究仍然沒有針對情景和模型的不確定性進行評價。除此之外，前人的多層級不確定性分析框架只能評估系統中模型的敏感性，無法對模型中各種過程的概念或數學模型的敏感性進行分析。在地下水建模中，通常存在多個概念或數學模型用于刻畫同一個地下水流動過程，目前少有研究者對模型中各個過程的概念模型和數學模型進行敏感性分析。

　　為此，筆者提出了一種改進的層級制敏感性分析方法，量化地下水模型中情景、模型和參數的不確定性，并利用參數分組的方法在層級制不確定性分析框架中評估地下水模型涉及的各個過程及其概念或數學模型的不確定性對模型輸出不確定性的貢獻，并將提出的新方法應用于理想的地下水污染運移模擬模型以說明該方法的運用。

　　二、改進的層級制全局敏感性分析方法

　　2.1 基于方差分析的層級制敏感性分析

　　將不同來源的不確定性置于多層級框架的不同層，然后利用方差分析將輸出結果的不確定性按照層級制框架分配到不同的不確定性輸入上，計算各部分的方差在總方差中的占比，并把這種百分比定義為敏感度系數來精確量化不同來源的不確定性的重要性。

　　按照情景 - 模型 - 參數的多層級框架，利用總方差定律對模型輸出結果的總方差進行分解。總方差共被分解為 4 項：情景貢獻的方差、模型貢獻的方差、參數貢獻的方差和觀測數據貢獻的方差。由于觀測數據貢獻的方差獨立于情景、模型和參數貢獻的方差，故本研究不考慮觀測數據的不確定性。使用其余 3 部分的方差與總方差的比值定義情景的敏感度系數、模型的敏感度系數和參數的敏感度系數，3 種敏感度系數在數學上的和為 “1”，該條件也是該方法在應用時的檢驗標準。

　　2.2 過程子模型的敏感度系數

　　一個復雜的地下水系統通常由多個過程耦合組成，例如：地下水流動過程、溶質運移過程、能量傳輸過程等。研究者們通常根據有限的觀測信息利用數學方程簡化水流運動、污染(溶質)運移等真實過程，這就造成每種過程的數學模型并不唯一，即各種過程的概念或數學模型存在不確定性。目前地下水領域少有研究針對此類不確定性。本研究利用參數分組的方法擴展參數的敏感度系數，在多層級不確定性分析框架中評估參數組的敏感性，從而評估模型中各種過程的概念或數學模型的不確定性對輸出結果不確定性的貢獻。

　　將描述某一過程的概念或數學模型包含的所有參數分為一組作為條件對參數貢獻的方差進一步分解，等式右邊第 1 項代表該參數組，即該數學模型或概念模型(稱為過程子模型)的不確定性貢獻的方差，右邊第 2 項代表除該參數組以外其余參數產生的方差。定義過程子模型的敏感度系數為該參數組不確定性貢獻的方差與參數貢獻的方差的比值。

　　三、理想的地下水污染運移模擬模型介紹

　　通過理想的地下水污染運移模擬模型來說明上述方法的運用。寬度 L=10000 m，地下水流和溶質運移考慮為一維穩定流溶質運移情況。均勻穩定入滲的條件下，考慮 3 種情景，以情景 1 作為基準值，其降雨強度 P 為 1524 mm/a，較干燥和濕潤的情景下，降雨強度分別是基準值的 80% 和 180%。兩側河流設定為定水頭邊界，h₁=330 m，河水位 h₂受到融雪過程的影響。在潛水含水層剖面的中間位置 x₀=5000 m 處存在一個污染源，污染物為乙烯(ETH)。

　　對于降雨入滲補給過程，考慮 2 種數學模型來描述(R₁和 R₂)，入滲補給數學模型 R₁中考慮參數取值的不確定性，參數 a 設定為服從均值為 16.88、標準差為 1 的正態分布;入滲補給數學模型 R₂中參數 b 設定為服從均勻分布的隨機參數，分布區間為 [0.1,0.2]。

　　潛水含水層的地質情況共考慮 2 種數學模型來描述，模型 G₁代表整個河間地塊為均質，使用統一的滲透系數 K 來描述;模型 G₂代表河間地塊為非均質，以 x=7000 m 為分界處，左側區域滲透系數為 K₁，右側滲透系數為 K₂，滲透系數均設定為服從正態分布的隨機參數。

　　河流 2 的水頭 h₂受融雪徑流過程的影響，河流 2 水位(h₂)和流量 Q 的關系如下，假設河流徑流量 Q(m³/s)受上游融雪量控制，計算為 Q=C₁×C_sn×S×SVC×A，使用 2 種數學模型來評估日融雪量 S，即度日模型一般形式(SN₁)和融入輻射變量后的度日模型(SN₂)，度日因子 f₁和 f₂分別設定為服從正態分布的隨機參數，數學模型 SN₂中太陽短波輻射或者凈輻射 Rₙ=80 W/m²，系數 r 設定為服從正態分布的隨機參數。

　　綜上，整個潛水含水層一維穩定流溶質運移模型由 3 部分組成：①降雨入滲補給過程可由 2 種數學模型刻畫(過程子模型 R₁和 R₂);②地質過程考慮均質和非均質 2 種過程子模型(G₁和 G₂);③融雪模型可由 2 種數學模型刻畫(過程子模型 SN₁和 SN₂)。通過組合以上 3 種過程的所有子模型得到 8 種模型用于描述該系統，為每個模型分配相等的先驗權重，所考慮的輸出(A)是 5100～5700 m 的所有位置(每 100 m 設置 1 個計算點)的乙烯質量濃度，在 MATLAB 軟件中采用蒙特卡羅模擬進行編程計算。

　　四、結果與討論

　　4.1 5400 m 處乙烯質量濃度預測結果的敏感性分析

　　4.1.1 預測結果的不確定性

　　不同情景下不同模型的預測結果不同，表明預測存在較大的不確定性。其中，包含均質地質過程的模型(均質模型)與包含非均質地質過程的模型(非均質模型)預測結果差異最為明顯，非均質模型預測結果的均值明顯小于均質模型的均值，且非均質模型的不確定性范圍明顯大于均質模型的不確定性范圍。不同情景之間的預測結果差異并不明顯，即在此案例中，情景對結果預測的影響程度并不大。

　　4.1.2 基于方差分析的層級制敏感性分析結果

　　模型的敏感度系數(Sₘ)顯著大于情景和參數的敏感度系數(S_sc、S_pa)，即模型的不確定性是輸出結果不確定性最重要的來源;而情景的敏感度系數最小，意味著情景的不確定性可以忽略不計。僅執行 1000 次運算，敏感度系數就已經全部收斂，相對誤差小于 0.01，模型運行時長僅 2.88 s，相較于前人對該案例的參數敏感性分析，計算效率大幅提升。

　　將總方差分解為情景間方差和情景內方差時，情景內方差幾乎是結果方差的全部來源，情景間方差可以忽略不計。將情景內方差進一步分解為模型間方差和模型內方差時，模型間方差是情景內方差的主要來源，也是總方差的主要來源。模型內方差是 3 種情景下模型內方差的平均結果，非均質模型中所有參數的不確定性是模型內方差的主要來源。

　　基于方差分析的層級制敏感性分析結果顯示，模型間的不確定性是輸出結果的不確定性最重要的來源，這提醒建模者為了顯著降低系統的不確定性，應該把更多的精力和資源放在模型結構上。模型內的不確定性即所有參數貢獻的不確定性明顯小于模型間的不確定性，這也提醒建模者針對該案例，如果花費過多的精力和成本在參數反演上是不明智的選擇。另外，盡管非均質模型的方差對模型內方差的貢獻更大，但這并不代表均質模型可以在建模的時候被忽略。雖然均質模型的方差貢獻小，但值得注意的是，若將該系統的模型確定為均質模型中的一個，系統的總方差將明顯減小，不確定性顯著降低，這意味著均質模型在這樣的系統中應該慎重選擇，合理取舍。

　　4.2 5100～5700 m 范圍內所有位置預測結果的敏感性分析

　　4.2.1 預測結果的不確定性

　　3 種情景下，各個模型預測的均值不盡相同，表明預測存在不確定性。比較不同情景下所有模型的結果可以發現，隨著降雨強度的變化，模型的預測結果并未發生明顯的變化，這意味著在 5100～5700 m 的空間范圍中，情景間的不確定性很小。但所有模型的預測均值曲線在 5300 m 以后分離程度明顯，這表明模型間的預測結果差異很大，即模型間的不確定性很大。

　　在該空間中大多數位置，均質模型的濃度均值明顯高于非均質模型的均值。這是因為非均質模型的等效滲透系數相較于均質模型的滲透系數更小，模擬時間相同時，乙烯在均質模型中運移的更遠，污染范圍更大，故在靠近污染源的位置，非均質模型的預測濃度均值更大;而在遠離污染源的位置，均質模型預測濃度均值更大。

　　4.2.2 基于方差分析的層級制敏感性分析結果

　　在 5100～5700 m 的空間范圍內，情景的敏感度系數(S_sc)始終遠小于其他 2 種系數且幾乎接近于 0，這意味著情景的變化對輸出結果的不確定性并不重要，即可以忽略情景的不確定性。在 5300 m 后的所有位置，對輸出結果的方差貢獻最大的是模型間的不確定性(Sₘ)。在 5200 m 處，模型的敏感度系數(Sₘ)明顯小于參數的敏感度系數(S_pa)，這意味著在該位置系統總方差的主要來源是模型內方差，即所有參數的不確定性貢獻的方差。

　　將總方差分解為情景內方差和情景間方差，所有位置的結果均顯示，95% 以上的總方差來自于情景內方差。將情景內方差進一步分解為模型間方差和模型內方差，可知模型間方差始終是情景內方差的主要來源。各部分方差總體上隨著距離的增大先增大后減小，在 5400 m 處達到峰值。將 3 種情景下的模型內方差進一步分解為 8 個模型的方差，在 5400 m 及其之前的所有位置，非均質模型預測結果的方差是模型內方差的主要貢獻者，而在 5400 m 以后的位置，均質模型預測結果的方差成為模型內方差的主要來源。

　　4.3 過程子模型的敏感性分析

　　為了加深對模型結構的理解，使用本研究提出的過程子模型的敏感度系數精確評估了不同地質模型(G₁，G₂)、入滲補給模型(R₁，R₂)和融雪模型(SN₁，SN₂)對輸出結果不確定性的貢獻。結果顯示，所評估的位置不同，敏感度系數的排名不同。總體而言，地質模型始終是最重要的模型，入滲補給模型的敏感性幾乎接近于 0。這是因為地質模型中的滲透系數直接影響著污染物乙烯的運移速度，而入滲補給模型涉及的參數通過影響地下水位，間接影響乙烯的運移速度。在 5400 m 及其之前的所有位置，地質模型中的非均質子模型 G₂是所有過程子模型中最重要的模型，其敏感度系數顯著大于其余模型的敏感度系數;而在 5500 m 及其以后的位置，均質子模型 G₁是最重要的模型。

　　綜上所述，建議將模型校準的主要精力和資源集中在地質過程的數學模型上。為了提高建模中參數反演效率，可以考慮將入滲補給模型和融雪模型中的不確定性參數固定為某個經驗值，優先反演地質過程的數學模型中包含的參數;或者也可以舍棄不敏感的過程子模型，選取更重要的數學模型來描述過程。

　　五、結語

　　改進的層級制敏感性分析方法能夠量化并比較情景、模型、參數 3 種不同類型的輸入不確定性對輸出結果不確定性的貢獻，同時評估地下水模型中各個過程的概念或數學模型的不確定性對輸出結果不確定性的影響程度，為地下水模型提供精確全面的敏感性分析。相較于以往的參數敏感性分析，改進的層級制敏感性分析方法的計算效率大幅提升，并且能夠提供更有用的敏感性信息。

　　本研究使用一個理想的地下水污染運移模擬模型展示了該方法的運用，結果表明，模型的不確定性是該案例中乙烯濃度預測不確定性的主要來源。對所有過程子模型進行分析后發現，地質模型是模型輸出結果不確定性最重要的來源，其中非均質子模型相較于均質子模型對結果的不確定性貢獻更大。這些信息有助于指導建模者利用有限資源去優先減少該案例的模型不確定性，尤其是描述地質過程的數學模型的不確定性，從而最優地減少預測濃度結果的不確定性。

劉玉姣;戴恒;李躍東;崔節波;文章,中國地質大學(武漢);長江流域環境水科學湖北省重點實驗室;生物地質與環境地質國家重點實驗室;生態環境部華南環境科學研究所,202405