劉亞磊，梁 杏，2，朱常坤，李 靜

(1.中國地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北武漢 430074;2.中國地質(zhì)大學(xué)(武漢)生物地質(zhì)與環(huán)境地質(zhì)國家重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430074)

分形理論是由 Mandelbrot［1］提出，用來研究無規(guī)則圖形以及復(fù)雜結(jié)構(gòu)特征的方法。分形是度量上述土壤物理因素的一種手段，分形維數(shù)則是其具體的表現(xiàn)形式［2］。Bartoli［3］和 Tyler等人［4］發(fā)現(xiàn):反映土壤質(zhì)地和結(jié)構(gòu)的物理因素，如容重、粒徑分布、孔隙度及孔隙的連通狀況等，都表現(xiàn)出分形特征，這些物理因素對(duì)土壤的水力參數(shù)存在直接或者間接的影響。因此，國內(nèi)外學(xué)者［4～7］開展了利用土壤相關(guān)物理因素的分形維數(shù)研究土壤水力參數(shù)的工作，并在前人的基礎(chǔ)上不斷提出新的看法與改進(jìn)手段:Tyler等［4］首先將分形理論成功地應(yīng)用于預(yù)測(cè)土壤水分特征曲線，提出了相應(yīng)分形模型并將其擴(kuò)展到三維空間。Kravchenko［8］對(duì)Tyler-Wheatcraft模型進(jìn)行改進(jìn)，提出了一種根據(jù)土壤顆粒分析數(shù)據(jù)分段計(jì)算孔隙表面分形維數(shù)的方法，并分段估計(jì)土壤水分特征曲線，取得了較好的結(jié)果;Toledo等［9］利用分形方法建立水力傳導(dǎo)率模型;在上述基礎(chǔ)上，劉建立等［10］由土壤顆粒質(zhì)量分布曲線計(jì)算得出孔隙表面分形維數(shù)，并利用Burdine模型和Mualem模型預(yù)測(cè)非飽和水力傳導(dǎo)度，預(yù)測(cè)精度較高。分形理論的引入使得間接計(jì)算土壤水力參數(shù)多了一種簡(jiǎn)便而且準(zhǔn)確的方法。

但是，王國梁等［11］認(rèn)為由不同粒徑的土壤顆粒質(zhì)量計(jì)算分形維數(shù)存在不合理的假設(shè)，提出由土壤顆粒體積的大小和數(shù)量來計(jì)算體積分形維數(shù);楊金玲等［12］采用激光衍射法與吸管法實(shí)測(cè)了土壤的粒徑分布，計(jì)算并比較顆粒的質(zhì)量分形維數(shù)和體積分形維數(shù)，發(fā)現(xiàn)二者間存在一定的線性相關(guān)關(guān)系。隨著科技的進(jìn)步，激光衍射可以快速獲取土壤顆粒體積累積曲線，獨(dú)立于顆粒質(zhì)量法，可以得到土壤任意兩粒徑之間的體積百分含量。利用顆粒體積分布曲線計(jì)算的孔隙表面分形維數(shù)，可以避免斯托克斯公式中不合理的假設(shè)，探討此法預(yù)測(cè)土壤水分特征曲線在國內(nèi)相關(guān)研究中并未見報(bào)道。

本文采用土壤顆粒體積分布計(jì)算土壤孔隙表面分形維數(shù)，結(jié)合de Gennes分形模型，預(yù)測(cè)土壤水分特征曲線，通過與實(shí)測(cè)數(shù)據(jù)的對(duì)比分析，探討采用土壤顆粒體積分布曲線計(jì)算表面分形維數(shù)預(yù)測(cè)土壤水分特征曲線的合理性。

1 確定表面分形維數(shù)的原理

土壤的孔隙大小分布狀況對(duì)土壤的水分特征曲線有重要影響，而孔隙表面分形維數(shù)則是描述三維空間內(nèi)土壤孔隙表面不規(guī)則性的一種量度。為了確定此分維，劉建立［10］等基于 Kravchenko 和 Zhang［8］提出的一種由土壤粒徑分布曲線計(jì)算孔隙表面分形維數(shù)的方法，在孔隙體積與孔隙半徑的對(duì)應(yīng)關(guān)系下，建立了土壤顆粒質(zhì)量累積曲線與孔隙表面分形維數(shù)之間的函數(shù)關(guān)系。本文參考劉建立［10］等的推導(dǎo)過程進(jìn)行公式推演，具體步驟如下:

Pachepsky［13］認(rèn)為孔隙體積增量 dVp(≤r)和 dVp(＞r)與孔隙半徑r之間存在如下關(guān)系:

而Perrier［14］等根據(jù)孔隙體積與孔隙半徑之間的關(guān)系提出孔隙體積增量的一種變化形式:

式中:E——?dú)W氏幾何中的拓?fù)渚S數(shù)，E=3;

β——常數(shù);

DS——孔隙表面分形維數(shù)。

假設(shè)土壤孔隙為圓柱狀，土壤孔隙體積可表示為:

此處認(rèn)為土壤孔隙最小半徑為0，結(jié)合式(2)、(3)可求出孔隙半徑小于等于r的孔隙長度L(≤r)為:

Tyler等［5］認(rèn)為N個(gè)等價(jià)半徑為R的固體顆粒構(gòu)成的毛細(xì)管長度為l(R):

式中:R——顆粒粒徑;

N——顆粒粒徑為R時(shí)的顆粒數(shù)量;

D——彎曲毛細(xì)管的孔隙分形維數(shù)，Kravchenko等［8］認(rèn)為 D=DS-1。

Kravchenko等［8］假設(shè)土粒密度 ρs不變，認(rèn)為個(gè)數(shù)N可以由半徑為R的顆粒的質(zhì)量W(R)計(jì)算出來，即:

本文采用土壤顆粒體積V(R)計(jì)算相應(yīng)粒徑級(jí)別下的顆粒數(shù)量，避免了上式中關(guān)于土粒密度不變的假設(shè)，得到個(gè)數(shù)N的計(jì)算式:

式中:S——土壤顆粒形狀的系數(shù);

e——孔隙比，即孔隙體積與顆粒體積之比:

r是顆粒半徑為R時(shí)組成的土壤孔隙半徑。假設(shè)R最小為0，對(duì)式(5)進(jìn)行積分可得由半徑小于等于R的顆粒組成的孔隙半徑小于等于r的孔隙累計(jì)長度:

綜合式(4)、(8)和(9)，可得半徑為R的顆粒累積V(R)關(guān)于土壤孔隙表面分形維數(shù)之間的關(guān)系式:

對(duì)上式進(jìn)行積分則可得半徑小于等于R的土壤顆粒的累積體積，假設(shè)最小半徑為0:

式中:c——常數(shù)項(xiàng)。

對(duì)式(11)兩邊取對(duì)數(shù)得:

式(12)與Kravchenko和劉建立等人建立的方程是相似的，只是將土壤累積質(zhì)量W(≤R)替換為累積體積V(≤R)，此式的優(yōu)點(diǎn)是不必假設(shè)各級(jí)土壤顆粒密度不變，累積體積則可由激光衍射法測(cè)試得出。本文中所用顆粒半徑為相應(yīng)粒徑分級(jí)的上限值與下限值的算術(shù)平均值，通過上式與土壤顆粒體積的累積曲線，可計(jì)算土壤的表面分形維數(shù)Ds。

2 土壤水分特征曲線的分形模型

利用土壤的分形維數(shù)預(yù)測(cè)土壤水分特征曲線，眾多學(xué)者按照不同的理論提出多種預(yù)測(cè)模型。

de Gennes［16］根據(jù)土壤孔隙表面是由自相似性的孔洞互相嵌套組成或者團(tuán)聚體連接而成的兩種模式下分別導(dǎo)出的預(yù)測(cè)模型如下:

式中:Ψ——負(fù)壓絕對(duì)值;

Ψa——進(jìn)氣值;

θs——飽和含水量。

Perrier等［14］根據(jù)孔隙體積增量的變化形式，即式(2)，提出一種水分特征曲線的分形模型:

式中:V0——孔隙半徑為0時(shí)的孔隙體積;

V——總體積。

劉建立等［10］人假設(shè) V0/V=θs- θr(θr為殘余含水量)，代入式(14)即與Brooks-Corey模型一樣:

由于 θr對(duì)上式的影響不敏感，θr≈0，代入式(15)，其形式與式(13)完全相同。因此，選取式(13)作為本文預(yù)測(cè)土壤水分特征曲線的分形模型，采用Brooks-Corey模型作為擬合實(shí)測(cè)曲線的參考模型。

3 應(yīng)用實(shí)例

3．1 土壤樣品與粒徑分級(jí)

本文中所用到的材料為華北平原辛集地區(qū)土樣(表1)。

采用張力計(jì)法進(jìn)行土壤水分特征曲線的實(shí)測(cè);通過激光粒度儀對(duì)測(cè)試土樣進(jìn)行顆粒分析，獲取土壤顆粒體積累積曲線。

本文所用樣品主要為細(xì)粒土，土壤顆粒最大半徑小于 1mm(表 1)。結(jié)合 Kravchenko［8］與郭中領(lǐng)［16］等人的研究成果，基于美國土壤質(zhì)地分類系統(tǒng)，進(jìn)行以下幾種粒徑分級(jí)(表2)。對(duì)不同的分級(jí)方法分別進(jìn)行表面分形維數(shù)的計(jì)算。

3．2 表面分形維數(shù)計(jì)算

由所獲取的顆粒體積累積曲線，分別按照表2所列的四種不同分級(jí)方法，以lg(Rave)、lg［V(≤R)］為X、Y軸(圖1)，基于最小二乘法的線性回歸，得出二者關(guān)系曲線的斜率K。

表1 土壤樣品顆粒分級(jí)情況Table 1 Particle-size class of samples

表2 不同粒徑分級(jí)方法Table 2 Different methods for grading particle-size class

結(jié)合式(12)有K與Ds的關(guān)系:

對(duì)式(16)進(jìn)行求解，由于為三維空間，取兩個(gè)計(jì)算結(jié)果中較大的數(shù)值為DS，結(jié)果見表3。

由表3可以發(fā)現(xiàn)，Ds處于2.37～2.91之間。不同分級(jí)方法中，顆粒分級(jí)界限的最大半徑越小，累積體積對(duì)數(shù)與粒徑半徑對(duì)數(shù)之間的相關(guān)系數(shù)越大，即:D1＞D2＞傳統(tǒng)7級(jí)，說明利用顆粒體積分布計(jì)算表面分形維數(shù)同樣存在無標(biāo)度區(qū)間。計(jì)算結(jié)果顯示隨著樣品粘粒含量的減少，質(zhì)地越粗，表面分形維數(shù)逐漸降低的規(guī)律，與其他學(xué)者計(jì)算的分形維數(shù)規(guī)律是一致的。四種分級(jí)計(jì)算結(jié)果相比較，D3計(jì)算得到的表面分形維數(shù)處于D2與傳統(tǒng)7級(jí)之間，且與D2計(jì)算的結(jié)果間差距較小，說明增大分級(jí)密度對(duì)計(jì)算表面分形維數(shù)影響不大;采用不同粒徑分級(jí)方法計(jì)算會(huì)影響表面分形維數(shù)的大小，但不影響其反映的物理意義。

圖1 不同分級(jí)方法下K的擬合圖(部分)Fig．1 Slope coefficients of different particle-size grading methods

表3 計(jì)算得DS值與相關(guān)系數(shù)Table 3 DSand related coefficient of samples

3．3 實(shí)驗(yàn)法獲取SWRC

本節(jié)選取辛集1～4號(hào)樣品進(jìn)行土壤水分特征曲線的參數(shù)擬合與預(yù)測(cè)分析。首先將1～4號(hào)原狀樣通過張力計(jì)與稱重方式獲取試樣的土壤水分特征曲線，通過RETC軟件中的Brooks-Corey模型進(jìn)行擬合求參，獲取進(jìn)氣值Ψa，并與預(yù)測(cè)擬合結(jié)果對(duì)比。表4為所選4個(gè)樣品的參數(shù)擬合結(jié)果。

表4 土壤水分特征曲線擬合參數(shù)Table 4 The fitted parameters of SWRC

3．4 利用分形維數(shù)進(jìn)行SWRC的預(yù)測(cè)

分別將D1和D2兩種分級(jí)計(jì)算的表面分形維數(shù)(表3)與表4中所列參數(shù)代入de Gennes預(yù)測(cè)模型(式14)，繪制預(yù)測(cè)土壤水分特征曲線(圖2);并采用均方根誤差RMSE作為衡量預(yù)測(cè)模型與實(shí)測(cè)值之間準(zhǔn)確性的計(jì)算標(biāo)準(zhǔn)(以含水量作為評(píng)價(jià)值)，見式(17):

式中:Xobs——實(shí)測(cè)值;

Xi——預(yù)測(cè)值，i=1，2，3，…，N。

由圖2可以很直觀地發(fā)現(xiàn)，采用D1分級(jí)方法計(jì)算出的表面分形維數(shù)預(yù)測(cè)的土壤水分特征曲線與實(shí)測(cè)值擬合得好，而通過D2進(jìn)行的計(jì)算結(jié)果相對(duì)較差(表5)。

圖2 兩種分形維數(shù)下土壤水分特征曲線的預(yù)測(cè)值與實(shí)測(cè)值Fig．2 Prediction and observed points of SWRC by using two kinds of surface fractal dimension

表5 預(yù)測(cè)與實(shí)測(cè)的RMSE值Table 5 The RMSE of prediction and observed points

對(duì)D1分級(jí)方法的預(yù)測(cè)結(jié)果進(jìn)行分析:最大的均方根誤差出現(xiàn)在4號(hào)土樣，為0.0105cm3/cm3，而其他三個(gè)土樣的RMSE均小于0.01cm3/cm3，也就是說，采用此種分級(jí)方法計(jì)算得出的表面分形維數(shù)對(duì)預(yù)測(cè)土壤水分特征曲線具有較高的適用性，盡管3號(hào)土樣在高吸力段的預(yù)測(cè)值與實(shí)測(cè)值略有偏差，但其他試樣的預(yù)測(cè)結(jié)果在整個(gè)吸力范圍內(nèi)與實(shí)測(cè)值較匹配。D2的預(yù)測(cè)結(jié)果與實(shí)測(cè)值只有在低吸力段匹配較好，隨著吸力的增加，預(yù)測(cè)值與實(shí)測(cè)值出現(xiàn)較大偏差，土樣預(yù)測(cè)值與實(shí)測(cè)值之間的RMSE≥2.11E-02cm3/cm3。

4 結(jié)論

(1)利用土壤顆粒體積累積曲線計(jì)算表面分形維數(shù)，避免了土壤顆粒密度不變的假設(shè)，改善了計(jì)算過程中的缺陷。

(2)相同試樣由不同顆粒分級(jí)方法計(jì)算得出的表面分形維數(shù)大小不同;相同分級(jí)方法下，試樣的表面分形維數(shù)隨著土壤粘粒的減少而變小，符合一般規(guī)律;即顆粒分級(jí)方法的不同不影響其反映的物理意義。

(3)按照D1計(jì)算的表面分形維數(shù)預(yù)測(cè)土樣吸力大于進(jìn)氣值的曲線趨勢(shì)精度高，誤差小，這是由于土壤負(fù)壓超過進(jìn)氣值后，土壤水分特征曲線主要受由半徑小的顆粒組成的細(xì)小孔隙控制，從側(cè)面說明土壤孔隙對(duì)水力參數(shù)的影響。因此，此法可以用于預(yù)測(cè)更高吸力值情況下的曲線變化趨勢(shì);而D2的情況，其預(yù)測(cè)的效果較差，誤差偏大，但可以滿足田間粗略估計(jì)的需要;應(yīng)用土壤顆粒體積累積曲線在無標(biāo)度區(qū)間內(nèi)計(jì)算的表面分形維數(shù)更適合以de Gennes模型來預(yù)測(cè)細(xì)粒土的土壤水分特征曲線。

(4)筆者認(rèn)為，利用顆粒體積累積曲線計(jì)算表面分形維數(shù)的方法比較適用于相對(duì)均質(zhì)的土樣，可以避免由于參與顆粒分析的樣品量少而不具有代表性的不利影響。在無標(biāo)度區(qū)間計(jì)算的表面分形維數(shù)預(yù)測(cè)結(jié)果較好，這可能是由于大于0.1mm的土壤顆粒組成的孔隙并不具有很強(qiáng)的自相似性從而影響分形維數(shù)的計(jì)算。因此，利用土壤顆粒體積累積曲線按照無標(biāo)度區(qū)間計(jì)算表面分形維數(shù)，從而預(yù)測(cè)土壤水分特征曲線是可行的，此法具有較明確的物理意義。

［1］ Mandelbrot B B.The fractal geometry of nature［M］.San Francisco:W H Freeman，1983.

［2］ 任權(quán)，王家鼎，袁中夏，等.高速鐵路地基黃土微結(jié)構(gòu)的分形研究［J］.水文地質(zhì)工程地質(zhì)，2007，34(6):76-82.［Ren Q，Wang J D，Yuan Z X，et al.Fractal study on loess microstructure of one high-speed railway foundation［J］.Hydrogeology ＆ Engineering Geology，2007，34(6):76-82.(in Chinese)］

［3］ Bartoli F.Structure and self-similarity in silty and sandy:the fractal approach［J］.Soil Sci，1991，42:167-185.

［4］ Tyler S W，Wheatcraft S W.Fractal scaling of soil particle-size distribution:Analysisand limitations［J］.Soil Sci.Soc.Am.1992，56:362-369.

［5］ Tyler S W，Wheatcraft S W.Application of fractal mathematics to soil water retention estimation［J］.Soil Sci Soc Am，1989，53:987-99.

［6］ 楊培嶺，羅遠(yuǎn)培，石元春.用粒徑的重量分布表征的土壤分形特征［J］.科學(xué)通報(bào)，1993，38(20):1896-1899.［YANG P L，LUO Y P，SHI Y C.Fractal feature ofsoilon expression by weight distribution of particle size［J］.Chinese Science Bulletin，1993， 38(20):1896-1899.(in Chinese)］

［7］ 楊建，陳家軍，楊周喜，等.松散砂?？紫督Y(jié)構(gòu)、孔隙分形特征及滲透率研究［J］.水文地質(zhì)工程地質(zhì)，2008，35(3):93-98.［YANG J，CHEN J J，YANG Z X，et al.A study of pore structure，pore fractal feature and permeability of unconsolidated sand［J］.Hydrogeology＆ Engineering geology，2008，35(3):93-98.(in Chinese)］

［8］ Kravchenko A，Zhang R D.Estimating the soil water retention from particle 2 size distributions:A fractal approach［J］.Soil Science，1998，163:171-179.

［9］ Toledo P G，Novy R A，Davis H T，et al.Hydraulic conductivity of porous media at low water content［J］.Soil Science Society of America Journal，1990，54:673-679.

［10］ 劉建立，徐紹輝，劉慧，等.確定田間土壤水力傳導(dǎo)率的分形方法［J］.水科學(xué)進(jìn)展，2003，14(4):464-469.［LIU J L，XU S H，LIU H，et al.Fractal method for estimating field hydraulic conductivity function from soilparticle-size distribution［J］.Advances In Water Science，2003，14(4):464-469.(in Chinese)］

［11］ 王國梁，周生路，趙其國.土壤顆粒的體積分形維數(shù)及其在土地利用中的應(yīng)用［J］.土壤學(xué)報(bào)，2005，42(4):545-550.［WANG G L，ZHOU S L，ZHAO Q G.Volume fractal dimension of soil particles and its applications to land use［J］.Acta Pedologica Sinica，2005，42(4):545-550.(in Chinese)］

［12］ 楊金玲，李德成，張甘霖，等.土壤顆粒粒徑分布質(zhì)量分形維數(shù)和體積分形維數(shù)的對(duì)比［J］.土壤學(xué)報(bào)，2008，45(3):415-419.［YANG J L，LI D C，ZHANG G L，et al.Comparison of mass and volume fractal dimensions of soil particle-size distributions［J］.Acta Pedologica Sinica，2008，45(3):415-419.(in Chinese)］

［13］ Pachepsky Y A，Shcherbakov R A，Korsunskaya L P.Scaling of soil water retention using a fractal model［J］.Soil Science，1995，159:99-104.

［14］ Perrier E，Rieu M，Sposito G，et al.Models of water retention curve for soils with a fractal pore size distribution［J］.Water Resources Research，1996，32:3025-3031.

［15］ de Gennes P G.Partial filling of a fractal structure by a wetting fluid［C］//Adler D.Physics of Disordered Materials.New York:Plenum Press，1985:227-241.

［16］ 郭中領(lǐng)，符素華，張學(xué)會(huì)，等.土壤粒徑重量分布分形特征的無標(biāo)度區(qū)間［J］.土壤通報(bào)，2010，41(3):537-541.［GUO Z L，F(xiàn)U S H，ZHANG X H，et al.Scale-ree Domain of Fractal Characteristic of the Soil Particle-size Distribution［J］.Chinese Journal of Soil Science，2010，41(3):537-541.(in Chinese)］