張亦弛，呂明明，關(guān) 濤，王佳俊，余 佳，任炳昱

（天津大學(xué) 水利工程仿真與安全國家重點實驗室，天津 300072）

1 研究背景

對于水電工程壩基裂隙巖體，探明其內(nèi)部的裂隙發(fā)育情況能夠為壩基工程地質(zhì)分析與決策提供關(guān)鍵依據(jù)，對于保障水電工程的安全穩(wěn)定具有重要的理論與現(xiàn)實意義［1-3］。裂隙是指巖石受到構(gòu)造變形的破裂作用或物理成巖作用形成的沒有明顯位移的面狀不連續(xù)體［4］，在巖體內(nèi)部常呈現(xiàn)網(wǎng)狀分布。離散裂隙網(wǎng)絡(luò)（Discrete Fracture Networks，DFN）模型是由Baecher 等［5］提出的一種研究裂隙空間展布的有效手段，該方法將每一個裂隙表達為一個具有厚度的圓盤模型?；趶膸r體表面露頭與內(nèi)部鉆孔中觀測的裂隙先驗信息，通過示性點過程進行裂隙網(wǎng)絡(luò)隨機建模，其中點過程基于穩(wěn)態(tài)泊松過程確定裂隙中心點的位置［6］，示性過程確定裂隙的屬性，即圓盤裂隙的幾何參數(shù)［4］。

三維離散裂隙網(wǎng)絡(luò)模型本質(zhì)上是對巖體內(nèi)裂隙的傾向、傾角、直徑、開度等屬性特征聯(lián)合分布的幾何表達［4，6］，因此，裂隙網(wǎng)絡(luò)的多參數(shù)模擬是目前的研究趨勢［6-8］，其關(guān)鍵問題在于對裂隙的幾何參數(shù)之間的多維聯(lián)合分布進行概率密度估計與采樣?，F(xiàn)有的DFN 建模方法通常假設(shè)裂隙的幾何參數(shù)之間互相獨立，并采用典型分布對每個參數(shù)的邊緣分布進行概率密度估計（Probability Density Estima?tion）與采樣（Sampling），例如：描述直徑與跡長的Gamma 分布與對數(shù)正態(tài)分布、描述開度的負指數(shù)分布與對數(shù)正態(tài)分布、描述走向的Von-Mises 分布等［4，6，9］。Mendoza-Torres 等［10］的研究指出現(xiàn)有DFN建模研究通常采用的單參數(shù)模擬方法忽略了裂隙幾何參數(shù)之間可能存在的相關(guān)性［10］，當(dāng)實測裂隙參數(shù)之間不獨立時，分別從各邊緣分布中抽樣的結(jié)果不一定會服從聯(lián)合分布，相反，直接對聯(lián)合分布進行估計的方法自然地考慮了參數(shù)之間的相關(guān)性［10］。然而，傳統(tǒng)的DFN 建模方法僅針對傾向與傾角參數(shù)考慮了聯(lián)合分布估計的問題，且仍然采用典型分布假設(shè)，包括：Fisher 分布、雙正態(tài)分布、Bing?ham 分布、Kent 分布等［6，11］。事實上，由于巖體內(nèi)部裂隙發(fā)育的內(nèi)在規(guī)律十分復(fù)雜，工程實測數(shù)據(jù)中三個以上裂隙幾何參數(shù)的聯(lián)合分布一般不符合上述幾種典型分布中的任何一種［10］，此時基于典型分布的估計方法缺乏足夠的擬合能力［10］，且嚴(yán)格的分布類型假設(shè)可能帶來偏倚與誤差［10-11］。綜上，目前的DFN 建模方法面臨著難以有效擬合多維幾何參數(shù)的聯(lián)合分布的難題，亟待提出一種能夠靈活且準(zhǔn)確地對實測裂隙的多參數(shù)聯(lián)合分布進行概率密度估計與采樣的方法，從而實現(xiàn)DFN 的多參數(shù)模擬。

深度生成模型屬于一種深度學(xué)習(xí)模型［12-13］，其優(yōu)勢在于能夠?qū)W習(xí)分布類型未知的高維聯(lián)合分布并從中采集高質(zhì)量的樣本，能夠克服傳統(tǒng)統(tǒng)計學(xué)模型難以擴展到高維問題的不足［14-15］，對于突破DFN多參數(shù)模擬的瓶頸具有較大的潛力。深度生成模型近幾年開始被應(yīng)用于建立更加真實的大尺度地質(zhì)構(gòu)造三維模型［16］，但尚未被引入DFN 建模的研究中［17-18］。

2017年，Papamakarios 等［19］與Kingma 等［20］提出了一種新的深度生成模型——自回歸流（Autore?gressive Flow），該模型通過建立易解的概率密度（Tractable Density）實現(xiàn)對于復(fù)雜分布下對數(shù)似然函數(shù)值的精確計算［21］，克服了其他基于變分推斷（Variational Inference）的深度生成模型需要對概率密度進行近似估計的不足［12］，具備能夠?qū)崿F(xiàn)更加準(zhǔn)確的概率密度極大似然估計的優(yōu)勢［14-15］，因而成為了目前深度生成模型領(lǐng)域的研究熱點，應(yīng)用于概率密度估計［19］、變分推斷［20］、圖像生成［22］、語音合成［23］、文本生成［24］等方向的研究。本文將自回歸流模型引入DFN 建模中。裂隙網(wǎng)絡(luò)的幾何參數(shù)通常具備優(yōu)勢分組的現(xiàn)象，即巖體在多次構(gòu)造作用下產(chǎn)生的裂隙網(wǎng)絡(luò)的幾何參數(shù)聚集于幾個峰值附近［6，8］，導(dǎo)致其幾何參數(shù)的分布呈現(xiàn)出多峰的特點。然而，自回歸流模型的初始分布通常采用一個單峰的高斯分布，對于多峰分布的估計能力有待提高。因此，本文提出基于改進自回歸流模型的裂隙多參數(shù)模擬方法，將自回歸流模型的標(biāo)準(zhǔn)化特征空間的高斯分布改進為高斯混合分布并結(jié)合DensityPeak 聚類算法［25］提出了密度峰值聚類自回歸流（Density Peak Clustering Autoregressive Flow，DPCAF）模型。其中，密度峰值聚類（Density Peak Clustering）算法于2014年發(fā)表于Science［25］，是目前最受關(guān)注的聚類算法之一，具有能夠從任意維度任意形狀的數(shù)據(jù)中快速確定簇的數(shù)量與聚類中心的優(yōu)勢［26］，適用于裂隙幾何參數(shù)的聚類問題。本文提出DPCAF 模型以期能夠有效擬合裂隙幾何參數(shù)的多維多峰聯(lián)合分布，有效提高裂隙網(wǎng)絡(luò)多參數(shù)模擬的精度，并使得DFN 模型更加真實地反映巖體內(nèi)部的裂隙分布。

2 研究框架

提出的基于密度峰值聚類自回歸流（Density Peak Clustering Autoregressive Flow，DPCAF）模型的三維裂隙網(wǎng)絡(luò)建模方法的研究框架如圖1所示，具體如下：

圖1 研究框架

（1）建立壩基裂隙數(shù)據(jù)集。從壩基裂隙巖體的鉆孔圖像數(shù)據(jù)中提取出巖體內(nèi)部裂隙的傾向、傾角、開度三個維度的幾何信息［6，27-28］，由于鉆孔圖像中未包含裂隙直徑信息，從建基面測窗數(shù)據(jù)中統(tǒng)計裂隙跡長的分布，用于輔助推斷裂隙直徑參數(shù)；

（2）建立DPCAF 深度生成模型。將自回歸流模型標(biāo)準(zhǔn)化特征空間的基礎(chǔ)分布由標(biāo)準(zhǔn)高斯分布改進為高斯混合分布，并結(jié)合DensityPeak 聚類算法確定裂隙分組數(shù)量與聚類中心，建立改進的自回歸流模型DPCAF，利用產(chǎn)狀與開度的實測數(shù)據(jù)集訓(xùn)練DPCAF 模型參數(shù)；

（3）多參數(shù)聯(lián)合分布擬合。利用訓(xùn)練完成的DPCAF 模型對傾向、傾角、開度參數(shù)的聯(lián)合分布進行概率密度估計，同時實現(xiàn)多參數(shù)優(yōu)勢分組；

（4）DFN 模型建立?；诹严抖鄥?shù)優(yōu)勢分組結(jié)果，聯(lián)合測窗數(shù)據(jù)與鉆孔圖像數(shù)據(jù)對每組裂隙的體密度（單位體積內(nèi)的裂隙數(shù)量）進行估計，并通過泊松過程得到各組裂隙中每一個裂隙面的中心點坐標(biāo)參數(shù)；利用DPCAF 模型對傾向、傾角、開度參數(shù)隨機采樣；基于測窗數(shù)據(jù)中的裂隙跡長分布優(yōu)化裂隙直徑參數(shù)；最后，基于圓盤裂隙模型，建立三維離散裂隙網(wǎng)絡(luò)。

3 基于密度峰值聚類自回歸流的DFN 多參數(shù)模擬方法

3.1 密度峰值聚類自回歸流模型由于裂隙網(wǎng)絡(luò)的幾何參數(shù)普遍存在優(yōu)勢分組的現(xiàn)象，其分布呈現(xiàn)出多峰的特點。然而，自回歸流模型的初始分布采用一個單峰的高斯分布，估計多峰的概率密度的能力有待提高。針對以上不足，本文提出的DPCAF 模型，將自回歸流模型的標(biāo)準(zhǔn)化空間的高斯分布改進為一個由均值、方差和各分支的混合比例所參數(shù)化的高斯混合分布，并結(jié)合DensityPeak 算法，以無監(jiān)督的形式實現(xiàn)裂隙多參數(shù)優(yōu)勢分組與多參數(shù)聯(lián)合概率密度估計的多任務(wù)學(xué)習(xí)，提高對實測裂隙數(shù)據(jù)多峰分布的擬合精度。

自回歸流模型屬于一種標(biāo)準(zhǔn)化流（Normalizing Flow）［14-15］，如圖2所示，其核心思想是將一系列簡單的可逆變換函數(shù)作為組件以“流”的形式依次嵌套為一個可逆的復(fù)合函數(shù)，構(gòu)造出一個具有更強擬合能力的非線性雙射，通過在實際數(shù)據(jù)的目標(biāo)分布與已知的基礎(chǔ)分布之間互相轉(zhuǎn)換實現(xiàn)多維聯(lián)合分布的估計。

圖2 自回歸流模型網(wǎng)絡(luò)結(jié)構(gòu)［19］

令隨機變量Z ∈RD服從已知的基礎(chǔ)分布，隨機變量服從目標(biāo)分布。利用N 個雙射函數(shù) f1，f2，…，fN的復(fù)合函數(shù)f 構(gòu)建可逆變換，將基礎(chǔ)分布變換為一個未知的目標(biāo)分布，即，其逆變換為雙射函數(shù) g1，g2，…，gN的復(fù)合函數(shù)g，即。 根據(jù)Change-of-Variables 公式［14］，令J 為雅克比行列式算子，則目標(biāo)分布的概率密度函數(shù)如下：

變換f 的方向稱作生成方向（Generative Direction），是進行采樣時的數(shù)據(jù)流向，逆變換g 的方向稱作標(biāo)準(zhǔn)化方向（Normalizing Direction），通常將復(fù)雜不規(guī)則的目標(biāo)分布變換為多元標(biāo)準(zhǔn)高斯分布［19-20］，是評估模型概率密度值的數(shù)據(jù)流向，如圖2所示。為了保證模型訓(xùn)練和預(yù)測的計算效率，要求可逆變換f 具有易解的雅克比行列式（Tractable Jacobian Determinant）。因此，引入另一種深度生成模型——自回歸模型（Auto-Regressive Model，AR）［29］，該模型基于概率鏈?zhǔn)椒▌t將聯(lián)合概率密度分解為每一個維數(shù)上的條件概率的乘積［29］，并逐個對每一維數(shù)據(jù)的條件概率分布進行“自回歸”計算。以自回歸模型作為雙射函數(shù) fi的標(biāo)準(zhǔn)化流模型稱作自回歸流模型（Autoregressive Flow）［19-20］，其中雙射函數(shù) fi具有如下的自回歸形式：

式中：τ稱作變換函數(shù)（transformer）［14］；ci稱作第i 個調(diào)節(jié)函數(shù)（conditioner）；。

該公式的自回歸特性體現(xiàn)在：第i 個調(diào)節(jié)函數(shù)ci只能以第1 個至第i-1 個維度的變量z1:i-1為輸入，其雅克比矩陣為下三角陣。將公式（2）作為Normalizing 方向的逆變換g［19］，并選擇仿射變換形式的變換函數(shù)τ［14］，如下：

式（3）中的變換函數(shù)τ是可逆的，且形式簡潔且引入的參數(shù)量較少，在實測裂隙數(shù)據(jù)的小樣本數(shù)據(jù)集上不易過擬合［19-20］，同時對于裂隙數(shù)據(jù)具備足夠的擬合能力。此時雅克比行列式的絕對值為：

用掩模自編碼器（Masked AutoEncoder for Density Estimation，MADE）［29］構(gòu)建調(diào)節(jié)函數(shù)ci，對多參數(shù)之間的相關(guān)性進行建模，逐個計算每個參數(shù)的條件概率并依據(jù)鏈?zhǔn)椒▌t求得聯(lián)合概率密度。在訓(xùn)練自回歸流模型的過程中，直接將負對數(shù)似然函數(shù)作為損失函數(shù)，通過梯度下降對每一個雙射函數(shù) fi內(nèi)部的仿射變換參數(shù)與掩模自編碼器參數(shù)進行優(yōu)化。

本文提出的DPCAF 模型中，基于高斯混合分布與密度峰值聚類改進自回歸流模型的方法如下：

首先，建立混合高斯分布作為DPCAF 模型的基礎(chǔ)分布。高斯混合分布是一種混合分布，記作：

式中：Gi為混合分支，服從高斯分布；為支密度；為混合比例。

利用DensityPeak 算法搜索密度峰值點作為聚類中心，將簇的數(shù)量作為高斯混合分布中分量的數(shù)量，定義DPCAF 模型的基礎(chǔ)分布，其中每一個分量對應(yīng)裂隙網(wǎng)絡(luò)的一個優(yōu)勢分組。將深度生成模型的先驗分布由單峰分布改進為混合分布是提高模型對于多峰分布擬合能力的有效手段［30-32］。以赤平極射投影圖上的裂隙產(chǎn)狀數(shù)據(jù)為例，對比自回歸流模型與DPCAF 模型的基礎(chǔ)分布，如圖3所示。

圖3 DPCAF 模型改進自回歸流模型示意

DensityPeak 算法所定義的密度峰值點包含兩方面特征［25］：（1）該點的局部密度大于其鄰近數(shù)據(jù)點；（2）該點與比它密度更高的數(shù)據(jù)點之間的距離較遠。

式中：ρi為局部密度，1/mm3；dc為截斷距離，mm；dij為數(shù)據(jù)點i 至數(shù)據(jù)點j 的距離，mm。當(dāng)a<0時當(dāng)a ≥0 時。

式中δi為與密度更高點的距離，mm。

其次，在以局部密度ρi為橫坐標(biāo)，以距離δi為縱坐標(biāo)，繪制出DensityPeak 決策圖上，自動提取出聚類中心附近的少量樣本點，作為已知分組標(biāo)簽的訓(xùn)練數(shù)據(jù)。搜索標(biāo)簽數(shù)據(jù)的規(guī)則如下：（1）樣本點按照其在決策圖上的局部密度從大到小的順序篩選；（2）樣本點與具有更高密度的最近點的距離大于設(shè)定的δ0距離閾值，mm。

最后，利用分布擬合空間的實測數(shù)據(jù)與DensityPeak 算法提取的部分?jǐn)?shù)據(jù)標(biāo)簽訓(xùn)練DPCAF 模型。將自回歸流的基礎(chǔ)分布定義為：

式中：高斯混合分布的混合分支數(shù)量為M，其中第k 個是均值為μk標(biāo)準(zhǔn)差為σk的高斯分布。

在最大似然估計的過程中對式（9）中的聯(lián)合似然函數(shù)進行最大化。根據(jù)貝葉斯定理，利用上述概率密度估計的結(jié)果計算出已知樣本點x 條件下標(biāo)簽l 的條件分布，并訓(xùn)練貝葉斯分類器實現(xiàn)裂隙分組。

3.2 基于DPCAF 模型的壩基巖體裂隙網(wǎng)絡(luò)多參數(shù)模擬流程為了更好地揭示鉆孔圖像數(shù)據(jù)中傾角、傾向、開度三個裂隙幾何參數(shù)的分布，本文建立了一個可逆的球坐標(biāo)系轉(zhuǎn)換，將實測數(shù)據(jù)點映射到分布擬合空間的點，如圖4所示。其中：α 為裂隙面傾向，（°）；β 為裂隙面傾角，（°）；δ 為裂隙開度，又稱作隙寬，mm；x，y，z 為分布擬合空間中的直角坐標(biāo)系下的坐標(biāo)。其中對于裂隙開度這一維度，本文使用數(shù)據(jù)點到球坐標(biāo)中心點的徑向距離對其進行表征，建立球面半徑與開度的轉(zhuǎn)換關(guān)系為，單位統(tǒng)一采用mm，其中δ0為開度的平均水平，實際工程的裂隙開度一般位于毫米至厘米量級，因此取δ0為10 mm，R0是開度等于δ0時的球面半徑，其數(shù)值也取為10 mm，λ 為縮放系數(shù)，本文取為0.1。

圖4 裂隙傾向、傾角、開度參數(shù)的球坐標(biāo)系轉(zhuǎn)換過程

DFN 多參數(shù)模擬的具體步驟如下：

（1）將傾角、傾向、開度的原始實測數(shù)據(jù)通過上述球坐標(biāo)系轉(zhuǎn)換為分布擬合空間的數(shù)據(jù)，并對DPCAF 模型進行訓(xùn)練，將DPCAF 模型輸出的分布擬合空間每個數(shù)據(jù)點的標(biāo)簽映射回原始數(shù)據(jù)，得到三個參數(shù)的優(yōu)勢分組結(jié)果；

（2）DPCAF 模型輸出分布擬合空間的聯(lián)合概率密度函數(shù)，對從中采集的樣本進行球坐標(biāo)系轉(zhuǎn)換的逆變換，從而得到傾角、傾向、開度參數(shù)的隨機模擬結(jié)果；

（3）基于裂隙網(wǎng)絡(luò)參數(shù)化建模計算模擬跡長分布下實測跡長數(shù)據(jù)平均對數(shù)似然，以裂隙直徑Gam?ma 分布的形狀參數(shù)和逆尺度參數(shù)為控制變量，以最大化似然函數(shù)為優(yōu)化目標(biāo)，利用粒子群優(yōu)化（Par?ticle Swarm Optimization，PSO）算法［33］進行求解，從而得到直徑參數(shù)的隨機模擬結(jié)果。

4 工程實例

以我國西南某水電工程的壩基巖體為研究對象，基于所提出的基于DPCAF 模型的DFN 多參數(shù)模擬方法，針對河床高程范圍內(nèi)的建基面下方50m×30m×30m（長×寬×深）區(qū)域范圍內(nèi)的裂隙網(wǎng)絡(luò)進行三維建模。鉆孔圖像數(shù)據(jù)包含309 條實測裂隙，開挖面測窗數(shù)據(jù)包含231 條實測裂隙，如圖5所示。

圖5 壩基巖體裂隙數(shù)據(jù)

4.1 基于DPCAF 模型的壩基DFN 多參數(shù)模擬從鉆孔圖像中提取的傾向、傾角與開度實測數(shù)據(jù)通過球坐標(biāo)系轉(zhuǎn)換至分布擬合空間的結(jié)果如圖6所示，圖中黑色的點代表轉(zhuǎn)換后得到的數(shù)據(jù)，半透明的半球面是開度等于10 mm 的等參面。

首先，利用DensityPeak 聚類算法確定聚類的分組數(shù)量，以及每一組的聚類中心。可以明顯地發(fā)現(xiàn)DensityPeak 決策圖的右上角出現(xiàn)三個聚類中心點，而其余點均分布于靠近橫軸的區(qū)域，因此優(yōu)勢分組數(shù)量為3 組。將決策圖中的點映射回三維坐標(biāo)系內(nèi)可以定位三個聚類中心的位置，如圖6所示。在DensityPeak 決策圖上自動搜索聚類中心附近的標(biāo)簽數(shù)據(jù)，設(shè)定標(biāo)簽數(shù)據(jù)在決策圖中縱坐標(biāo)的閾值為圖6 中的紅色虛線，最終得到的每一組的標(biāo)簽數(shù)據(jù)為20 個。

圖6 DensityPeak 算法搜索聚類中心

其次，利用所有實測數(shù)據(jù)點以及DensityPeak算法自動搜索得到的標(biāo)簽信息對DPCAF 模型進行訓(xùn)練。圖7 對訓(xùn)練完成的DPCAF 模型從分布擬合空間的數(shù)據(jù)中學(xué)習(xí)的三維聯(lián)合概率密度的形狀進行了可視化，其中對于每組裂隙繪制了兩個半透明的等概率密度曲面。

圖7 三維聯(lián)合分布概率密度估計結(jié)果

再者，結(jié)合鉆孔圖像與開挖面測窗中的裂隙數(shù)據(jù)，通過體密度估計得到三組裂隙網(wǎng)絡(luò)中裂隙數(shù)量分別為530、460 與470，并采用泊松過程對圓盤裂隙的中心點位置進行隨機模擬。經(jīng)過擬合優(yōu)度檢驗，三組裂隙的實測跡長都服從對數(shù)正態(tài)分布。假設(shè)圓盤裂隙直徑服從Gamma 分布，通過動態(tài)優(yōu)化得到第一組的形狀參數(shù)為3.804，逆尺度參數(shù)為0.827，第二組的形狀參數(shù)為3.075，逆尺度參數(shù)為1.235，第三組的形狀參數(shù)為4.982，逆尺度參數(shù)為1.400。

根據(jù)所建立的概率分布可以生成多個三維離散裂隙網(wǎng)絡(luò)，對其中一次隨機模擬生成的范例進行可視化如圖8所示。進一步對該裂隙網(wǎng)絡(luò)的跡線模擬結(jié)果進行數(shù)值檢驗，采用單尾F 檢驗與雙尾t 檢驗，顯著度設(shè)置為0.05，如表1所示。F 檢驗與t 檢驗的P 值均大于顯著度，說明模擬跡長的方差和均值都與實測結(jié)果無顯著差異。提出的DFN 建模方法基于深度生成模型充分挖掘了壩基裂隙巖體內(nèi)部鉆孔與開挖面測窗數(shù)據(jù)所蘊含的裂隙幾何信息，多元數(shù)據(jù)互相印證，提高了DFN 模型的可靠性，為水電工程壩基裂隙巖體的工程地質(zhì)分析與決策提供了重要支撐。

圖8 壩基三維離散裂隙網(wǎng)絡(luò)隨機模擬范例

表1 跡線隨機模擬范例數(shù)值檢驗

4.2 對比分析與討論首先，通過DPCAF 多參數(shù)模擬與傳統(tǒng)單參數(shù)模擬結(jié)果的對比分析，以及基于典型分布假設(shè)方法與DPCAF 方法的傾向與傾角二參數(shù)聯(lián)合分布估計的對比分析，驗證DPCAF 模型的準(zhǔn)確性；其次，為了驗證DPCAF 模型對于多維參數(shù)聯(lián)合分布估計的有效性，對基于GMM 模型與DP?CAF 模型的傾向、傾角、開度多參數(shù)模擬結(jié)果進行對比分析與討論。

4.2.1 裂隙跡線數(shù)據(jù)參數(shù)相關(guān)性的單參數(shù)模擬與DPCAF 模擬對比分析 本案例采用Mendoza-Torres等［10］的文章中給出的裂隙跡線數(shù)據(jù)進行對比分析，如圖9（a）所示，該裂隙網(wǎng)絡(luò)為一組由張拉裂隙與剪切裂隙構(gòu)成的共軛裂隙系統(tǒng)，包含兩組裂隙，其中一組優(yōu)勢裂隙沿南北方向延伸，具備較大的跡長但裂隙數(shù)量較少，另一組裂隙沿東西方向延伸，跡長相對較小但裂隙數(shù)量較多，因此，該裂隙系統(tǒng)的走向與跡長具備特定的依賴關(guān)系［10］，如圖9（a）所示，其中跡長為前10%的裂隙跡線標(biāo)記為黑色，其余裂隙跡線標(biāo)記為藍色。跡長與走向分別進行單參數(shù)隨機模擬生成的一個范例如圖9（b）所示，采用DPCAF 模型對跡長與走向進行多參數(shù)隨機模擬生成的一個范例如圖9（c）所示。顯然，實測數(shù)據(jù)中跡長較大的裂隙主要集中于南北走向，然而，采用傳統(tǒng)的單參數(shù)方法模擬的大跡長裂隙沒有呈現(xiàn)出明顯的優(yōu)勢走向，甚至有更多的大跡長裂隙為東西走向。因此，單參數(shù)模擬方法丟失了裂隙跡長與走向之間的依賴關(guān)系。相反，圖9（c）中的模擬結(jié)果與圖9（a）具備良好的一致性，因此DPCAF 模型在進行裂隙幾何參數(shù)聯(lián)合分布估計的過程中能夠良好地捕捉跡長與走向之間的依賴關(guān)系。

圖9 2D 裂隙網(wǎng)絡(luò)隨機模擬范例對比

為了將走向與跡長之間的關(guān)聯(lián)結(jié)構(gòu)進行可視化，將圖9 中裂隙的幾何參數(shù)繪制于圖10（a）—（c）的散點圖中，并進一步利用Mendoza-Torres 等［10］的研究中所采用的方法，將散點圖轉(zhuǎn)換為分位數(shù)分布圖，如圖10（d）—（f）所示。分位數(shù)分布圖展示的是分位數(shù)的概率，例如，對于散點圖中橫坐標(biāo)所表示的裂隙走向處于中位數(shù)的點在分位數(shù)分布圖中的橫坐標(biāo)u=0.5，而散點圖中橫坐標(biāo)處于四分位數(shù)的點在分位數(shù)分布圖中的橫坐標(biāo)u=0.25，對于縱坐標(biāo)的轉(zhuǎn)換方法同理。實測數(shù)據(jù)中兩個變量的秩相關(guān)系數(shù)為0.635，說明走向數(shù)據(jù)與跡長數(shù)據(jù)不是獨立的［10］。分位數(shù)分布圖10（d）與圖10（f）中的點都呈現(xiàn)出相同的分布規(guī)律，然而圖10（e）中的點趨向于離散的均勻分布，因此DPCAF 多參數(shù)模擬的結(jié)果有效地還原了走向與跡長兩個變量之間的二元結(jié)構(gòu)，能夠克服傳統(tǒng)單參數(shù)方法難以模擬裂隙參數(shù)之間關(guān)聯(lián)關(guān)系的不足。

圖10 裂隙走向與跡長隨機模擬范例的散點圖與分位數(shù)分布圖對比

4.2.2 裂隙產(chǎn)狀數(shù)據(jù)二參數(shù)聯(lián)合分布的典型分布模擬與DPCAF 模擬對比分析 目前離散裂隙網(wǎng)絡(luò)模擬方法中，涉及到聯(lián)合分布估計的參數(shù)只有傾向和傾角參數(shù)，通常假設(shè)傾向與傾角服從雙正態(tài)分布或Fisher 分布。Shanley-Mahtab 數(shù)據(jù)集［34］為Shanley 與Mahtab 從美國亞利桑那州San Manual 銅礦收集的286 個裂隙產(chǎn)狀測量值。該數(shù)據(jù)集被廣泛地應(yīng)用于關(guān)于產(chǎn)狀分布的研究中［11］。因此，本文將所提出的DPCAF 模型應(yīng)用于Shanley-Mahtab 數(shù)據(jù)集中第一組與第二組產(chǎn)狀數(shù)據(jù)的分布估計，并與傳統(tǒng)方法進行對比，結(jié)果如表2所示。Fisher 分布的估計精度大于雙正態(tài)分布，但兩者的對數(shù)似然函數(shù)都小于0，而DPCAF 模型對兩組產(chǎn)狀數(shù)據(jù)分布估計的對數(shù)似然函數(shù)都在8 以上。可見，實測產(chǎn)狀數(shù)據(jù)的分布可能并不完全服從雙正態(tài)分布與Fisher 分布的假設(shè)，而DPCAF 模型的概率密度估計結(jié)果更加接近數(shù)據(jù)的實際分布。

表2 產(chǎn)狀數(shù)據(jù)概率密度估計精度對比

圖12 采樣模擬裂隙數(shù)據(jù)范例可視化對比

4.2.3 鉆孔圖像數(shù)據(jù)多參數(shù)聯(lián)合分布的GMM 模擬與DPCAF 模擬對比分析 為了對傾向、傾角與開度數(shù)據(jù)的概率密度估計的精度進行評價，將提出的DPCAF 模型與目前主流的有限混合模型（Finite Mixture Model）——高斯混合模型（Gaussian Mixture Model，GMM）進行對比，如圖11所示。

三維的GMM 模型的每一個分支對應(yīng)的等概率密度曲面都為“置信橢球”，但由于實測裂隙開度參數(shù)的分布較為集中，對應(yīng)于球坐標(biāo)系下的三維空間點基本分布于開度等參球面附近，而不是呈橢球狀分布。GMM 用“置信橢球”擬合實際數(shù)據(jù)時表現(xiàn)出了一些系統(tǒng)性的偏倚，例如，圖11（b）中綠色的第二組和圖11（c）中紫色的第三組裂隙所估計的置信橢球的一端明顯向開度等參球面外側(cè)翹起。相反，DPCAF 模型可以建立任意復(fù)雜不規(guī)則形狀的分布，所估計的概率密度形狀良好地擬合了真實數(shù)據(jù)點的分布。采用評價概率密度估計的精度的常用評價指標(biāo)——平均對數(shù)似然（Mean Log-Likelihood）［19］進行量化對比。其中，DPCAF 模型的平均對數(shù)似然都在21 至22左右，而GMM 模型的平均對數(shù)似然都在-6 以下，可見DPCAF 的概率密度估計精度高于GMM。

進一步對比分析GMM 與DPCAF 模型的采樣模擬的誤差，控制采集樣本的數(shù)量相同，采用模擬數(shù)據(jù)分布與實測數(shù)據(jù)分布之間的Wasserstein 距離［35］作為模擬誤差的評價指標(biāo)，如表3所示。

表3 采樣模擬裂隙數(shù)據(jù)范例誤差對比

圖12 以第三組裂隙為例，分別采用球坐標(biāo)可視化與DFN 可視化方法對三維數(shù)據(jù)分布建立形象直觀的表達。在DFN 可視化方法中，固定每個圓盤結(jié)構(gòu)面的中心點位置和半徑，將采樣數(shù)據(jù)所包含的傾向、傾角和開度幾何參數(shù)分別賦值給不同的圓盤模型，得到虛擬的離散網(wǎng)絡(luò)模型，其中為了提升可視化效果將開度數(shù)據(jù)乘以500 作為圓盤的厚度，如圖12所示?？梢园l(fā)現(xiàn)，DPCAF 采樣的每一組裂隙的分布以及三維裂隙的整體分布都比GMM 的采樣結(jié)果更加接近實測數(shù)據(jù)，從而證明了本文提出的DPCAF 模型的模擬精度更高。

5 結(jié)論

裂隙幾何參數(shù)分布的估計與擬合是三維離散裂隙網(wǎng)絡(luò)隨機模擬的核心環(huán)節(jié)。然而，現(xiàn)有的DFN建模方法基于Gamma 分布、對數(shù)正態(tài)分布、Fisher 分布等典型分布對單參數(shù)的邊緣分布或傾向傾角的聯(lián)合分布進行估計，面臨著無法拓展到更多維裂隙參數(shù)的復(fù)雜聯(lián)合分布的難題。針對上述問題，本文提出了基于改進自回歸流模型的DFN 多參數(shù)模擬方法，取得了如下成果：

（1）考慮裂隙幾何參數(shù)具有優(yōu)勢分組的特點改進自回歸流模型，將標(biāo)準(zhǔn)化特征空間的高斯分布改進為高斯混合分布，并結(jié)合DensityPeak 算法搜索聚類中心，提出了DPCAF（Density Peak Clustering Autoregressive Flow）模型，克服了自回歸流模型對于裂隙數(shù)據(jù)多峰分布擬合能力的不足；

（2）將自回歸流模型中基于可逆變換擬合復(fù)雜多維聯(lián)合分布的思想引入DFN 建模中，有效解決了DFN 多參數(shù)模擬所面臨的多維聯(lián)合分布的概率密度估計與采樣的核心難題，在DFN 建模的過程中同時完成多參數(shù)優(yōu)勢分組與多參數(shù)概率密度估計；

（3）工程應(yīng)用表明，對于實測裂隙數(shù)據(jù)的復(fù)雜聯(lián)合分布，本文提出的DPCAF 模型相比與傳統(tǒng)基于典型分布的估計方法具備更強的擬合能力，并且能夠有效還原多參數(shù)之間的關(guān)聯(lián)關(guān)系；DPCAF 模型相比于雙正態(tài)分布與Fisher 分布對于產(chǎn)狀數(shù)據(jù)的概率密度估計精度更高，相比于GMM 模型對于鉆孔圖像中的傾向、傾角、開度數(shù)據(jù)聯(lián)合分布的概率密度估計的平均精度更高，且采樣結(jié)果的誤差縮小了7.84%。因此，DPCAF 模型能夠有效提高多參數(shù)擬合的準(zhǔn)確性，保證DFN 建模的精度。

本文提出的改進自回歸流模型能夠為工程分析過程中涉及到的多維聯(lián)合概率密度估計問題提供一種新的技術(shù)手段［36-37］。此外，隨著水電工程地質(zhì)勘測與感知技術(shù)的不斷進步［38］，未來將可能獲取更多種類的裂隙數(shù)據(jù)，可以將其引入DFN 多參數(shù)模擬中進一步提升裂隙網(wǎng)絡(luò)模型的可靠性。