土體網(wǎng)狀裂隙分形插值法生成技術

王媛 ，馮迪 ，陳尚星

(1.河海大學 巖土力學與堤壩工程教育部重點實驗室，江蘇 南京，210098；2.河海大學 隧道與城市軌道工程研究所，江蘇 南京，210098)

由于構造運動、自然氣候和人類活動等作用，自然界的土壤表面可能會發(fā)育著大量的裂隙。裂隙的存在對土體的力學性質(zhì)[1?2]和滲透特性[3]產(chǎn)生了重要的影響，可能會導致土體的工程性質(zhì)變差，裂土地區(qū)的各類工程建設常因此而遭受破壞，造成巨大的經(jīng)濟損失。此外，裂隙土體中水分的流動和溶質(zhì)的運移在農(nóng)業(yè)工程[4]和環(huán)境工程[5]中也得到廣泛關注。土體裂隙網(wǎng)絡的研究[6?9]是開展裂隙土體工作的基礎，只有建立了與實際情況近似的網(wǎng)絡，其力學和滲流分析結果才是可靠的。許多學者借鑒巖體裂隙網(wǎng)絡的生成方法，將裂隙簡化為線段，通過蒙特卡洛方法來生成土體的裂隙網(wǎng)絡[10?11]。然而由此得到的裂隙網(wǎng)絡與實際網(wǎng)絡形態(tài)差別比較大，其主要原因是土體中裂隙形態(tài)與巖體中裂隙不同，這種裂隙形態(tài)的不同會引起其性質(zhì)的不同，例如連通率、滲透張量以及表征單元體積等。巖體中的裂隙多為 1組(或幾組)沿某優(yōu)勢方向的裂隙的隨機分布，而土體中的裂隙多為互相聯(lián)通的多邊形，如果用傳統(tǒng)的生成巖體裂隙的方法，無法模擬土體中裂隙的形態(tài)。本文作者基于分形插值理論，建立了土體裂隙網(wǎng)絡的生成方法，并編制了可視化程序，得到了與實際裂隙網(wǎng)絡形態(tài)較為接近的二維裂隙網(wǎng)絡，為研究裂隙土體的力學性質(zhì)和滲流特性建立了數(shù)值模型。

1 裂隙網(wǎng)絡幾何特征的調(diào)查與統(tǒng)計

1.1 土體裂隙網(wǎng)絡的現(xiàn)場調(diào)查

采用改進的測窗法對裂隙網(wǎng)絡進行測量與調(diào)查[12]。首先，利用數(shù)碼相機進行裂隙圖片的拍攝；然后對拍攝的裂隙圖片進行處理，得到二元化的裂隙灰度圖片，最后將處理過的裂隙圖片導入AutoCAD中，對裂隙網(wǎng)絡進行素描，得到裂隙空間分布圖像，并獲取裂隙網(wǎng)絡的信息。具體步驟如下：

(1) 選擇好場地以后，對拍攝區(qū)域進行選擇，要求拍攝區(qū)域平整不存在凸起和凹陷，同時拍攝區(qū)域內(nèi)的土體顏色要均一、少有雜色物質(zhì)。在拍攝之前要對拍攝區(qū)域進行清理，去除裂隙中的充填雜物，使裂隙充分暴露出來；去除土體表面上的其他雜物，剔除它們可能對影像造成的干擾。在拍攝區(qū)域邊緣放置標度尺，作為確定分析區(qū)域實際尺寸的依據(jù)和比例換算的標準。拍攝時鏡頭主光軸與拍攝面保持垂直，這樣能夠真實反映裂隙的形態(tài)，提高計算準確性。

(2) 將拍攝好的照片轉儲到計算機上，根據(jù)圖像中的標度尺計算確定圖像與實物的比例尺。剪掉圖像中的標度尺部分，挑選圖像質(zhì)量較好并且沒有干擾的區(qū)域作為計算區(qū)域，計算區(qū)域應該盡量處于照片的中央位置。將圖像調(diào)整為灰度模式。對圖像進行增益除噪處理，使裂隙與其他圖像信息最大程度地分離。設置適當?shù)拈y值，將灰度圖轉換為只存在純黑和純白 2種顏色的二元位圖，即實際裂隙網(wǎng)絡圖(如圖1所示)，以此作為土體裂隙幾何學特征研究的素材。

(3) 將處理過的土體裂隙圖像導入到AutoCAD軟件中，對照裂隙圖像，用線段素描得到裂隙網(wǎng)絡幾何特征圖(如圖 2所示)。裂隙網(wǎng)絡圖是進行裂隙分組、統(tǒng)計、計算和裂隙網(wǎng)絡幾何特征分析的基礎。根據(jù)圖像與實物的比例尺，就可以獲取裂隙幾何參數(shù)的真實信息，本文中分析區(qū)域的實際尺寸(長×寬)為 275 mm×220 mm。

圖1 實際裂隙網(wǎng)絡圖Fig.1 Graph of real crack network

圖2 裂隙網(wǎng)絡幾何特征圖Fig.2 Geometric features of crack network

1.2 裂隙網(wǎng)絡的幾何特征統(tǒng)計

通過調(diào)查發(fā)現(xiàn)裂隙土的裂隙大多首尾相接，最終形成由四邊形、五邊形、六邊形及七邊形組成的網(wǎng)狀裂隙網(wǎng)絡。為了能夠在一定程度上反映裂隙的發(fā)育，需要對裂隙網(wǎng)絡的幾何特征進行概化。定義節(jié)點為裂隙與其他裂隙相交的點(如圖 2中圓圈所示)，兩節(jié)點間即為一條裂隙；拐點為裂隙延展方向發(fā)生明顯改變的點(如圖 2中三角形所示)。整個裂隙網(wǎng)絡的幾何結構由節(jié)點、拐點和裂隙段構成，可以通過以下參數(shù)來進行描述：

(1) 節(jié)點密度：單位面積的節(jié)點個數(shù)。

(2) 節(jié)點的度：節(jié)點與其他節(jié)點連接的個數(shù)，它反映了裂隙網(wǎng)絡節(jié)點間的連接程度。如圖 3所示，A點的度為1，B點的度為2，C點的度為3。

(3) 拐點密度：每條裂隙的平均拐點個數(shù)。

(4) 拐點變幅：每條裂隙拐點相對于節(jié)點連線的垂直距離(如圖 4所示)，它反映了每條裂隙的細部結構形態(tài)。

圖3 節(jié)點的度Fig.3 Degree of node

圖4 拐點變幅Fig.4 Amplitude of inflection

(5) 裂隙寬度：是指裂隙張開的寬度，是描述裂隙結構特征的一個重要指標。

裂隙網(wǎng)絡節(jié)點度的統(tǒng)計結果見表 1。可見：裂隙網(wǎng)絡共有節(jié)點41個，其中36個節(jié)點的度為3。裂隙網(wǎng)絡拐點統(tǒng)計結果見表 2?？梢姡毫严毒W(wǎng)絡共有拐點72個，裂隙共有49條，拐點個數(shù)為2的裂隙共有25條，其次為拐點個數(shù)為0和1的裂隙，分別為11條和9條。拐點的最大變幅為17.24 mm，最小變幅為0.06 mm，平均變幅為3.64 mm。裂隙寬度最大為4.12 mm，最小為1.09 mm，平均為2.57 mm。

表1 裂隙網(wǎng)絡節(jié)點統(tǒng)計結果Table 1 Statistical results of crack network nodes

表2 裂隙網(wǎng)絡拐點統(tǒng)計結果Table 2 Statistical results of crack network inflections

2 分形插值法生成土體的裂隙網(wǎng)絡

裂隙土體中的網(wǎng)狀裂隙形態(tài)復雜，裂隙延展方向及隙寬變化較大，如用直線段模擬，就會與實際形態(tài)產(chǎn)生較大的差別。為了更接近實際裂隙網(wǎng)絡，需用另一種方法來進行模擬。易順民等[13?14]在研究膨脹土裂隙時，發(fā)現(xiàn)裂隙面的起伏形態(tài)也具有統(tǒng)計意義上的自相似性符合分形特征，依據(jù)這一結果，本文建立了一種利用分形插值來生成土體裂隙網(wǎng)絡的方法。

2.1 迭代函數(shù)系與分形插值

迭代函數(shù)系(IFS)是一種仿射映射變換方法，其實質(zhì)是通過仿射映射變換，將整體的形態(tài)變換到局部，應用這一特性，人們模擬重構出了大自然中許多具有自相似形態(tài)的圖案[15]。

并滿足如下條件：

則可以確定以下方程組：

對于每個仿射變換ωi有ai，ci，di，ei，fi5 個參數(shù)，如果要確定這5個參數(shù)，就必須增加1個方程或固定其中一個參數(shù)。在這里選取參數(shù)d代表變換中垂直方向變化的參數(shù)，稱其為垂直比例因子。令d為任意一個取定的實數(shù)，聯(lián)解方程可得：

2.2 土體裂隙網(wǎng)絡的生成方法

土體裂隙網(wǎng)絡的生成中采用的參數(shù)是節(jié)點、拐點數(shù)和分維數(shù)。我們通過統(tǒng)計實際土體單位面積的平均節(jié)點數(shù)來生成模擬區(qū)內(nèi)裂隙的節(jié)點數(shù)，然后根據(jù)節(jié)點之間位置關系生成拐點。裂隙土的裂隙符合分形特征，根據(jù)這一結論，裂隙土裂隙的形態(tài)應該具有自相似特征，我們可以根據(jù)節(jié)點、拐點的位置特征，插值出節(jié)點、拐點間的裂隙形態(tài)特征，獲得的裂隙和實際裂隙具有相同的分形特征。

分形插值法生成土體裂隙網(wǎng)絡的過程如下：

(1) 生成域和分析域。考慮到邊界周圍裂隙的節(jié)點可能會在域外，實際裂隙網(wǎng)絡生成域應比分析域大一些，然后用分析區(qū)域邊界去切割生成的裂隙網(wǎng)絡，就可以得到分析域的裂隙網(wǎng)絡圖。

(2) 節(jié)點的生成。根據(jù)節(jié)點密度，計算生成域內(nèi)節(jié)點數(shù)，按均勻分布隨機生成節(jié)點。

(3) 確定節(jié)點間的連接關系。節(jié)點間的連接需要遵從一定的關系，即裂隙之間的拓撲關系。通過對實際裂隙網(wǎng)絡觀察發(fā)現(xiàn)：節(jié)點總是與周圍距離相近的幾個節(jié)點相連。節(jié)點的度k表示節(jié)點與其他節(jié)點的連接數(shù)，因此本文假定每個節(jié)點與其距離最近的k個節(jié)點相連，由此可以確定整個區(qū)域中所有節(jié)點間的連接關系。

(4) 拐點生成。統(tǒng)計節(jié)點間的拐點數(shù)，以平均值作為模擬時節(jié)點間生成的拐點數(shù)。量測每條裂隙拐點的變幅，然后對整個區(qū)域的拐點變幅值進行平均，將其作為裂隙模擬時的變幅值。設兩節(jié)點的坐標分別為(x0,y0), (x1,y1)，則第i個拐點的坐標(xi,yi)為：

式中：m為拐點數(shù)；為拐點變幅平均值；α為該裂隙與x軸的夾角。

(5) 裂隙寬度。測量每條裂隙的寬度，然后對整個區(qū)域的裂隙寬度進行統(tǒng)計，根據(jù)下式確定第i條裂隙的寬度bi：

式中:bmax為最大隙寬；bmin為最小隙寬；random(bmax?bmin)表示在 0～(bmax?bmin)之間的均勻隨機數(shù)。

(6) 分形插值。根據(jù)節(jié)點坐標和節(jié)點間的拐點坐標，由式(5)獲得迭代函數(shù)的參數(shù)，通過下式進行迭代計算：

獲得插值點的坐標，經(jīng)過L次插值迭代后，插值點足夠多便連成線，表現(xiàn)出裂隙的形態(tài)變化和粗細變化。

2.3 裂隙網(wǎng)絡的生成結果與分析

2.3.1 參數(shù)的選取

裂隙網(wǎng)絡生成域為344 mm×275 mm(長×寬)，分析域為275 mm×220 mm(長×寬)，根據(jù)實際裂隙網(wǎng)絡幾何參數(shù)的統(tǒng)計結果可知：相同節(jié)點密度下，生成域內(nèi)節(jié)點數(shù)n為51，節(jié)點的度k平均值為3，節(jié)點間的拐點數(shù)m均值為2，拐點變幅h均值為3.64 mm，最大隙寬為4.12 mm，最小隙寬為1.09 mm，插值迭代次數(shù)L為6。改變節(jié)點坐標隨機數(shù)產(chǎn)生種子可以得到不同的裂隙網(wǎng)絡樣本，如圖5所示。

圖5 不同的裂隙網(wǎng)絡樣本Fig.5 Different samples of crack network

2.3.2 節(jié)點數(shù)n的影響

生成域內(nèi)節(jié)點數(shù)n對分形插值網(wǎng)絡結果的影響如圖6所示。由分維數(shù)計算方法[14]可以得到裂隙網(wǎng)絡的維數(shù)，從圖6可以看出：隨著n的增大，裂隙網(wǎng)絡的裂隙密度和分維數(shù)也在不斷的增大。

圖6 參數(shù)n對分形插值網(wǎng)絡結果的影響Fig.6 Effect of n on crack network generated by fractal interpolation

2.3.3 垂直比例因子d的影響

圖7 參數(shù)d對分形插值網(wǎng)絡結果的影響Fig.7 Effect of d on crack network generated by fractal interpolation

垂直比例因子d對分形插值網(wǎng)絡結果的影響如圖7所示，從圖7可以看出：d雖然是個自由參數(shù)，但其對分形插值的結果影響很大。當d=0時，分形插值退化為線性插值，此時得到的裂隙網(wǎng)絡不能反映土體裂隙的細部形態(tài)；隨著d的不斷增大，土體裂隙的細部起伏形態(tài)表現(xiàn)的越來越明顯，裂隙網(wǎng)絡的分維數(shù)逐漸增大；但是，當d=0.6時，裂隙已呈現(xiàn)出不光滑形態(tài)。因此，在裂隙網(wǎng)絡生成過程中可以通過改變垂直比例因子d來控制土體裂隙網(wǎng)絡的精細程度。

2.3.4 關于裂隙網(wǎng)絡生成結果的討論

在裂隙網(wǎng)絡生成程序中嵌入分維數(shù)計算程序，當?shù)玫降姆中尾逯稻W(wǎng)絡,分維數(shù)等于實際裂隙網(wǎng)絡圖像的分維數(shù)時，插值結束，此時得到的分形插值網(wǎng)絡如圖8所示。

圖8 分形插值網(wǎng)絡(d=0.281，分維數(shù)1.466)Fig.8 Crack network generated by fractal interpolation

目前，土體裂隙網(wǎng)絡的生成主要是通過蒙特卡洛方法來生成[9]，按照裂隙網(wǎng)絡幾何參數(shù)統(tǒng)計結果，得到裂隙網(wǎng)絡如圖9所示。該方法主要存在2個缺點：一是不能反映土體裂隙的細部起伏形態(tài)；二是與土體裂隙多為互相聯(lián)通多邊形的實際形態(tài)不符。將圖8和圖9與土體實際的裂隙網(wǎng)絡(圖1)對比可以發(fā)現(xiàn)：分形插值法得到的裂隙網(wǎng)絡顯然更符合實際形態(tài)。

圖9 蒙特卡洛法生成的裂隙網(wǎng)絡Fig.9 Crack network generated by Monte-Carlo method

3 土體裂隙網(wǎng)絡的表征單元體積

由于土體中裂隙的存在，使得原來的連續(xù)介質(zhì)假定變得不合理，需要用非連續(xù)介質(zhì)的方法進行研究。在現(xiàn)階段的研究中較為成熟并行之有效的方法是等效連續(xù)介質(zhì)模型，采用等效連續(xù)介質(zhì)模型可以將由于裂隙引起的力學或滲流的集中效應等效平均到整個土體中，將裂隙土體模擬為連續(xù)的土體，然后應用經(jīng)典的連續(xù)介質(zhì)理論進行研究分析。在應用等效連續(xù)介質(zhì)進行分析時存在有效性問題，即研究對象必須存在表征單元體積(REV)[18]。表征單元體積(REV)是指當試件的體積大于某一特定值V0時，由試驗獲得的參數(shù)不在隨試件的體積變化而變化，該體積就稱為表征單元體積(如圖10所示)。為此，本文引入表征單元體積的概念來對生成的裂隙網(wǎng)絡質(zhì)量做進一步地評價和討論。

圖10 表征單元體積示意圖Fig.10 Schematic diagram of REV

在多孔介質(zhì)中，一般根據(jù)空隙率來定義表征單元體積[19]。對于裂隙土體，其不連續(xù)性主要是由于裂隙引起的，求表征單元體積的過程就是將裂隙影響效應均一化的過程，因此可以用裂隙率來定義表征單元體積。當裂隙率不隨深度變化而變化時，體裂隙率就等于面裂隙率，面裂隙率定義為：

式中：nc為面裂隙率；Ac為裂隙面積；A為裂隙土體面積。

用一系列不同尺寸的窗口(如圖11所示)去切割裂隙網(wǎng)絡，計算出每個尺寸窗口對應的裂隙率，繪制裂隙率與計算窗口面積之間的關系曲線圖。當裂隙率變化連續(xù)小于誤差范圍且無明顯變化趨勢時，即認為裂隙率穩(wěn)定，確定不出現(xiàn)大的波動的那個最小計算窗口面積，即為表征單元體積所對應的面積。

選取裂隙網(wǎng)絡中間的220 mm×220 mm(長×寬)區(qū)域作為表征單元體積的計算區(qū)域，用步長為 13.75 mm×13.75 mm(長×寬)的16個窗口去切割裂隙網(wǎng)絡，分形插值網(wǎng)絡與實際裂隙網(wǎng)絡的裂隙率隨計算窗口的變化如圖12所示。

圖11 一系列不同尺寸的窗口Fig.11 A series of windows with different scales

圖12 裂隙率隨計算窗口變化圖Fig.12 Variation of crack ratio with different computing windows

從圖12可以看出：實際裂隙網(wǎng)絡的表征單元面積大約為第12個計算窗口，即165 mm×165 mm(長×寬)；模擬得到的分形插值裂隙網(wǎng)絡的表征單元體積大約為第10個計算窗口，即137.5 mm×137.5 mm(長×寬)；從第12個計算窗口開始，兩者的裂隙率趨于相同。結果表明：生成網(wǎng)絡與實際網(wǎng)絡的表征單元體積相差不大，兩者的整體分布特征較為接近。當然，兩者之間還是存在一定的誤差，說明用分形插值法所得到的裂隙網(wǎng)絡還需要進一步的優(yōu)化，使其更加符合實際裂隙發(fā)育的情況。用拐點數(shù)、拐點分布特征及分維值控制生成裂隙還存在一些缺點，以后若從裂隙形態(tài)的頻譜特征去統(tǒng)計和模擬裂隙會改善單條裂隙的相似性，同時建立裂隙間的拓撲關系函數(shù)，會使整個裂隙網(wǎng)絡更趨于合理。

4 結論

(1) 采用改進的測窗法對裂隙網(wǎng)絡進行調(diào)查較為簡便，統(tǒng)計結果表明土體裂隙網(wǎng)絡中節(jié)點的度以3居多，裂隙間的拐點個數(shù)一般為2，其次為0和1。

(2) 分形插值法生成的裂隙網(wǎng)絡能很好地反映出土體裂隙的細部變化特征，基本符合其變化形態(tài)。裂隙網(wǎng)絡的裂隙密度和分維數(shù)隨著節(jié)點密度的增大而增大；垂直比例因子d對結果的影響很大，隨著d的不斷增大，土體裂隙的細部起伏形態(tài)表現(xiàn)得越來越明顯，且裂隙網(wǎng)絡分維數(shù)逐漸增大。

(3) 分形插值生成的裂隙網(wǎng)絡與土體實際裂隙網(wǎng)絡的表征單元體積和穩(wěn)定裂隙率相差不大，表明兩者的整體分布特征相似。

(4) 基于分形插值理論的土體裂隙網(wǎng)絡生成技術切實可行，生成的二維裂隙網(wǎng)絡與土體實際裂隙網(wǎng)絡形態(tài)較為接近。不可忽略裂隙網(wǎng)絡的隨機性和裂隙深度對于土體性質(zhì)的影響，如何全面考慮參數(shù)的隨機性并建立土體三維裂隙網(wǎng)絡模型需要進一步的研究。

[1] 胡卸文, 李群豐, 趙澤三, 等.裂隙性粘土的力學特性[J].巖土工程學報, 1994, 16(4): 81?88.HU Xiewen, LI Qunfeng, ZHAO Zesan, et al.Mechanical properties of fissured clay[J].Chinese Journal of Geotechnical Engineering, 1994, 16(4): 81?88.

[2] 包承鋼.非飽和土的性狀及膨脹土邊坡穩(wěn)定問題[J].巖土工程學報, 2004, 26(1): 1?15.BAO Chenggang.Behavior of unsaturated soil and stability of expansive soil slope[J].Chinese Journal of Geotechnical Engineering, 2004, 26(1): 1?15.

[3] 馬佳.裂土優(yōu)勢流與邊坡穩(wěn)定分析方法[D].武漢: 中國科學院武漢巖土力學研究所, 2007: 6?9.MA Jia.Preferential flow and stability analysis method for fissure clay slopes[D].Wuhan: Institute of Rock and Soil Mechanics.The Chinese Academy of Sciences, 2007: 6?9.

[4] Tuong T P, Cabangon R J, Wopereis M C S.Quantifying flow processes during land soaking of cracked rice soils[J].Soil Science Society of America Journal, 1996, 60(3): 872?879.

[5] Alonso E E, Springman S M, Ng C W W.Monitoring large-scale tests for nuclear waste disposal[J].Geotechnical and Geological Engineering, 2008, 26(6): 817?826.

[6] Velde B.Structure of surface cracks in soil and muds[J].Geoderma, 1999, 93: 101?124.

[7] 盧再華, 陳正漢, 蒲毅彬.膨脹土干濕循環(huán)脹縮裂隙演化的CT試驗研究[J].巖土力學, 2002, 23(4): 417?422.LU Zaihua, CHEN Zhenghan PU Yibin.A CT study on the crack evolution of expansive soil during drying and wetting cycles[J].Rock and Soil Mechanics, 2002, 23(4): 417?422.

[8] Bear J U, Kent T F, Anderson S H.Image analysis and fractal geometry to characterize soil desiccation cracks[J].Geoderma,2009, 154: 153?163.

[9] Li J H, Zhang L M.Geometric parameters and REV of a crack network in soil[J].Computers and Geotechnics, 2010, 37:466?475.

[10] 袁俊平.非飽和膨脹土的裂隙概化模型與邊坡穩(wěn)定研究[D].南京: 河海大學巖土工程研究所, 2003: 22?25.YUAN Junping.Generalized Model of Fissures Distribution &Slope Stability Analysis for Unsaturated Expansive Soils[D].Nanjing: Hohai University.Institute of Geotechnical, 2003:22?25.

[11] Li J H, Zhang L M, Wang Y.Permeability tensor and representative elementary volume of saturated cracked soil[J].Canadian Geotechnical Journal, 2009, 46(8): 928?942.

[12] 李亞萍.甘肅北山花崗巖裂隙幾何特征研究[D].北京: 中國地震局地質(zhì)研究所, 2005: 6?9.LI Yaping.Geometrical Characteristics of Fractures in Granite in the Beishan Area, Gansu Province[D].Beijing: China Earthquake Administration.Institute of Geology, 2005: 6?9.

[13] 易順民, 黎志恒, 張延中.膨脹土裂隙結構的分形特征及其意義[J].巖土工程學報, 1999, 21(3): 294?298.YI Shunmin, LI Zhiheng, ZHANG Yanzhong.The fractal characteristics of fractures in expansion soil and its significance[J].Chinese Journal of Geotechnical Engineering,1999, 21(3): 294?298.

[14] 陳尚星.基于分形理論的土體裂隙網(wǎng)絡研究[D].南京: 河海大學巖土工程研究所, 2006: 44?52.CHEN Shangxing.Reseach on crack network of soil by fractal dimension theory[D].Nanjing: Hohai University.Institute of Geotechnical, 2006: 44?52.

[15] 楊杰.分形插值及其分形維數(shù)研究[J].武漢工業(yè)學院學報,2006, 25(1): 9?11.YANG Jie.Study on fractal interpolated and its fractal dimension[J].Journal of Wuhan polytechnic University, 2006,25(1): 9?11.

[16] 王興元, 梁慶永, 沈立新.真彩色IFS吸引子的計算機構造[J].大連理工大學學報, 2005, 45(1): 142?146.WANG Xingyuan, LIANG Qingyong, SHEN Lixin.True color IFS attractors constructed on computer[J].Journal of Dalian University of Technology, 2005, 45(1): 142?146.

[17] 陳慧琴.分形插值及分形維數(shù)的圖解法[J].江西科學, 2010,28(2): 167?169.CHEN Huiqin.Study on fractal interpolation and its fractal dimension by graphic method solving[J].Jiangxi Science, 2010,28(2): 167?169.

[18] 向文飛, 周創(chuàng)兵.裂隙巖體表征單元體研究進展[J].巖石力學與工程學報, 2005, 24(增2): 5686?5692.XIANG Wenfei, ZHOU Chuangbing.The advances in investigation of representative elementary volume for fractures rock mass[J].Chinese Journal of Rock Mechanics and Engineering, 2005, 24(supp2): 5686?5692.

[19] Bear J.Dynamics of fluids in porous media[M].New York:Elsevier, 2007: 15?27.