基于貝葉斯機器學習的大渡河上游巖石物理力學參數概率分布研究

鄧 欽 宣, 張 韶 鵬, 邵 偉, 楊 洋

(1.四川足木足河流域水電開發有限公司,四川 成都 610041;2.四川大學水利水電學院,四川 成都 610200)

0 引 言

在水利水電工程中,大壩壩基及邊坡穩定是保證水利水電工程安全至關重要的影響因素之一[1]。針對水利水電工程結構,《水利水電工程結構可靠度設計統一標準》[2]要求采用基于可靠度理論的分項系數設計方法進行結構設計,而可靠度設計要求提供巖石物理力學參數的分布信息。巖石在形成過程中經歷了長時間跨度的復雜地質過程,巖體的材料和工程參數具有不同的分布特征且通常具有相關性,其聯合分布具有復雜特征[3-4]。若各參數的聯合概率分布估計不當,或者簡單地忽略它們的相關性,會導致錯估結構失效概率,導致有偏差的設計結果,從而直接影響工程安全[5-6]。

以準確估計巖石物理力學概率分布為目標,國內外諸多學者進行了大量研究,經歷了從確定性分析到不確定性分析,從簡單概率統計到概率分布擬合,從單參數邊緣分布到多參數聯合分布的研究歷程,發展到了考慮不確定性的各參數不同邊緣分布類型及不同相關結構的概率分布構建階段[7-12]。然而,因巖石土體材料的場地或區域特異性、材料天然變異性和形成過程的復雜性,其聯合分布常有悖于經典概率分布,且需要事先進行分布類型和相關結構的假設,而現有方法無法恰當解決該難題[13]。為此,本文將基于貝葉斯機器學習框架的高斯混合模型(Gaussian Mixture Model, GMM)應用于大渡河上游巖石物理力學參數概率分布構建任務中,在考慮統計不確定條件下精準刻畫聯合分布特征及各參數間的相關性,構建區域性多元巖石物理力學參數概率分布,為該區域的后續工程建設和設計提供參考。

1 高斯混合模型

高斯混合模型是混合分量為高斯分布的一種特殊混合模型,于1894年由生物統計學家Karl Person首次提出并應用于生物統計學中的偏態數據分析[14-15],并在之后的百年間獲得了持續而深入的發展,已成為機器學習領域不可或缺的模型之一。高斯混合模型不僅具有極強的靈活性,而且繼承了高斯分布的數學便利性,已被成功應用于多學科、多領域的研究,包括航空航天、醫學、信號處理與分析、經濟學以及社會科學等諸多方面[16]。

p(X|ω)=π1N(X|μ1,Σ1)+…+πkN(X|μk,Σk)+πKN(X|μK,ΣK)

(1)

一般形式為:

(2)

式中K為高斯分量的個數;N(X|μk,Σk)為1到K個高斯分量中第k個分量的概率密度函數;πk為第k個高斯分量的權重參數;向量π=[π1,…,πK]也被稱為權重分布。為使GMM模型的概率密度函數在其定義域上積分為1,權重πk必須滿足以下約束條件:

(3)

式中μk為第k個高斯分量的均值向量;Σk為第k個高斯分量的協方差矩陣,由各維度變量的方差和各變量之間的相關系數矩陣表示如下:

(4)

而在研究模型復雜度時,協方差矩陣Σk作為d(d+1)/2個參數進行計算。因此,多維高斯混合模型共有K(d+d(d+1)/2+1)-1個模型參數,模型參數求解涉及高維求解問題。為直觀展示其參數規模,圖1展示了不同變量維度條件下模型參數個數Np與高斯分量個數K的對應關系。

圖1 不同變量維度條件下模型參數個數Np與高斯分量個數K的對應關系

由圖1可知,高斯混合模型的模型參數個數Np隨著高斯分量個數K的增加而線性增加,10維高斯混合分布在高斯分量個數為3時,即有高達198個參數。高斯模型復雜度在為概率分布擬合帶來極大便利的同時,也給模型學習帶來了挑戰,因此,本文采用貝葉斯機器學習方法解決模型學習難題。

2 貝葉斯機器學習方法

除去上文提到的模型復雜度難題外,高斯混合模型的學習還因其為隱變量模型的特性存在標簽切換問題,而在考慮不確定性分析的前提下,被廣泛應用于求解GMM模型的EM算法失效[17]。因此,本節應用筆者之前提出的貝葉斯學習框架解決模型學習問題。

2.1 貝葉斯參數學習

(5)

(6)

(7)

式中p(xi|Si,ωK,MK)為觀測數據xi屬于第Si個高斯分量的似然函數;p(Si|ωK,Mg)為在給定模型參數ωK和模型MK的情況下,觀測數據xi屬于第Si個高斯分量的概率。一旦隱變量SK被確定,則觀測數據xi所歸屬的高斯分量即被確定。為定量描述觀測數據xi的歸屬情況,引入變量nk,k= 1, 2, …,K,其代表了屬于第k個高斯分量的觀測數據樣本的個數。

(8)

式中xi,k代表了屬于第k個高斯分量的第i個觀測數據;p(xi,k|μk,∑k為第i個觀測數據歸屬于第k個高斯分量概率密度函數(PDF)值,可表示為:

(9)

GMM參數ωK(例如:π,μk, Σk,k= 1, 2, …,K)的共軛先驗分布p(ωk|Mk)的設置如下[16-19]:

π～Dirichlet(α1,…,αk)

(10)

μk～N(bk,∑k/Bk)

(11)

(12)

式中 權重參數π= [π1,π2, …,πK] 服從狄利克雷分布(Dirichlet distribution),其概率密度函數PDF為:

(13)

式中a1,a2, …,aK為迪利克雷分布的分布參數。

均值參數μk的共軛先驗分布為以bk為均值,∑k/Bk為方差的多維高斯分布,其中Bk為事先給定的比例因子。

(14)

式中 Гd(·)為d維的Gamma函數,其表達式為:

(15)

在多維高斯分布的求解算法中,Wishart分布常被用來構造精度矩陣的共軛先驗分布。概率密度函數表達式中Ck為Wishart分布的尺度參數,其為一個d×d的對稱非奇異方陣。ck為Wishart分布的自由度參數,其與分布的自由度vk存在線性關系(即ck=vk/2)。

至此,表征GMM參數和隱變量的聯合后驗分布公式(5)的分子部分已經完全給出。但是,由于歸一化常數仍然未知,仍舊無法得到參數和隱變量的后驗分布樣本。為了在考慮標簽切換的前提下解決此問題,本文采用隨機吉布斯抽樣方法進行模型參數求解,此處不贅述。

2.2 貝葉斯模型比選

在應用數學模型進行數據分析建模時,模型不確定性是無法避免的,而對于高斯混合模型而言,其模型不確定性的最大來源就是高斯分量個數K的不確定性。確定高斯混合模型中的高斯分量個數會不可避免地涉及過擬合混合模型的學習問題,而由于標簽切換和過擬合問題的存在,此時的高斯混合模型經常是不可識別的[18]。

隨著計算科學的發展,科研工作者們提出了許多不同的模型選擇方法來嘗試考慮模型選擇的不確定性。其中基于模型證據的貝葉斯模型比選方法同信息準則等方法比較具有優越性,因此,本文應用貝葉斯模型比選方法,權衡模型復雜度和擬合優度,選擇最優備選模型。

(16)

(17)

然而,當高斯分量個數大于1時,高斯混合模型的模型證據并沒有解析解。在算法實踐中,模型證據需要通過合適的后處理手段(postprocessing manner)求得。在后處理算法中,基于隨機模擬方法的模型證據求解方法(Simulation-based approximations)的應用最為成功。本研究主要應用隨機模擬方法中的橋采樣(Bridge Sampling, BS)方法求解GMM的模型證據。

3 大渡河上游巖石物理力學參數概率分布構建

3.1 區域、數據概況

大渡河位于四川西部,是岷江水系最大的支流,年徑流量470億m3,干流全長1 062 km,天然落差4 175.0 m,為我國重要的水電能源基地之一。如今,大渡河干流形成了以下爾呷為龍頭的28級開發方案,其中上游共規劃有3級水電站,自上而下分別為下爾呷、巴拉和達維水電站,其中巴拉水電站在建,下爾呷、達維正處于項目前期階段。

大渡河上游巖石主要以變質砂巖和板巖為主,巴拉水電站區域分布有花崗巖侵入區,本文收集了大渡河上游3級水電站的巖石物理力學參數試驗數據,以巖體干密度ρd、飽和吸水率w、飽和抗壓強度Rw為主要研究對象,考慮統計不確定性和參數相關性,通過所提方法構建大渡河上游巖石物理力學參數概率分布,驗證所提方法有效性,為大渡河上游區域水電工程可靠度設計提供依據。大渡河上游巖石物理力學參數基本統計信息見表1。

表1 大渡河上游巖石物理力學參數基本統計信息

如表1所示,所選的三個參數的統計特征之間存在較大差異,其中干密度ρd的分布較集中,變異系數僅為0.01,而飽和吸水率w及飽和抗壓強度Rw的變異性較大,其中飽和吸水率的變異系數甚至高達0.91,這無疑為參數概率分布特征的準確表征帶來了困難。

大渡河上游巖石物理力學參數試驗數據的二維散點矩陣(圖2),展示了各參數的頻率統計直方圖以及各參數之間的二維數據分布散點。

圖2 大渡河上游巖石物理力學參數試驗數據的二維散點圖矩陣(n=222)

如圖2所示,參數并不嚴格服從正態分布且具有多模態特征,且各參數之間明顯具有非線性相關關系,很難通過現有的方法對這種復雜相關的特征進行描述。為了恰當表述巖石力學參數的多模態特征,同時處理參數之間的非線性相關關系,接下來應用所提方法對大渡河上游巖石物理力學參數試驗數據進行分析。

3.2 GMM模型選擇結果

首先假設共有5個備選模型M1,M2,M3,M4,M5,即高斯分量個數K的取值范圍為從1到5的正整數(Kmax= 5),模型下標的數字即代表了其具有的高斯分量個數。備選模型的權重參數,各高斯分量的均值和精度矩陣的共軛先驗分布,分別取為狄利克雷分布、正態分布和Wishart分布。然后應用RGS-GS方法學習高斯混合模型,設置RGS方法后驗樣本數、burn-in樣本數和BS方法中的重要性抽樣樣本數均為10 000。

貝葉斯模型選擇結果見圖3。圖中展示了5個備選模型的模型概率和模型證據的對數值,其中模型證據對數值用帶有實心方形標記的實線表示,模型概率用直方圖表示。由圖3可知,隨著模型分量個數的增加,模型證據的對數值在K= 4時達到最大值-489.08,之后下降至-489.68。因此,最優模型高斯分量個數為4的高斯混合模型GMM4,其對應的模型概率為0.59。

圖3 貝葉斯模型選擇結果

3.3 GMM建模結果

在進行最優模型結構進行選擇后,應用所提貝葉斯學習框架對模型參數進行學習。圖4(a)到(c)分別展示了各巖石力學參數的邊緣概率分布擬合結果。

圖4 各巖石力學參數的邊緣概率分布擬合結果

圖4中,黑色實線代表了PDF的最可能值,由黑色虛線圍成的區域代表了PDF的95%置信區間,反映了模型統計不確定性的大小,直方圖為歸一化的頻率分布直方圖,展示了數據本身的分布特征及模型統計不確定性的大小。由圖4可知,歸一化的頻率分布直方圖分布于95%置信區間內,且與MPV值十分接近,這表明高斯混合模型GMM4的邊緣概率分布不僅能很好地擬合觀測數據的實際邊緣分布特征,還能正確地表征其不確定性。

通過二維變量聯合概率密度對數值等值線圖,展示多維高斯混合模型對多維巖土體參數的聯合分布的表征能力和對參數間相關性及相關結構的刻畫能力,各巖石力學參數的聯合概率分布擬合結果見圖5。

圖5 各巖石力學參數的邊緣概率分布擬合結果

圖5展示了各變量間二維聯合概率密度函數對數值的等值線圖,圖中用黑色空心圓形標記代表實測巖石物理力學參數數據;黑色實心正方形標記代表了識別出的各高斯分量的均值,也可認為是各高斯分量的中心;黑色虛線代表了各參數的聯合概率密度對數值等值線,為增強結果的可視性,填充顏色從淺到深的演變對應概率密度函數對數值的從大到小變化。由各子圖可知,觀測樣本點集中分布在顏色較淺的高概率密度區域,隨著填充顏色從淺到深,觀測樣本點的分布密度也逐漸下降。此現象說明學習出的高斯混合模型能夠有效刻畫參數空間中的概率密度變化情況。各子圖中的概率密度函數等值線并不是規則的橢圓形,而是呈現出隨著數據點密度變化的不規則圖形。這說明GMM模型表征的相關結構并不是高斯型或者其他傳統類型的相關結構,而是由數據的特征決定的“數據驅動”相關結構,可以有效表述巖石物理力學參數分布特征。

GMM4的參數的最可能值如表2所示。

表2 GMM4參數最可能值

4 結 語

本文將基于貝葉斯機器學習框架的高斯混合模型應用于大渡河上游巖石物理力學參數概率分布構建任務中,在考慮統計不確定條件下精準刻畫聯合分布特征及各參數間相關性,構建了區域性多元巖石物理力學參數概率分布,得出以下結論:

(1)大渡河上游巖體干密度ρd、飽和吸水率w及飽和抗壓強度Rw具有較大變異性,且分布具有多峰、多模態特征,具有復雜相關結構。

(2)所提方法和模型打破了現有方法必須事先假設各參數概率分布類型及相關結構類型的假設,能夠在考慮統計不確定條件下精準刻畫聯合分布特征及各參數間相關性,有效表述了大渡河上游巖石物理力學參數分布特征。

(3)本文給出了所選參數聯合概率分布的GMM模型,明確了模型參數,可直接用于后續工程設計工作,為后續工程提供了參考。