水文地質(zhì)學(xué)中的反演問題

地下水反演是一個綜合性問題，研究重點(diǎn)是對含水層特征進(jìn)行刻畫。許多建模人員已經(jīng)對概念模型的不確定性（顯著的時間和空間變化）、尺度依賴性、許多類型的未知參數(shù)（導(dǎo)水系數(shù)、補(bǔ)給和邊界條件等）、非線性以及對含水層性質(zhì)經(jīng)常表現(xiàn)出低靈敏性的狀態(tài)變量（特別是水頭和濃度）進(jìn)行了研究。由于這些難題的存在，因此不可能把模型校準(zhǔn)與模擬過程分割開來，而是應(yīng)當(dāng)將其作為認(rèn)識含水層特征的重要一步。事實(shí)上，真實(shí)參數(shù)的估計方法盡管在計算細(xì)節(jié)上會有一些差異，但在本質(zhì)上沒有什么區(qū)別。就是否還有充足的空間能夠改進(jìn)地下水反演問題存在爭議：開發(fā)用戶友好的代碼、通過地質(zhì)統(tǒng)計學(xué)方法調(diào)整變異性、將地質(zhì)學(xué)信息和其它類型的數(shù)據(jù)（溫度、同位素濃度和地下水年齡等）相結(jié)合、正確計算不確定性等。盡管如此，即使利用現(xiàn)有的代碼進(jìn)行自動校準(zhǔn)，也仍然極大地推動了模擬工作的發(fā)展。因此，主張將這一方法的應(yīng)用進(jìn)行標(biāo)準(zhǔn)化。一、概述廣義上，反演模擬是指利用模擬參數(shù)的測量值來收集模型信息的一個過程，其中包括兩個相關(guān)的概念：模型識別和參數(shù)估計，在本文中，后者與校準(zhǔn)具有相同的含義。“模型識別”是應(yīng)用一些方法來尋找模型的特征，如控制方程、邊界條件、時間范圍或非均質(zhì)性等。而“參數(shù)估計”是對那些刻畫特征的性質(zhì)進(jìn)行賦值。地下水模擬中的“反演”與上述定義大致相同。事實(shí)上，一些水文地質(zhì)學(xué)反演問題已經(jīng)建立了標(biāo)準(zhǔn)規(guī)程，例如，水文地質(zhì)學(xué)家采用抽水試驗(yàn)來確定含水層的特征。解釋這些試驗(yàn)結(jié)果的過程包括：首先確定最合適的模型（根據(jù)實(shí)際的測量值和地質(zhì)條件），然后根據(jù)曲線擬合來估計模型參數(shù)，利用自動程序可以很容易完成這項工作。在解析解很少的領(lǐng)域中，如非飽和帶，自動化的校準(zhǔn)程序已經(jīng)標(biāo)準(zhǔn)化了（Kool和Parker，1988；Hollenbeck和Jensen，1998）。在通過適當(dāng)?shù)臏y試刻畫空間變化的領(lǐng)域中，如醫(yī)學(xué)（掃描或X線斷層攝影術(shù)）（Rudin等，1999；Rhli等，2002）和地球物理學(xué)（Sambridge和Mosegaard，2000）等等，自校準(zhǔn)已經(jīng)成為標(biāo)準(zhǔn)。然而，地下水模擬顯示出了能夠區(qū)別于其它模擬的幾個特性：

成本地下水模擬費(fèi)用相對較高。模型每次運(yùn)行時，都要求解大型方程組（以獲取系統(tǒng)準(zhǔn)確而真實(shí)的特征）。另外，計算機(jī)科學(xué)的發(fā)展可以大大提高模型的性能；

時間依賴性狀態(tài)變量，如水位和濃度具有時間依賴性。許多流動問題本質(zhì)上都處于穩(wěn)定流條件下，但即使如此，仍然會包含水位隨時間波動的信息，這樣就需要進(jìn)行非穩(wěn)定流的模擬；

非均質(zhì)性水力傳導(dǎo)系數(shù)K通常是最重要的水力學(xué)指標(biāo)，其值會相差幾個數(shù)量級；導(dǎo)水系數(shù)也是如此，相當(dāng)于二維模型中的K（在下文中，這兩個參數(shù)將會交替使用）；

不同類型的參數(shù)。當(dāng)重點(diǎn)考慮導(dǎo)水系數(shù)時，其它參數(shù)（補(bǔ)給、邊界流量等）同樣具有不確定性和相關(guān)性；

尺度依賴性野外測量的參數(shù)通常只能表示含水層的一小部分，這樣，在性質(zhì)和數(shù)量上就與模型所需的參數(shù)具有很大的差異；

模型不確定性含水層的幾何特征和非均質(zhì)性受地質(zhì)條件的控制，而地質(zhì)條件很難得到確切的認(rèn)識。

靈敏度低在不同的問題中，狀態(tài)變量可能會對模擬參數(shù)表現(xiàn)出很低的靈敏度（即信息量不足）。特別是水頭（最常用的，有時是唯一的測量指標(biāo)）有時幾乎不包含關(guān)于水力傳導(dǎo)系數(shù)的信息。由于以上特征，含水層模型預(yù)測具有極高的不確定性。另外，參數(shù)估計不像其它領(lǐng)域那樣可以明確地表達(dá)出來，在解釋抽水試驗(yàn)結(jié)果時，就可能遇到這種情況。研究人員會將試驗(yàn)獲得的水頭值直接輸入到程序中，經(jīng)過適當(dāng)?shù)亩ㄐ阅Ｐ头治龊?，得到參?shù)值。反演是含水層模擬過程中的一步重要工作。地下水的姐妹學(xué)科——地表水文學(xué)也具有許多類似的特征，其中的很多問題已經(jīng)得到解決。地下水反演的另一個重要特點(diǎn)是其相對獨(dú)立性。目前，已經(jīng)開發(fā)出了很多獨(dú)立于其它領(lǐng)域的反演方法。最早的方法是將簡單的等效水頭（假定已知）代入流動方程中，這樣可以得到關(guān)于導(dǎo)水系數(shù)的一階偏微分方程（Stallman，1956）。Neuman（1973）將這一方法稱為“直接法”，這種方法相對簡單、易于理解，在Nelson（1960，1961）之后，得到了廣泛應(yīng)用。事實(shí)上，根據(jù)測量水頭畫出流網(wǎng)后，可以得到導(dǎo)水系數(shù)（Bennet和Meyer，1952）。然而，這一方法存在幾個不足之處：首先，整個研究范圍內(nèi)（時間和空間）的水頭值（以及補(bǔ)給、貯水系數(shù)和邊界條件）必須是已知的，這樣才能進(jìn)行插值工作，但同時也會引入誤差，使最終估計的導(dǎo)水系數(shù)受到了一些人為影響；其次，這一方法不穩(wěn)定，水頭微小的誤差會造成導(dǎo)水系數(shù)很大的誤差。為了解決第一個問題，最近大部分的反演方法采用了Neuman（1973）所稱為的“間接法”，這一方法承認(rèn)測量存在誤差，并尋找可以將這些誤差最小化的水力性質(zhì)。也就是說，通過對目標(biāo)函數(shù)的最小化來尋找參數(shù)，這一工作的計算量會很大。為了克服不穩(wěn)定性，采取了一系列方法：在目標(biāo)函數(shù)中增加一個正則項（調(diào)整項）以消除無根據(jù)的擺動、考慮其它類型的數(shù)據(jù)、減少被估參數(shù)的數(shù)量（但不丟失表現(xiàn)空間變化的能力）等。這些問題（計算、穩(wěn)定、不同類型的數(shù)據(jù)和空間變化等）一直是許多研究的焦點(diǎn)，如Yeh（1986）、Carrera（1987）、Kool等（1987）、McLaughlin和Townley（1996）以及deMarsily（1999）等人的研究。本文的目的是介紹含水層特征的反演模擬研究現(xiàn)狀，重點(diǎn)是尋找現(xiàn)有方法的相似點(diǎn)，探討如何去應(yīng)用，而非是對這些方法的簡單論述。校準(zhǔn)是否應(yīng)當(dāng)作為模擬過程的必要步驟尚存在爭議，一般的模擬過程包括：第一步，收集真實(shí)系統(tǒng)所有可用的信息，建立合理的概念模型；第二步，采用控制方程來表示概念模型；第三步，離散并求解控制方程，同時將一些未知的水力參數(shù)（即面狀補(bǔ)給）離散化作為模型參數(shù)的函數(shù)；第四步，優(yōu)化模型參數(shù)；第五步，如果需要的話，應(yīng)當(dāng)進(jìn)行誤差分析，根據(jù)新的試驗(yàn)結(jié)果將一些最初的概念模型舍棄，保留有效的模型。本文將在以下章節(jié)中進(jìn)一步分析知識和數(shù)據(jù)在確定概念模型時的重要性，通常認(rèn)為這是地下水模擬過程中最常見的誤差源（Bredehoeft，2004）。在“需要估計什么”一節(jié)中介紹如何利用模型參數(shù)來表示含水層性質(zhì)，下一步工作是優(yōu)化（校準(zhǔn)）模型參數(shù)。在“如何估計模型參數(shù)”和“代碼如何工作”的章節(jié)中介紹了校準(zhǔn)方法和程序。獲得參數(shù)之后，就需要對模型的質(zhì)量進(jìn)行檢驗(yàn)，這是“模型性能如何”一節(jié)的主要內(nèi)容。在“實(shí)際做什么”一節(jié)中，將討論實(shí)際的應(yīng)用趨勢。文末討論了將來所面臨的挑戰(zhàn)。二、概念模型：知識和數(shù)據(jù)知識和數(shù)據(jù)是建立概念模型的基礎(chǔ)。已經(jīng)通過許多方法對“知識和數(shù)據(jù)”進(jìn)行了定義。通常認(rèn)為“知識”是對真實(shí)事物的認(rèn)識?？茖W(xué)和哲學(xué)所關(guān)注的是這一認(rèn)識是否客觀而真實(shí)。這樣，只有客觀而真實(shí)的認(rèn)識才被認(rèn)為是“知識”。此處的“數(shù)據(jù)”是指關(guān)于含水層的所有信息，而不僅僅是野外測量的數(shù)據(jù)。地下水建模的基礎(chǔ)是對地下水行為的科學(xué)認(rèn)識和詳細(xì)的場地資料，用數(shù)學(xué)方程來表達(dá)對研究對象的認(rèn)識，通常會將這些方程處理成一系列可進(jìn)行數(shù)值求解的離散方程。數(shù)據(jù)既可以定性又可以定量地定義概念模型。表1是一些數(shù)據(jù)類型及其在模擬過程中的作用?？梢愿鶕?jù)不同的標(biāo)準(zhǔn)對數(shù)據(jù)進(jìn)行分類。其中一種分類方式是將數(shù)據(jù)劃分為定性（“軟”）和定量（“硬”）數(shù)據(jù)；另一種分類方式則考慮了數(shù)據(jù)的時間行為：靜態(tài)數(shù)據(jù)指的是不隨時間變化的量，而動態(tài)數(shù)據(jù)指的是隨時間變化的量；第三種分類方式是將數(shù)據(jù)劃分為被估參數(shù)和狀態(tài)變量（水頭和流量等）的測量值。表1水文地質(zhì)學(xué)常用的數(shù)據(jù)分類類型實(shí)例應(yīng)用知識達(dá)西定律和維度等控制方程軟數(shù)據(jù)定性和定量（間接）地質(zhì)（非均質(zhì)性）、地球物理學(xué)、遙感、粒度分析、測量誤差分布用參數(shù)表示、校準(zhǔn)硬數(shù)據(jù)定量抽水試驗(yàn)、示蹤試驗(yàn)、水頭、化學(xué)組分的濃度、其它測量驗(yàn)前信息、校準(zhǔn)、確認(rèn)通過測量獲得的定量數(shù)據(jù)會存在誤差。因此，應(yīng)當(dāng)采用概率密度函數(shù)（下文稱之為pdf）而非單一值來定義測量。在研究過程中，經(jīng)常會采用高斯分布，因?yàn)橐炎C實(shí)它適合用于描述很多變量（如導(dǎo)水系數(shù)的對數(shù)，下文簡稱log-T），另外，高斯分布可以采用平均值和標(biāo)準(zhǔn)差來描述。假設(shè)這一分布易于利用，經(jīng)常要對不服從高斯pdf的數(shù)據(jù)進(jìn)行轉(zhuǎn)換，得到的pdf接近于高斯分布（如導(dǎo)水系數(shù)）。這是一種常見的誤差源，因?yàn)閘og-T的點(diǎn)值服從高斯分布并不表示空間分布也服從多高斯分布（Gómez-Hernández和Wen，1998）。事實(shí)上，關(guān)于這一問題有很多爭議（Meier等，1999），但Log-T的多高斯特征是水文地質(zhì)學(xué)中最常見的假設(shè)。利用地質(zhì)統(tǒng)計學(xué)軟件來產(chǎn)生另一種可選擇的空間分布，例如，采用馬爾可夫鏈（Markov-chain）法（Weissmann等，1999）。狀態(tài)變量原始測量值的轉(zhuǎn)換也有助于增加數(shù)據(jù)對參數(shù)的靈敏度。在地表水水文學(xué)中，這一問題已得到了很多關(guān)注（Meixner等，1999）。在地下水研究工作中，相關(guān)的一個實(shí)例是模擬穿透曲線時，采用排泄點(diǎn)處的觀測總量而非單個濃度值。也可以利用峰值濃度和峰值時間（如Woodbury和Rubin，2000）。另外，采用不同類型的數(shù)據(jù)或利用權(quán)重方法也可增加數(shù)據(jù)的靈敏度（如Wagne和Gorelick，1987；Anderman和Hill，1999）。在利用不同類型的數(shù)據(jù)時，沒有通用的規(guī)則來確定數(shù)據(jù)的信息量（Larocque等，1999）。溫度資料可以反映關(guān)于垂直流量的信息（Woodbury等，1987）。根據(jù)環(huán)境同位素和測年資料，可以了解區(qū)域流動趨勢（Varni和Carrera，1998）；Chen等（2003）采用沉降率來幫助確定隔水層的滲透性。Hill（1992）采用河川徑流的收益和損失來進(jìn)行模型校準(zhǔn)。總之，有大量的狀態(tài)變量資料可以利用?？梢圆捎貌煌绞綄?shù)據(jù)應(yīng)用到反演過程中，通常考慮用包含狀態(tài)變量的觀測值h和模型參數(shù)p的矢量ω來表示大部分?jǐn)?shù)據(jù)：

（1）

式中，h和p表示實(shí)際值，h*和p*表示測量值（或適當(dāng)?shù)霓D(zhuǎn)換值，p*通常稱為“驗(yàn)前估計”），和表示誤差。狀態(tài)變量測量值的先驗(yàn)誤差獨(dú)立于那些參數(shù)的誤差，這樣的假設(shè)比較合理。事實(shí)上，所有測量誤差都是獨(dú)立的，假設(shè)通常都比較合理。上文提及的轉(zhuǎn)換會造成數(shù)據(jù)之間的依賴性增加。對于模型參數(shù)更是如此，測量值與模型所需的數(shù)據(jù)絕不會完全吻合。為了獲得p*，需要利用所有可用的參數(shù)資料，如參數(shù)的點(diǎn)測量值和軟數(shù)據(jù)（如地球物理圖像）等。用一個數(shù)據(jù)矢量zp包括所有這些測量值，之后利用zp，將p的驗(yàn)前估計和協(xié)方差矩陣用線性條件期望方程的一般形式來表示。\o"點(diǎn)擊圖片看全圖"（2）\o"點(diǎn)擊圖片看全圖"（3）式中，Qpp、Qzze和Qpze分別表示模型參數(shù)的協(xié)方差、zp中的測量值、模型參數(shù)和zp中測量值的交叉協(xié)方差；E[·]表示期望值。實(shí)際上，評價這些矩陣的方法與p和zp的性質(zhì)有關(guān)，如果zp表示點(diǎn)測量值，則可以通過協(xié)方差函數(shù)獲得Qzze（即Qzze是簡單克里金矩陣）。在這種情況下，根據(jù)協(xié)方差函數(shù)也可以求得Qpze，如果p表示log-T的塊值，則Qpze可能是平均值。如果p的先驗(yàn)期望E(p)、zp和未知，則需要擴(kuò)展這些矩陣，來加強(qiáng)無偏條件。另外，如果zp不僅包括點(diǎn)測量結(jié)果，還包括軟數(shù)據(jù)，則在Qzze中應(yīng)當(dāng)包括這些數(shù)據(jù)之間的相關(guān)性，Qzze會成為一個協(xié)克里金矩陣。由于交叉協(xié)方差很難估計，因此，通常采用普通克里金法、外部漂移克里金法或其它非平穩(wěn)的克里金法。顯然，計算這些矩陣非常復(fù)雜，詳見一些地質(zhì)統(tǒng)計學(xué)的文章，如Kitanidis（1997）和Rubin（2003）。Carrera等（1993）對如何處理大量不同類型的數(shù)據(jù)進(jìn)行了逐步介紹。然而，無論如何構(gòu)建這些矩陣，都必須明確以下兩個問題：除了嚴(yán)格定性的數(shù)據(jù)（如地質(zhì)描述），還需要提供關(guān)于狀態(tài)變量和模型參數(shù)的所有類型的數(shù)據(jù)；對這些數(shù)據(jù)進(jìn)行統(tǒng)計描述，以確定它們在反演過程中的權(quán)重。三、需要估計什么控制方程的系數(shù)是水力性質(zhì)，需要對這些值進(jìn)行估計。除了具有梯度或張量性質(zhì)，所有這些水力性質(zhì)都具有空間變異性，其中一些也隨時間而變化。獲得連續(xù)模型范圍內(nèi)每一點(diǎn)的水力性質(zhì)是不可能的。因此，需要采用離散表示法。根據(jù)少量的模型參數(shù)來表示水力性質(zhì)的方法稱為參數(shù)化。選擇模型參數(shù)并不是件容易的事，建模者往往采用很多自由度來對變異性進(jìn)行準(zhǔn)確描述。然而，如果模型參數(shù)量很大，那么反演可能會成為不適定問題（Hadamard，1902）。這樣就需要折衷，這是參數(shù)化方法的動機(jī)。許多研究人員將含水層性質(zhì)表示為未知模型參數(shù)的線性組合，這樣，任何水力性質(zhì)q（如面狀補(bǔ)給）都可以表示為：\o"點(diǎn)擊圖片看全圖"

（4）

式中，N是模型參數(shù)pi的數(shù)量，q0(x)是水力性質(zhì)的附加因子（初始模型參數(shù)為0時的q的初始值），是將面狀補(bǔ)給在時間和空間上進(jìn)行參數(shù)化的內(nèi)插函數(shù)。通過離散化可以將式（4）改寫為矩陣形式：q=q0+Aqp（5）式中，Aq是面狀補(bǔ)給內(nèi)插函數(shù)矩陣，p是模型參數(shù)矢量。這樣定義參數(shù)化方法的根據(jù)是：（1）對模型參數(shù)p的估計；（2）計算矩陣A的內(nèi)插方法；（3）水力性質(zhì)q0的初始分布。參數(shù)化方法不同于定義這些矢量和矩陣的方法，最常用的方法如下：（一）分區(qū)法將研究區(qū)分為若干個子區(qū)進(jìn)行參數(shù)化。一般而言，模型參數(shù)的每個矢量組成pi與一個子區(qū)相關(guān)。這樣在每個子區(qū)內(nèi)，假定屬性q(x)為常數(shù)或者規(guī)定以事先確定的方式變化，如果點(diǎn)x落在所研究的子區(qū)之外，則認(rèn)為式（4）中內(nèi)插函數(shù)的值為0。通常時間和空間變化是分開的（Carrera和Neuman，1986），因此：q=(x,t)=pifx(x)ft(t)（6）式中，fx是空間函數(shù)（例如，如果pi表示子區(qū)的水力傳導(dǎo)系統(tǒng)，q表示導(dǎo)水系數(shù)，則fx可以表示含水層厚度）；ft是時間函數(shù)（例如，如果q是面狀補(bǔ)給，ft可以表示根據(jù)土壤質(zhì)量平衡獲得的隨時間變化的補(bǔ)給量，pi是未知系數(shù)）。分區(qū)的主要優(yōu)點(diǎn)是在分析地質(zhì)信息時的一般性和靈活性，例如，分區(qū)可以表示全部或部分地質(zhì)單元等。然而，需要強(qiáng)調(diào)的是，分區(qū)并不排斥地質(zhì)統(tǒng)計學(xué)方法的使用。事實(shí)上，Clifton和Neuman（1982）綜合應(yīng)用了分區(qū)法和克里金法。當(dāng)子區(qū)不必很大時，在保證地質(zhì)單元連貫性的同時，分區(qū)的目的就是降低p的維數(shù)（Stallman，1956）。事實(shí)上，Carrera等（1993）認(rèn)為，關(guān)于參數(shù)變化的可利用的地質(zhì)信息非常具有說服力（可以歸類為確定性數(shù)據(jù)）時，比常規(guī)地質(zhì)統(tǒng)計學(xué)方法要更勝一籌。有些學(xué)者認(rèn)為分區(qū)法過于嚴(yán)格，因此，提出了對子區(qū)形狀進(jìn)行優(yōu)化的研究。特別引人關(guān)注的是在校準(zhǔn)過程中得到的子區(qū)形狀的“幾何變形”（Roggero和Hu，1998）。（二）點(diǎn)估計子區(qū)大小趨于0（事實(shí)上是趨于元素或柵格大小）的分區(qū)方式可作為分區(qū)法的特例，此時仍可以采用式（4）（如Meier等，2001）。然而，參數(shù)的空間維數(shù)就會很大，很難找到更合適的替換公式（Kitanidis和Vomvoris，1983；Dagan，1985；McLaughlin和Townley，1996）。下文會作進(jìn)一步說明。（三）啟發(fā)式內(nèi)插函數(shù)可以任意選擇式（4）中的內(nèi)插函數(shù)αi。選擇不同的類型，包括有限元（Yeh和Yoon，1981）和嶺函數(shù)（Mantoglou，2003）等。這些方法非常靈活，但是在如何定義模型參數(shù)的先驗(yàn)信息方面不是很清楚。（四）試驗(yàn)點(diǎn)法在這種情況下，p表示在一組試驗(yàn)點(diǎn)中，屬性q的未知值。這一方法最早由deMarsily等（1984）提出，現(xiàn)在已經(jīng)得到廣泛應(yīng)用（Ramarao等，1995；Vesselinov等，2001；Hernandez等2003），已經(jīng)成為非線性地質(zhì)統(tǒng)計反演的標(biāo)準(zhǔn)方法。在試驗(yàn)點(diǎn)測量值不確定的情況下，可以考慮將式（4）進(jìn)行一般化：\o"點(diǎn)擊圖片看全圖"（7）

式中，和分別表示測量值和試驗(yàn)點(diǎn)的克里金加權(quán)系數(shù)。對式（7）和式（4）進(jìn)行比較，顯然，q0是等式右邊的第一項，\o"點(diǎn)擊圖片看全圖"與相等。這一方法的主要優(yōu)點(diǎn)是相對簡單而且靈活，然而，直到2004年才認(rèn)識到需要考慮參數(shù)的先驗(yàn)信息。傳統(tǒng)上，式（7）的第一項是通過直接測量獲得的，但是在反演過程中不加以考慮。Doherty（2003）針對模型參數(shù)的不均勻性制定了調(diào)整標(biāo)準(zhǔn)，但沒有用到先驗(yàn)信息。

（五）條件模擬到目前為止，所提出的方法都是基于尋求某種最佳估計。在以下章節(jié)可以看出，有時更傾向于尋求q(x,t)的可能模擬結(jié)果。Neuman和Wierenga（2003）提出了一種全面的模擬方法。Sahuquillo等（1992）、Gómez-Hernández等（1997）和Capilla等（1998）提出一種自校準(zhǔn)方法，其中用q0(x)表示根據(jù)所有可用信息（包括軟數(shù)據(jù)和硬數(shù)據(jù)）得到的log-T模擬結(jié)果，αi(x,t)pi表示由一組主點(diǎn)masterpoint造成的擾動。另一種方法是將q(x)作為隨機(jī)函數(shù)的線性組合（Roggero和Hu，1998；Hu，2002），其中表示條件模擬，模型參數(shù)只是那些模擬的加權(quán)。總之，首先，需要表示水力性質(zhì)的變化；其次，最常用的參數(shù)化方法的基本表達(dá)式如式（4）?，F(xiàn)在的關(guān)鍵問題是如何找到模型參數(shù)p。四、如何估計模型參數(shù)p估計問題涉及到“最佳”模型參數(shù)的概念。問題是“最佳”意味著什么？關(guān)于這一問題沒有理想的答案。事實(shí)上，可能沒有哪一組模型參數(shù)可以代表“最佳”的現(xiàn)實(shí)情況，這也是條件模擬方法的動力——尋找可能的參數(shù)集合。下文將介紹一些最廣泛使用的方法。（一）優(yōu)化方法在這些方法中，將參數(shù)集合定義為目標(biāo)函數(shù)的最小值。確定目標(biāo)函數(shù)的最初目標(biāo)是保證計算值和測量值之間可以很好地擬合。擬合度可以由計算水頭和測量水頭之間的方差來表示，在數(shù)學(xué)術(shù)語中，這個方差是一個范數(shù)（稱為L2范數(shù)）。然而，也可以利用其它一些范數(shù)，而且有一些已經(jīng)應(yīng)用到地下水研究工作中。例如，可以用測量值和計算值之差的絕對值（L1范數(shù)）來量化誤差。每個范數(shù)都有其優(yōu)點(diǎn)：L2范數(shù)比L1范數(shù)更敏感，L1范數(shù)比L2范數(shù)更難應(yīng)用（Woodbury等，1987；Xiang等，1992）。定義性能指標(biāo)其實(shí)也就是一個最小化問題，即尋求使所選范數(shù)最小的模型參數(shù)。采用L2范數(shù)，將最小化目標(biāo)函數(shù)寫為：\o"點(diǎn)擊圖片看全圖"（8）式中，是權(quán)重矩陣。當(dāng)Fh作為目標(biāo)函數(shù)時，問題經(jīng)常是不適定的。目標(biāo)函數(shù)的解也變得不穩(wěn)定，即參數(shù)集不同，但得到的Fh值卻非常相似。有時對初始參數(shù)表現(xiàn)出極度的敏感性，這樣模擬人員會認(rèn)為解是不唯一的。為了解決這些問題，一些研究人員選擇設(shè)置估計參數(shù)的上限和下限。事實(shí)上，這樣并沒有使問題得到解決。解決的只是在上下限之間波動，而并沒有提高可靠性。于是，Neuman（1973）在目標(biāo)函數(shù)中引入一項Fp：F=Fh+λFp（9）（10）該等式最初是根據(jù)擬合和穩(wěn)定性提出的。根據(jù)Tihonov（1993）的觀點(diǎn)，可以將Fp作為正則項（調(diào)整項）。在水文學(xué)中，Emselem和deMarsily（1971）采用該項來降低擺動。當(dāng)Ch作為測量誤差的協(xié)方差矩陣時，可以根據(jù)統(tǒng)計平均值得到式（9）。Gavalas等（1976）將模型參數(shù)的后驗(yàn)概率密度函數(shù)最大化（最大后驗(yàn)法，MAP）來推導(dǎo)式（9）。得到的目標(biāo)函數(shù)是λ＝1時的式（9）。Carrera和Neuman（1986）通過將參數(shù)可能性最大化（最大似然估計）推導(dǎo)式（9）。把式（9）放在統(tǒng)計學(xué)背景下進(jìn)行敘述的優(yōu)點(diǎn)是，得到的方法不僅能估計控制含水層性質(zhì)的參數(shù)，也能估計那些控制不確定性的參數(shù)（方差、變異函數(shù)和似然性）。事實(shí)上，后者和前者同樣重要（Zimmerman等，1998）。將目標(biāo)函數(shù)最小化是一項艱巨的工作，這是由于狀態(tài)變量和參數(shù)之間的關(guān)系通常是非線性的。這也是反演問題的公式劃分為線性和非線性的原因。（二）線性方法不考慮用公式表示，可以將問題線性化為：h(p)=h0+Jhp(p-p*)（11）式中，Jhp是雅可比矩陣行列式（靈敏度矩陣），包括關(guān)于p的h的導(dǎo)數(shù)。后面的章節(jié)會對其計算方式進(jìn)行討論（靈敏度方程）。假設(shè)這個方程有效，可以將式（2）進(jìn)行一般化，得到給定h條件下的p的條件期望（Dagan，1985；Rubin和Dagan，1987）。這樣，就需要確定數(shù)據(jù)和參數(shù)Qph的交叉協(xié)方差以及Qhh的協(xié)方差矩陣，假定：

(12)

HYPERLINK"/uploadfile/200731819154341.gif"

(13)則：

(14)通過采用與式（2）相似的方式，擴(kuò)展矩陣和，就可以將式（14）看作協(xié)克里金（Kitanidis和Vomvoris，1983；Hoeksema和Kitanidis，1984），是通過將估計參數(shù)的方差最小化得到的。最重要的是式（14）和克里金估計值與幾何參數(shù)無關(guān)，一旦確定Qhh，就可以從理論上估計任一點(diǎn)的p。事實(shí)上，Kitanidis和Vomvoris（1983）認(rèn)為需要估計的唯一參數(shù)就是表征log-T場的統(tǒng)計性質(zhì)的參數(shù)。另外，也需要知道測量誤差的協(xié)方差。（三）非線性方法根據(jù)式（11）引出一個問題：在參數(shù)先驗(yàn)估計時為什么要進(jìn)行線性化？事實(shí)上，Carrera和Glorioso（1991）證明圍繞估計參數(shù)本身線性化更好一些，這樣更接近于真實(shí)值。在實(shí)際工作中，可以通過反復(fù)迭代進(jìn)行線性化。Carrera和Glorioso（1991）研究表明，結(jié)果與將式（9）最小化的結(jié)果相同，式（14）的結(jié)果與將式（9）按照高斯－牛頓方法進(jìn)行一級迭代進(jìn)行最小化的結(jié)果相同。這樣就形成一個閉合環(huán)。對提及的估計方法的總結(jié)見表2：表2估算方法總結(jié)方法類型評估對象運(yùn)算法則（等式編號）研究人員最小二乘法非線性Fh（8）Cooley（1977）、Hill等（1998）最大似然估計非線性-2lnP(p/h)（9）CarreraandNeuman（1986）條件期望線性或非線性E(p/h)（14）Dagan（1985）、Carrera和Glorioso（1991）協(xié)克里金法線性或非線性方差估計（14）Kitanidis和Vomvoris（1983）、Carrera等（1993）最大后驗(yàn)法非線性P(p/h)（9）其中λ＝1Gavalas等（1976）各種方法的差異與統(tǒng)計參數(shù)的刻畫和對不確定性的評價有關(guān)。本文沒有討論這些方法如何處理上面兩個問題，而是著重分析它們的重要性。五、代碼如何工作？許多代碼基本上都是按照如下的迭代過程運(yùn)行的，如：（1）初始化：讀取輸入數(shù)據(jù)，設(shè)置迭代計數(shù)器i=0，初始化參數(shù)p0；（2）求解模擬問題h(pi)，計算目標(biāo)函數(shù)Fi及其梯度gi（假定是連續(xù)可微分的）和雅可比矩陣行列式Jhp；（3）計算校正矢量d，可能會利用到先前的迭代信息以及gi和Jhp；（4）更新參數(shù)，pi+1=pi+d；（5）如果已達(dá)到收斂，則停止計算，否則，取i=i+1返回到第（2）步。除了確定校正矢量d以及計算g和Jhp以外，其余步驟都非常簡單。最常用方法如下：（一）計算校正矢量dd的計算屬于優(yōu)化范疇。關(guān)于數(shù)值優(yōu)化方法的文獻(xiàn)很多，但是只有少量應(yīng)用于求解地下水反演問題（表3）。Cooley（1997）提出采用Marquardt方法（Marquardt，1963），這一方法的基礎(chǔ)是把目標(biāo)函數(shù)的二次逼近最小化，同時限制每次迭代的步長大小。這一方法效果非常好，通常經(jīng)過少量迭代后就會收斂，但是由于需要計算雅可比矩陣，因此比較耗時。Carrera和Neuman（1986）提出將擬牛頓方法和共軛梯度法結(jié)合使用。與Marquardt方法相比，這一方法的效果要差一些，但每次迭代的計算量要小得多，僅需要計算目標(biāo)函數(shù)的梯度。方法的選擇取決于研究問題。Cooley（1985）對這兩種方法進(jìn)行比較，認(rèn)為Marquardt方法要好一些。表3地下水反演模擬方法比較方法計算指標(biāo)收斂級收斂1研究人員高斯－牛頓，MarquardtF2、g3和J42局部Cooley（1977）共軛梯度F和g1局部Carrear和Neuman（1986）模擬退火遺傳算法F0全局Rao等（2003），Tsai等（2003）備注：1、全局收斂表示這種方法可以從理論上擺脫局部極?。痪植渴諗勘硎具@種方法很難擺脫一個局部極小值；2、目標(biāo)函數(shù)；3、梯度矢量；4、雅可比矩陣。大部分代碼都利用了這些方法，一般稱之為“下降法”，因?yàn)檫@些方法試圖在每次迭代中改進(jìn)目標(biāo)函數(shù)。因此，它們往往易收斂于局部極小?？朔植繕O小值是許多方法的研究動機(jī)，如模擬退火或遺傳算法（Tsai等，2003）。已證明在地表水水文學(xué)研究中一種非常有效的方法——SCE算法（ShuffledComplexEvolution）（Duan等，1992），但尚未應(yīng)用在地下水領(lǐng)域中。這些方法在實(shí)際工作中尚未得到廣泛應(yīng)用，原因是：當(dāng)參數(shù)數(shù)量很大時（大于20），方法的性能會降低，這在含水層模擬工作中經(jīng)常出現(xiàn)；另外，由于需要對目標(biāo)函數(shù)進(jìn)行大量評價，因此導(dǎo)致地下水模型的成本過高。比較不同優(yōu)化方法時，前述問題的復(fù)雜性具有一定的關(guān)聯(lián)。例如，如果直接問題是非線性的，則與計算梯度或雅各比矩陣相比，評價目標(biāo)函數(shù)的成本會更高，因?yàn)檫@些計算是無法迭代的。這樣在解決非線性問題時，傾向于采用二階收斂而非一階或零階收斂的方法。（二）靈敏度方程靈敏度是狀態(tài)變量對模型參數(shù)的導(dǎo)數(shù)。這一方程非常有用，主要有兩個原因：一是可以應(yīng)用于上文提及的一些優(yōu)化方法中，二是可以得到關(guān)于模型和估計參數(shù)可靠性的信息。另外，由于根據(jù)該方程可以確定哪些數(shù)據(jù)對哪些參數(shù)最靈敏，因此可以用于設(shè)計優(yōu)化網(wǎng)絡(luò)（見“模型性能如何”一節(jié)）。靈敏度一般可以通過以下三種方法獲得：直接求導(dǎo)法、伴隨狀態(tài)法和有限差分法。Carrera等（1990）、Carrera和Medina（1994）對這些方法進(jìn)行了比較。直接求導(dǎo)法是對模擬方程簡單求導(dǎo)，在穩(wěn)定條件下，可以得到：\o"點(diǎn)擊圖片看全圖"（15）需要注意，在式（15）中，可以計算節(jié)點(diǎn)水頭對參數(shù)的導(dǎo)數(shù)，但是采用有限元法（FEM）的內(nèi)插函數(shù)可以計算空間上每一點(diǎn)的導(dǎo)數(shù)。同時也可以計算與水頭有關(guān)的任何變量的導(dǎo)數(shù)。伴隨狀態(tài)法是根據(jù)將F的最小化作為關(guān)于h和p的優(yōu)化問題，同時將狀態(tài)方程作為等式約束。伴隨狀態(tài)矢量λ是聯(lián)合優(yōu)化問題的拉格朗日乘子的集合，可以根據(jù)下式獲得：\o"點(diǎn)擊圖片看全圖"（16）這樣目標(biāo)函數(shù)的梯度就可以寫為：HYPERLINK"/uploadfile/2007318192059990.gif"

（17）伴隨狀態(tài)方程也可以用于計算雅可比矩陣，當(dāng)觀測點(diǎn)數(shù)量小于參數(shù)數(shù)量時，這一方法具有一定優(yōu)勢，也可以準(zhǔn)確計算Hessian矩陣（目標(biāo)函數(shù)對模型參數(shù)的二階導(dǎo)數(shù)）（詳見Carrera和Medina，1994；Medina和Carrera，2003的研究成果）。有限差分是根據(jù)增量比取近似導(dǎo)數(shù)求得的：\o"點(diǎn)擊圖片看全圖"（18）采用有限差分時無需計算雅可比矩陣。直接問題的性質(zhì)也會影響導(dǎo)數(shù)計算方法的選擇。如果問題是非線性的（如非飽和流），則計算量會非常大，這樣采用有限差分法的成本要比精確方法高（見表4）。然而，有限差分法比較靈活，越來越多的代碼使用了這一方法，如UCODE和PEST（Poeter和Hill，1998；Doherty等，2002）。表4導(dǎo)數(shù)計算的方法比較方法優(yōu)點(diǎn)缺點(diǎn)每次迭代成本線性問題非線性問題直接求導(dǎo)法精確難編程10111101伴隨狀態(tài)法精確難編程211n/a有限差分法容易編程不精確1011010備注：1、如果可以利用這樣一個事實(shí)，即所有參數(shù)的系統(tǒng)矩陣是相同的，則可以減小精確求導(dǎo)和伴隨狀態(tài)法的計算成本；2、根據(jù)100個參數(shù)和20個觀測點(diǎn)假設(shè)問題的模擬運(yùn)行時間，估計高斯－牛頓迭代的成本；3、對于直接非線性問題（即非線性控制方程，如非飽和流問題），假定需要進(jìn)行10次迭代。

六、模型性能如何？根據(jù)模擬流程，在對參數(shù)進(jìn)行估計時，需要評價模型概念和估計參數(shù)的性能。模型質(zhì)量受3個因素影響（Beck，1987）：（1）模型結(jié)構(gòu)的不確定性；（2）模型結(jié)構(gòu)中出現(xiàn)的參數(shù)值的不確定性；（3）系統(tǒng)預(yù)測的不確定性。此處的不確定性不僅指誤差的隨機(jī)波動，也指系統(tǒng)偏差。（一）模型結(jié)構(gòu)的不確定性由于數(shù)據(jù)不完整和不一致，概念模型具有許多不確定性特征。建模人員被迫做簡化假設(shè)。在參數(shù)化、離散化和選擇邊界條件等過程中會引入誤差。解決這些不確定性經(jīng)常需要引入幾個概念模型。可以根據(jù)殘差分布、參數(shù)相關(guān)和參數(shù)的置信區(qū)間等對不同的模型進(jìn)行比較。理想情況下，一個好的模型應(yīng)當(dāng)與觀測結(jié)果、不相關(guān)殘差和合理的參數(shù)值相符?，F(xiàn)在有幾種模型符合這些標(biāo)準(zhǔn)，需要在其中選擇最合適的模型。在時間序列分析和地下水應(yīng)用研究中，已經(jīng)確定了幾種模型選擇標(biāo)準(zhǔn)（Akaike，1974，1977；Rissanen，1978；Schwarz，1978；Hannan，1980；Kashyap，1982）。Carrera和Neuman（1986）應(yīng)用這些標(biāo)準(zhǔn)中的其中4條進(jìn)行綜合測試研究，認(rèn)為Kashyap標(biāo)準(zhǔn)最佳。Carrera等（1990）以及Median和Carrera（1996）也得到類似的研究結(jié)果。然而，研究人員可能會質(zhì)疑只選擇一種模型是否合適。如果其它模型與有效的數(shù)據(jù)相容但不采用這些模型的話，會不符合邏輯。這種質(zhì)疑使得許多模型得到認(rèn)可，并用來確定模型預(yù)測的不確定性（Beven和Binley，1992；Beven和Freer，2001）。（二）校準(zhǔn)過程中的難點(diǎn)優(yōu)化問題中最重要的難題是多解性、不可識別性和不（18）采用有限差分時無需計算雅可比矩陣。直接問題的性質(zhì)也會影響導(dǎo)數(shù)計算方法的選擇。如果問題是非線性的（如非飽和流），則計算量會非常大，這樣采用有限差分法的成本要比精確方法高（見表4）。然而，有限差分法比較靈活，越來越多的代碼使用了這一方法，如UCODE和PEST（Poeter和Hill，1998；Doherty等，2002）。表4導(dǎo)數(shù)計算的方法比較方法優(yōu)點(diǎn)缺點(diǎn)每次迭代成本線性問題非線性問題直接求導(dǎo)法精確難編程10111101伴隨狀態(tài)法精確難編程211n/a有限差分法容易編程不精確1011010備注：1、如果可以利用這樣一個事實(shí)，即所有參數(shù)的系統(tǒng)矩陣是相同的，則可以減小精確求導(dǎo)和伴隨狀態(tài)法的計算成本；2、根據(jù)100個參數(shù)和20個觀測點(diǎn)假設(shè)問題的模擬運(yùn)行時間，估計高斯－牛頓迭代的成本；3、對于直接非線性問題（即非線性控制方程，如非飽和流問題），假定需要進(jìn)行10次迭代。穩(wěn)定性。當(dāng)不止一組參數(shù)能產(chǎn)生正問題的給定解時，會存在不可識別問題；當(dāng)不止一組參數(shù)可以使目標(biāo)函數(shù)最小化時，會存在多解性問題；當(dāng)觀測值的較小變化會引起估計參數(shù)的較大變化時，會存在不穩(wěn)定性問題，但是通常可以根據(jù)求解初始參數(shù)確定。另外，優(yōu)化算法，如Marquardt和共軛梯度法可能會陷入局部極小，不能表現(xiàn)出全局優(yōu)化。Carrera和Neuman（1986）對這些概念進(jìn)行了討論，顯示這些概念之間具有較強(qiáng)的相關(guān)性。事實(shí)上，無論是否是線性問題，都可以用不確定性來刻畫。不穩(wěn)定性和較大的不確定性是不同的概念，但是二者經(jīng)常與延伸的置信區(qū)間有關(guān)，可以通過后驗(yàn)協(xié)方差矩陣的特征值和特征向量來刻畫。該矩陣的特征向量構(gòu)成了一組正交矢量，每個矢量都與一個特征值相關(guān)。與最大特征值相關(guān)的矢量表示具有最大不確定性的參數(shù)的線性組合，而與最小特征值相關(guān)的矢量表示不確定性最小。如果特征值差別極大，則很可能會出現(xiàn)不穩(wěn)定性。估計參數(shù)之間具有某種聯(lián)系，這樣會造成參數(shù)的置信區(qū)間比假如參數(shù)是獨(dú)立時的置信區(qū)間要大很多。解決不穩(wěn)定性和大的不確定性的方法包括：

正則化（調(diào)整參數(shù)），這也正是在式（9）的目標(biāo)函數(shù)中引入Fp的原因。Weiss和Smith（1998）認(rèn)為先驗(yàn)信息是減少不確定性最有效的方法；

減少參數(shù)量，這是“需要估計什么”一節(jié)討論的優(yōu)化方案的動機(jī)；

增加數(shù)據(jù)量和數(shù)據(jù)類型，已在“概念模型：知識和數(shù)據(jù)”一節(jié)討論過；

優(yōu)化觀測方案，可以設(shè)計觀測網(wǎng)和試驗(yàn)來減小模型的不確定性，增加數(shù)據(jù)的可靠性，來區(qū)別對待不同的模型（Knopman和Voss，1989；Usunoff等，1992）。盡管如此，還是不可能得到問題的唯一解。這樣就促使一些研究人員不采用估計，而使用隨機(jī)模擬方法（如，Gómez-Hernández等，1997）。不確定性與所有模擬的總集，而非不確定性的統(tǒng)計特征有關(guān)。還有另一種可供選擇的方法，Yapo等（1998）采用稱為帕累托最優(yōu)（ParetoOptimum）的概念來研究不穩(wěn)定性，指出改進(jìn)目標(biāo)函數(shù)某一組分的參數(shù)矢量集合會造成另一個組分的失效。（三）預(yù)測系統(tǒng)未來行為的難點(diǎn)所有以上的難點(diǎn)會造成模型參數(shù)具有一些不確定性，因此模型預(yù)測也就遺傳了這種不確定性。多尺度概念模型可能會進(jìn)一步造成模型預(yù)測的不確定性。因此，評價預(yù)測的不確定性需要分析參數(shù)和模型不確定性的影響。后者通?？梢愿鶕?jù)模型模擬和評價預(yù)測范圍來分析（如Medina和Carrera，1996）。當(dāng)需要進(jìn)行系統(tǒng)化時，由于所建立的概念模型是非系統(tǒng)的，因此很難進(jìn)行系統(tǒng)化。這樣就需要重點(diǎn)評價參數(shù)不確定性的影響?？梢圆捎脦追N方法來量化預(yù)測的不確定性：線性近似、非線性近似和蒙特卡羅方法。線性方法相對簡單，根據(jù)后驗(yàn)協(xié)方差矩陣進(jìn)行計算：\o"點(diǎn)擊圖片看全圖"（19）式中，Jhp是雅可比矩陣，Ch是測量協(xié)方差矩陣，Cp是參數(shù)協(xié)方差矩陣。如果f是一個模型預(yù)測值（f是參數(shù)p的一個函數(shù)），可以根據(jù)下式確定方差的下限：\o"點(diǎn)擊圖片看全圖"（20)

式中，表示獨(dú)立于參數(shù)不確定性的模型誤差。該方程只是線性近似，因此實(shí)際不確定性可能遠(yuǎn)大于所得到的結(jié)果，這也是需要開發(fā)更為復(fù)雜方法的原因。然而，根據(jù)這一等式可以清楚地量化出不同的因素如何影響預(yù)測的不確定性。一般而言，預(yù)測不確定性與以下條件有關(guān)：（1）預(yù)測對模型參數(shù)和初始條件的靈敏度，顯然，只有當(dāng)預(yù)測對參數(shù)響應(yīng)靈敏時，這些參數(shù)才會引起重視，可以根據(jù)求得。這也是有時會采用靈敏度分析而非誤差分析的原因。（2）模型參數(shù)的不確定性（和/或初始狀態(tài)），可以根據(jù)協(xié)方差矩陣Σp求得。非線性方法應(yīng)用難度更大。Vecchia和Cooley（1987）提出計算非線性置信區(qū)間的方法，得到的結(jié)論如下：（1）線性和非線性置信區(qū)間通常是互相偏移或變化的，而且非線性置信區(qū)間通常要大一些；（2）非線性置信區(qū)間通常比線性置信區(qū)間的變化范圍大；（3）隨著區(qū)間的增大，線性和非線性置信區(qū)間的差異也隨之增大；（4）先驗(yàn)信息可以改變置信區(qū)間的大小。Christensen和Cooley（1999）提出了一種量化模型非線性的方法。目前PEST的預(yù)測分析包括非線性置信區(qū)間。蒙特卡羅方法是計算能力最強(qiáng)的方法，基于不同參數(shù)集的正問題評價。這一方法的主要優(yōu)點(diǎn)是易于理解，可以產(chǎn)生一個概率密度函數(shù)，無需進(jìn)行復(fù)雜的假定；主要問題是很難確定所需的模擬量。實(shí)例之一是普適似然不確定估計（GLUE）方法（Beven和Freer，2001）。七、實(shí)際做什么？（一）應(yīng)用趨勢科技期刊中的文獻(xiàn)往往集中在新方法和新解釋的發(fā)展上，因此不適用于分析應(yīng)用趨勢，只有內(nèi)部報告或?qū)ｎ}研討會才經(jīng)常有關(guān)于應(yīng)用的分析，因此，本節(jié)分析是基于作者的個人觀點(diǎn)，會有一些偏差?？梢枣溄親YPERLINK"/portal.cgi"/portal.cgi的“webofscience”，查詢主要期刊上已出版的所有相關(guān)文章。根據(jù)統(tǒng)計，關(guān)于反演模擬的論文數(shù)量基本穩(wěn)定（其中有5%是關(guān)于地下水模擬研究），而且在過去的13年里，關(guān)于反演模擬研究的論文數(shù)量在緩慢而穩(wěn)定地增加。這些論文覆蓋面很廣，在前文中已經(jīng)引用了一些。然而，很難確定應(yīng)用趨勢，唯一可以確定的是地質(zhì)統(tǒng)計方法反演往往用來分析相對小范圍的問題，如解釋水力試驗(yàn)，而大尺度模型趨向于分析分區(qū)問題?？梢詫⒋蠖鄶?shù)地下水應(yīng)用軟件劃分為這兩類，將在下文作進(jìn)一步討論。（二）地質(zhì)統(tǒng)計學(xué)反演在模擬含水層時，為了確定水力性質(zhì)（特別是導(dǎo)水系數(shù)）的變化，應(yīng)用到了地質(zhì)統(tǒng)計學(xué)方法。然而，可以從兩方面來認(rèn)識這種需要：（1）約束模型參數(shù)；（2）模型預(yù)測。盡管這兩點(diǎn)并不排斥，但很少放在一起使用。在方法上，地質(zhì)統(tǒng)計反演遵循Clifton和Neuman（1982）最初提出的步驟，也就是說，首先提出一個隨機(jī)模型（即不考慮Log-T是否穩(wěn)定，變異函數(shù)和平均值是什么等）；其次，根據(jù)方程（2）和（3），采用可以利用的資料來產(chǎn)生一個前驗(yàn)，log-T的最佳估計及其協(xié)方差；最后，通過最小化方程（9）或利用方程（14），在前兩步分析的基礎(chǔ)上獲得模型參數(shù)的估計。Kitanidis和Vomvoris（1983）、Carrera和Neuman（1986）認(rèn)為需要對統(tǒng)計參數(shù)（誤差方差和相關(guān)距離等）進(jìn)行最優(yōu)估計，因此對這一方法進(jìn)行了改進(jìn)。許多研究人員已經(jīng)認(rèn)識到這些參數(shù)的重要性。事實(shí)上，這也是Zimmerman（1998）對不同地質(zhì)統(tǒng)計學(xué)反演技術(shù)進(jìn)行比較之后得出的結(jié)論之一。野外數(shù)據(jù)的成功應(yīng)用也僅限于相對小尺度的問題，其中包括水力試驗(yàn)解釋（Yeh和Liu，2000；Meier等，2001；Vesselinov等，2001）、井捕獲區(qū)描繪（Vassolo等，1998；Kunstmann等，2002；Harrar等，2003）以及其它應(yīng)用（Barlebo等，2004）。大尺度含水層很少應(yīng)用這一方法，這一事實(shí)說明，很難用地質(zhì)統(tǒng)計學(xué)方法模擬大尺度含水層，同時也表明了地質(zhì)統(tǒng)計學(xué)方法在應(yīng)用過程中的兩個局限性：不能充分考慮地質(zhì)信息，不能重現(xiàn)真實(shí)的變化性。地質(zhì)信息通常用術(shù)語（如沉積類型、導(dǎo)水特征的方向性等）來表達(dá)，這在反演過程中很難說明。這類信息在大尺度范圍內(nèi)可能會很準(zhǔn)確，但在試驗(yàn)尺度內(nèi)，則很少是準(zhǔn)確的，在此過程中，唯一可以確定的是滲透性是變化的，這也是地質(zhì)統(tǒng)計學(xué)所公認(rèn)的。不過，在離散特征得到識別后，可以通過將其作為確定性的特征，而使反演方法得到明顯改進(jìn)，這一點(diǎn)也為Meier等人（2001）所認(rèn)識，他們采用剪切帶應(yīng)力狀態(tài)的相關(guān)信息來確定裂隙導(dǎo)水系數(shù)的各向異性。然而，地質(zhì)統(tǒng)計學(xué)估計的主要問題是不能重現(xiàn)真實(shí)的變化。正如在“如何估計模型參數(shù)p”一節(jié)所討論的，根據(jù)式（9）和（14）可以得到給定數(shù)據(jù)的log-T的條件期望，同樣的，估計不能反映真實(shí)的偏差，因?yàn)槠谕狄呀?jīng)將偏差給過濾掉了。當(dāng)資料充分時（通常需要經(jīng)過詳細(xì)的水力試驗(yàn)），可以從一定程度上確定變化類型（如Meier等，2001）；否則，只有估計協(xié)方差才能提供關(guān)于變異性和導(dǎo)水性質(zhì)好的區(qū)域（溶質(zhì)運(yùn)移通道）或?qū)再|(zhì)差的區(qū)域的信息。解決大尺度地質(zhì)統(tǒng)計學(xué)反演反題的另一個方法是利用地球物理數(shù)據(jù)（Hubbard和Rubin，2000）。Gómez-Hernández及其合作伙伴（Gómez-Hernández等，1997；Capilla等，1998）使用條件模擬對以上問題進(jìn)行了研究。理想情況下，所有這些模擬的平均值都等于條件估計值，但是每一個模擬都可以重現(xiàn)假定的變異性。這樣，在預(yù)測對變異性比較敏感的過程（如污染物運(yùn)移）時，采用模擬方法就適合得多。（三）基于地質(zhì)條件的反演含水層尺度模型和試驗(yàn)尺度模型的主要區(qū)別是：確定變異性和參數(shù)類型的多樣性時對地質(zhì)條件的依賴程度。關(guān)于變異性的地質(zhì)資料一般都不準(zhǔn)確，但是不能忽略。不同地層或同一個地層中不同的地質(zhì)單元，會具有不同的性質(zhì)。這樣，在描繪這些地質(zhì)單元時，就不能忽略這種信息的價值。然而，很難準(zhǔn)確認(rèn)識地質(zhì)單元的邊界條件。因此，當(dāng)認(rèn)識到幾何形狀非常重要時，就要進(jìn)行大量的工作對不同的幾何形狀進(jìn)行研究。這個過程包括確定幾何形狀，并通過可利用的資料進(jìn)行測試，重復(fù)進(jìn)行直到滿意為止。這一過程很單調(diào)，需要模擬人員和地質(zhì)人員的互動。事實(shí)上，這個工作沒有進(jìn)行很好的備案，因此，如果幾年后對模型進(jìn)行修改，就不能確定當(dāng)時為什么會選擇這種模型結(jié)構(gòu)。如果需要同時確定導(dǎo)水系數(shù)和其它類型的參數(shù)類型時，問題就會變得更復(fù)雜。根據(jù)作者的經(jīng)驗(yàn)，補(bǔ)給（補(bǔ)給量和補(bǔ)給時間）、邊界流量、抽水率（抽水率通常都會被低估）以及河流與含水層相互作用等方面經(jīng)常存在不確定性。Cooley等（1986）和Hill（1998）嘗試對這些問題進(jìn)行研究?？梢源_定下面的一些趨勢：

水力傳導(dǎo)系數(shù)和導(dǎo)水系數(shù)的點(diǎn)值易出現(xiàn)誤差。而且，當(dāng)模擬尺度比抽水試驗(yàn)尺度大得多時，這兩個參數(shù)用處不大。需要將這些測量值放在相關(guān)的地質(zhì)條件下進(jìn)行分析；

即使資料不完善，也要在模擬過程中考慮一些主要特征（導(dǎo)水裂隙帶和古河道等）；

含水層活動的許多信息包括在離散事件（洪水和強(qiáng)降雨）中。要充分利用這些資料就需要采用瞬時模擬；

模型校準(zhǔn)很少是唯一的（即不同的模型結(jié)構(gòu)可能都會與硬數(shù)據(jù)擬合很好）。在模型預(yù)測時，需要承認(rèn)這種不確定性。為了減小不確定性，就需要補(bǔ)充資料。這種方法具有幾個缺點(diǎn)：一方面，沒有正確考慮不確定性，模型參數(shù)的協(xié)方差矩陣不僅取決于硬數(shù)據(jù)，也取決于模擬人員的主觀判斷；另一方面，除校準(zhǔn)條件外，無法保證模型的有效性?？傊枰獙@個過程進(jìn)行系統(tǒng)化，這一工作可以通過地質(zhì)統(tǒng)計學(xué)來完成，因此，需要尋找能充分利用定性地質(zhì)資料的地質(zhì)統(tǒng)計學(xué)描述。八、下一步工作通過上述討論，顯而易見，在針對含水層刻畫和管理的地下水研究工