地下水污染物濃度與電阻率映射關系的構建*

趙少樺，董艷輝?，李學蘭，王禮恒

(1 中國科學院地質與地球物理研究所中國科學院頁巖氣與地質工程重點實驗室， 北京 100029； 2 中國科學院地球科學研究院， 北京 100029； 3 中國科學院大學， 北京 100049； 4 中南大學地球科學與信息物理學院， 長沙 410083)

地下水資源在保障社會經(jīng)濟發(fā)展和維持生態(tài)平衡中扮演著不可或缺的角色[1]。然而，由于盲目開發(fā)或保護不當造成的地下水污染問題日益嚴重[2]。為遏制地下水污染態(tài)勢進一步惡化，需要采取有效措施保護地下水安全或進行地下水污染治理，而開展這些工作的前提之一是全面了解地下水質量狀況[3]。目前查明地下水質量現(xiàn)狀的方法可劃分為直接方法和間接方法。直接方法通過進行水文地質鉆探、建立監(jiān)測井獲取地下水樣品進行成分分析。該方法雖然準確性高，但是樣本具有局限性且成本較高。間接方法則主要是通過地球物理手段查明相應的物理場(電場，磁場等)分布狀態(tài)，根據(jù)物理場異常的特征推測地下水質量狀況[4-7]。但是地球物理方法存在多解性的問題，這會造成對異常信息進行分析時所推測的結果可能不符合實際[8]。因此，應用地球物理方法探查地下水污染的關鍵是構建準確、可靠的地球物理信息與污染物濃度之間的映射關系。目前，國內外在此方面的研究并無成熟的方法，多數(shù)為直接利用Archie經(jīng)驗公式將視電阻率反演數(shù)據(jù)轉換為地下水污染物濃度數(shù)據(jù)[9]，但受限于Archie公式適用范圍，采用這種方式所估算的結果誤差較大。此外，一些學者通過電阻率成像技術(ERT，electrical resistance tomography)、數(shù)值實驗以及鹽示蹤實驗來研究視電阻率反演數(shù)據(jù)與濃度之間的關系[10-13]。

針對電阻率信息向污染物濃度信息轉換困難的問題，本文參考Singha和Moysey[8]提出的方法嘗試基于地下水流動與溶質運移模擬和地球物理正反演模擬來獲取含水層的電阻率數(shù)據(jù)與污染物濃度數(shù)據(jù)，并通過不同的統(tǒng)計方法構建其映射關系，從而達到通過視電阻率反演數(shù)據(jù)估計地下水污染物濃度的目的，其中本文將Singha所提方法中的線性擬合改進為多項式擬合。研究結果顯示該方法所估算的地下水污染物濃度相比于直接應用Archie公式能夠更加準確地獲取整體場域信息，同時多項式擬合建立的映射關系估算的結果是最佳的，所以合理地運用該方法將可成為鉆孔直接取樣監(jiān)測方法的有力補充。

1 研究方法與模型實例

1.1 研究方法

地電勘探中在野外獲取到的觀測數(shù)據(jù)為視電阻率數(shù)據(jù)，對觀測數(shù)據(jù)反演計算后可推斷地下物性結構。但是直接通過視電阻率反演剖面推斷出地下水污染物濃度及范圍難度大且不精確。而本研究通過建立視電阻率與污染物濃度的映射關系，應用該映射關系可實現(xiàn)對地下水污染物濃度的定量估計。

研究方法歸結為以下幾個步驟：1)通過建立理想的地下水水流和溶質運移模型模擬計算獲取地下水溶質參考濃度場，該參考濃度場相當于實際中存在的未知地下溶質的濃度場，而亟待解決的問題即是估計該濃度場的濃度。2)利用該濃度場的濃度數(shù)據(jù)建立地電模型進行正演計算得到類似于實際野外觀測數(shù)據(jù)。3)將獲取到的“觀測數(shù)據(jù)”進行反演計算得到視電阻率反演數(shù)據(jù)。4)通過統(tǒng)計擬合的方法建立視電阻率反演數(shù)據(jù)和濃度場的濃度映射關系。5)應用該映射關系根據(jù)視電阻率反演數(shù)據(jù)推算出一個估計的濃度場。通過對比參考濃度場和估計的濃度場，分析該映射關系的可靠性。

本文所建立的濃度和反演電阻率的映射關系并非一個經(jīng)驗公式，而是一種隨空間變化的關系。具體步驟如下(圖1):

圖1 建立濃度和電阻率映射關系的流程圖Fig.1 Flowchart of establishing the relationship between concentration and electrical resistivity

1)濃度場數(shù)據(jù)的獲?。杭俣ǖ叵滤诤畬又械倪\動滿足達西定律，據(jù)此可得地下水滲流連續(xù)性方程[14]：

(1)

式(1)描述地下水在多孔介質中的三維非穩(wěn)定運動，其中：K為滲透系數(shù)，h為水頭高度，SS為儲水率，t為時間，W為單位體積的源匯項。假定溶質在地下水中的遷移滿足菲克定律，可用如下對流彌散方程[15]描述：

(2)

式中：c為溶液中某種組分的濃度，μ為實際平均流速矢量，D為水動力彌散系數(shù)。聯(lián)立上述的兩個方程即可獲取溶質隨地下水遷移的時空分布數(shù)據(jù)。為了獲取多個濃度場數(shù)據(jù)以建立映射關系并進行統(tǒng)計檢驗，研究中利用隨機方法建立多個滲透結構作為模型輸入。

2)電阻率數(shù)據(jù)的獲取：電阻率為電導率的倒數(shù)，而電導率與溶質濃度呈相關性。模型模擬的濃度值在每一個網(wǎng)格中均代表其局部值，恰好能夠利用“局部適用”的Archie經(jīng)驗公式[9]直接將濃度轉換為電阻率：

ρb=F·ρf，

(3)

式中：ρb為電阻率(Ω·m)；ρf為溶液電阻率(Ω·m);F為轉換參數(shù)(無量綱)。

3)地球物理正反演模擬：地球物理正演是指由源的屬性導出到地球物理場的分布屬性；反演則是正演的逆過程，是由地球物理場的分布來推導源的屬性。利用Archie公式轉換而來的電阻率數(shù)據(jù)進行地球物理正演模擬，便可獲取類似于野外地球物理勘探工作中所觀測到的電阻率數(shù)據(jù)，再對正演數(shù)據(jù)進行反演，獲取視電阻率反演數(shù)據(jù)。

4)映射關系的建立：統(tǒng)計出不同模型中相應網(wǎng)格的濃度數(shù)據(jù)與視電阻率反演電阻數(shù)據(jù)，利用線性擬合和多項式擬合針對每個網(wǎng)格中的濃度數(shù)據(jù)和 視電阻率反演數(shù)據(jù)進行擬合。

5)應用上述建立的映射關系根據(jù)視電阻率反演數(shù)據(jù)推算出估計的濃度場，將估計的濃度場和參考濃度場進行對比，從不同角度分析該視電阻率反演數(shù)據(jù)和參考濃度數(shù)據(jù)映射關系的可靠性及該方法的可行性。

1.2 建立映射關系實例

聯(lián)合MODFLOW[16]與MT3DS[17]建立一個二維地下水流動與溶質運移模型：模型長200 m，深(高)45 m，采用1 m×1 m矩形網(wǎng)格對模擬區(qū)進行剖分；水流采用穩(wěn)定流模擬計算，溶質運移模擬期為20 d；模型兩側均為定水頭邊界，水力梯度設置為0.000 1，全場初始質量濃度為50 mg/L。通過改變滲透結構輸入獲取50個不同滲流場的地下水模型。不同滲透系數(shù)通過50次隨機實現(xiàn)，隨機過程平均值為117 m/d，標準差為1.58。在所獲取的50個滲流場投入點源污染物，污染源質量濃度為2 500 mg/L，最后獲取50個不同的濃度場分布。利用濃度值與電導率之間存在換算關系[9](1 mg/L=2 μS/cm)與Archie公式將濃度場分布轉換為電阻率場分布，本文中參數(shù)F取5[9]。

研究中采用高密度電阻率法進行地球物理正、反演計算，采用的軟件為Res2dmod與Res2dinv，該軟件是由M.H.Loke基于二維有限差分法和二維有限元方法開發(fā)的一套公認的高密度電阻率數(shù)據(jù)二維正、反演軟件[18]。正演計算過程中建立一個長200 m，深(高)45 m的地電模型，由于Res2dmod在剖分上不能將每一層的深度設置成相同參數(shù)，故僅在模型的敏感區(qū) (深度在0～15 m) 的剖分與地下水模型中的剖分一致，即在模型0～15 m之間采取1 m×1 m網(wǎng)格剖分，其余區(qū)域的剖分依次增大；水平剖分方面將電極距設置為1 m即可，完成模型的剖分之后，在不同的網(wǎng)格中賦以由相應濃度換算的電阻率值。正演計算之后利用Res2div針對正演計算的數(shù)據(jù)進行反演計算得到視電阻率反演數(shù)據(jù)。

最后，統(tǒng)計分析濃度數(shù)據(jù)和視電阻率反演數(shù)據(jù)，采用不同擬合方式擬合濃度值與視電阻率

反演數(shù)據(jù)之間的關系，進而確定二者映射關系。

2 不同擬合方式的應用結果分析

本研究一共分析5個不同滲透系數(shù)的模型進行不同擬合方式的應用結果。發(fā)現(xiàn)5個模型中有4個模型的應用結果相似，故在這4個模型中選取模型1代表該4個模型的應用結果，剩余的一個模型(模型2)的應用結果與其余4個有明顯的差異。所以本文只討論模型1和模型2的應用結果，其中模型1的滲透系數(shù)為75 m/d，模型2的滲透系數(shù)為28 m/d。

2.1 模型1不同擬合方式的應用結果分析

將MT3DMS模擬的濃度場作為參考濃度場(圖2(a))，Archie公式估算的濃度場(圖2(b))、線性擬合建立的映射關系估算的濃度場(圖2(c))以及多項式擬合建立的映射關系估算的濃度場(圖2(d))進行對比發(fā)現(xiàn)：對于污染暈形態(tài)的刻畫，Archie公式估算的濃度場中，污染暈呈現(xiàn)一個橢圓的形態(tài)，與參考濃度場差異明顯；線性擬合和多項式擬合建立的映射關系的應用結果均十分接近參考濃度場，顯著優(yōu)于直接應用Archie公式，但是二者的差別不明顯。因此，對于污染暈的刻畫，線性擬合與多項式擬合建立的映射關系估算的結果要顯著優(yōu)于Archie公式估算的結果。

圖2 模型1利用不同方式得到的濃度剖面圖Fig.2 The concentration maps of model 1 produced using different ways

對濃度值進行頻率統(tǒng)計：參考濃度場質量濃度多集中在0～100 mg/L(圖3(a))，存在高濃度值；Archie公式估算的質量濃度分布(圖3(b))集中在0～350 mg/L，未出現(xiàn)較高的濃度值(>400 mg/L)，這與參考濃度場明顯不同。線性擬合建立的映射關系估算的濃度分布(圖3(c))與多項式擬合建立的映射關系估算的濃度分布(圖3(d))均接近于參考濃度剖面濃度分布，主要集中在0～100 mg/L，也存在高濃度值。因此對于濃度分布而言，線性擬合與多項式擬合所建立的映射關系估算效果優(yōu)于直接利用Archie公式。

圖3 模型1利用不同方式得到的濃度剖面的濃度分布柱狀圖Fig.3 The concentration distribution histogram of concentration maps of model 1 produced using different ways

為進一步評估不同方法估算的污染物濃度值的可靠性，可將參考濃度作為橫坐標、計算濃度值作為縱坐標，評估散點是否分布在對角線附近，越靠近對角線說明二者越相似，表明該方法越可靠。Archie公式估算的散點偏離對角線，并表現(xiàn)出低濃度散點位于對角線上方，高濃度散點位于對角線下方，即在低濃度區(qū)域會高估濃度值，在高濃度區(qū)域會低估濃度值(圖4(a))；利用線性擬合建立的映射關系和多項式擬合建立的映射關系估算的濃度均集中分布在對角線附近(圖4(b)，4(c))，表明兩種映射關系估算的每個網(wǎng)格的濃度值都接近于參考濃度值。

綜上，在模型1中應用Archie公式法、線性擬合建立的映射關系及多項式擬合建立的映射關系分別估算地下水溶質濃度值，結果表明，線性擬合與多項式擬合建立的映射關系估算效果明顯優(yōu)于Archie公式法估算效果。

2.2 模型2不同擬合方式的應用結果分析

圖4 模型1不同方式計算的濃度與參考濃度散點圖Fig.4 The scatter plot of the reference concentration and the concentration produced using different ways

圖5 模型2利用不同方式得到的濃度剖面圖Fig.5 The concentration maps of model 2 produced using different ways

將模型2 的參考濃度場(圖5(a))與Archie公式估算的濃度場(圖5(b))、線性擬合建立的映射關系估算的濃度場(圖5(c))以及多項式擬合建立的映射關系估算的濃度場(圖5(d))進行對比發(fā)現(xiàn)：對于污染暈形態(tài)的刻畫，Archie公式估算的濃度場中污染暈的形態(tài)與參考濃度場差異明顯；線性擬合建立映射關系估算的污染暈形態(tài)雖然比起Archie公式要接近于參考濃度場，但在低濃度區(qū)域出現(xiàn)負值，這是由于假定場域中所有的點都符合同一個線性函數(shù)來建立映射關系所引起；多項式擬合建立的映射關系的應用結果解決線性擬合中出現(xiàn)負值的問題，也更接近于參考濃度場。因此對于污染暈的刻畫方面，多項式擬合建立的映射關系要優(yōu)于線性擬合建立的映射關系和直接應用Archie公式。

對濃度值進行頻率統(tǒng)計：參考濃度場質量濃度多集中在0～100 mg/L(圖6(a))，同樣存在高濃度值；Archie公式估算的質量濃度分布(圖6(b))在0～200 mg/L，未出現(xiàn)較高的濃度值(>400 mg/L)，這與參考濃度場明顯不符；線性擬合建立的映射關系估算的質量濃度分布(圖6(c))主要集中在0～100 mg/L，也存在高濃度值，但出現(xiàn)負濃度值的分布。多項式擬合建立的映射關系估算的濃度分布(圖6(d))更為接近于參考濃度剖面濃度分布，主要集中在0～100 mg/L，不存在負濃度值。因此進一步可以發(fā)現(xiàn)，線性擬合建立的映射關系在模型2上的應用是存在問題的，多項式擬合建立的映射關系要優(yōu)于線性擬合建立的映射關系和直接應用Archie公式。

同樣將參考濃度作為橫坐標、計算濃度值作為縱坐標，評估不同方法估算的污染物濃度值的可靠性。Archie公式估算的濃度仍然不準確，散點明顯地偏離對角線(圖7(a))；線性擬合建立的映射關系的散點(圖7(b))在低濃度部分出現(xiàn)異常，不集中分布于對角線；多項式擬合建立的映射關系的散點(圖7(c))更為集中對角線附近，即其估計的濃度值會更加接近于參考濃度值。

圖6 模型2利用不同方式得到的濃度剖面的濃度分布柱狀圖Fig.6 The concentration distribution histogram of concentration maps of model 2 produced using different ways

圖7 模型2不同方式計算濃度與參考濃度對比關系Fig.7 The scatter plot of the reference concentration and the concentration produced using different ways

綜上，在模型2中應用Archie公式法、線性擬合建立的映射關系及多項式擬合建立的映射關系分別估算地下水溶質濃度值，結果表明，多項式擬合建立的映射關系估算效果明顯優(yōu)于線性擬合建立的映射關系估算效果的與Archie公式法估算效果。

3 結論與展望

1)聯(lián)合地下水數(shù)值模擬與電阻率成像法建立地下水污染物濃度和電阻率之間的映射關系是可行的。利用所建立的映射關系，可從視電阻率反演剖面估計出地下水污染物的濃度，這較直接應用Archie公式有顯著的改善，因此可作為鉆井取樣直接監(jiān)測的有力補充。

2)地下污染物濃度和電阻率的映射關系具有空間特征，表現(xiàn)在每個網(wǎng)格中對應的映射關系都不相同；利用線性擬合和多項式擬合建立的映射關系在大多數(shù)模型上的應用效果差異不大，但都優(yōu)于直接應用Archie公式計算的結果；對于某些特殊的模型(如模型2)，利用線性擬合映射關系所計算的結果會出現(xiàn)負值，因此，其應用效果要劣于多項式擬合映射關系計算的效果。

3)很多因素都直接影響本文所提方法的準確性，包括地下水數(shù)值模型的準確性、Archie公式中的轉換參數(shù)、電阻率數(shù)據(jù)的準確性等。此外本文提出的方法還處于研究的初步階段，文章中所建立的地下水數(shù)值模型是理想模型，有別于實際的地下水運動與溶質遷移特征，因此需要在實際的場域中進行應用，進一步確定該方法的適用性。