周 洋

1 中國地震局地震研究所，武漢市洪山側(cè)路40號，430071 2 中國地震局地震大地測量重點(diǎn)實(shí)驗(yàn)室，武漢市洪山側(cè)路40號，430071 3 湖北省地震局，武漢市洪山側(cè)路48號，430071

地震預(yù)測的難度在于地球內(nèi)部的不可入性、地震的非頻發(fā)性及地震發(fā)生機(jī)理的復(fù)雜性。科學(xué)合理地進(jìn)行地震數(shù)據(jù)預(yù)測，對于救災(zāi)物資的有效配置及國計(jì)民生具有重要的指導(dǎo)意義[1-3]。隨著計(jì)算機(jī)硬件設(shè)備和智能算法的發(fā)展，利用計(jì)算機(jī)的高速計(jì)算能力和智能學(xué)習(xí)算法對地震活動進(jìn)行建模預(yù)測成為地震預(yù)測領(lǐng)域的主流研究方向?，F(xiàn)今的地震預(yù)測領(lǐng)域已成為橫跨多門基礎(chǔ)學(xué)科的綜合研究領(lǐng)域[4]，有學(xué)者利用機(jī)器學(xué)習(xí)算法從儲層預(yù)測、地震災(zāi)害信息預(yù)測、地震死亡人數(shù)評估等方面進(jìn)行研究[5-8]，并采用統(tǒng)計(jì)學(xué)方法和機(jī)器學(xué)習(xí)方法對股票價格進(jìn)行建模預(yù)測[9-10]。

目前地震數(shù)據(jù)預(yù)測仍面臨諸多問題，如當(dāng)預(yù)測時間和地理區(qū)域范圍較大時，特征指標(biāo)數(shù)據(jù)不能很好地描述地震的發(fā)生狀況，預(yù)測性能較差；神經(jīng)網(wǎng)絡(luò)等需要訓(xùn)練的算法，需要大量的異常數(shù)據(jù)集進(jìn)行訓(xùn)練，否則容易出現(xiàn)過擬合問題。針對這些問題，本文提出一種基于核混合效應(yīng)的回歸模型，并應(yīng)用于湖北地震臺站地球物理數(shù)據(jù)的預(yù)測仿真，結(jié)果表明，相較于傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)機(jī)器學(xué)習(xí)算法，該模型預(yù)測效果較好，預(yù)測精度較高[11-12]。

1 核混合效應(yīng)回歸模型

1.1 最小二乘支持向量機(jī)與線性混合模型

內(nèi)核提供了一種有效的方法來制定基于內(nèi)部乘積的線性算法的非線性概括，基于內(nèi)核的學(xué)習(xí)已越來越多地用于各個學(xué)科的數(shù)據(jù)分析程序中。支持向量機(jī)(SVM)是一種受監(jiān)督的基于內(nèi)核的算法，可用于分類和回歸應(yīng)用程序。SVM找到約束優(yōu)化問題的解，該解以原始形式用非線性特征映射函數(shù)φ(x)表示，通過引入拉格朗日乘數(shù)和滿足Mercer條件的核函數(shù)K(xi,xj)=φ(xi)Tφ(xj)，可將原始形式轉(zhuǎn)換為通常更易于求解的對偶形式。

最近引入了最小二乘支持向量機(jī)(LS-SVM)作為經(jīng)典SVM分類器的一種變體，其中不等式約束由等式約束代替，這種重構(gòu)允許通過將相應(yīng)的凸二次規(guī)劃問題簡化為線性方程組來解決SVM問題。除了提供簡單的軟件實(shí)現(xiàn)和增加數(shù)值穩(wěn)定性外，這種方法還可將經(jīng)典的SVM擴(kuò)展到解決數(shù)據(jù)分析和模式識別中的更多問題。

使用內(nèi)核可以通過在特征空間中執(zhí)行嶺回歸來解決回歸問題。給定代表N點(diǎn)坐標(biāo)的輸入矩陣X∈RN×d、d維空間中的x1,x2,…,xN和對應(yīng)的輸出向量y∈RN×1，并假設(shè)每個觀測值都與對應(yīng)的輸入有關(guān)。

yi=wTφ(xi)+b+ei,i=1,2,…,N

(1)

式中，yi為輸出向量觀測值，φ(xi)為輸入矩陣的坐標(biāo)，wT、b、ei為約束條件參數(shù)。

f(x)=wTφ(x)+b

(2)

在式(1)中給出的N個約束條件下，使用懲罰參數(shù)γ。格朗日函數(shù)寫為：

(3)

式中，αi(i=1,…,N)為拉格朗日乘數(shù)，e=[e1,e2,…,eN]T，α=[α1,α2,…,αN]T。通過定位靜止點(diǎn)找到w、b、ei、αi的最佳值，即為式(3)的拉格朗日函數(shù)點(diǎn)。使用S=[φ(x1),φ(x2),…,φ(xN)]T，最優(yōu)條件可總結(jié)為：

(4)

式中，Ip為大小為p的恒等矩陣，0p×q為0的p×q維矩陣，1p和0p分別為1和0的p×1維向量，k為特征空間的維數(shù)。消除向量w和e得到：

(5)

(6)

線性混合模型的觀測值可以寫為：

y=Xβ+Zu+e

(7)

式中，y為觀測向量，β為固定效應(yīng)向量(其中包括一個常數(shù)項(xiàng))，u～N(0,G)為未觀察到的隨機(jī)效應(yīng)向量，e～N(0,R)為隨機(jī)誤差的向量，X和Z分別為固定效應(yīng)和隨機(jī)效應(yīng)設(shè)計(jì)矩陣。向量u和e不相關(guān)，u和y的聯(lián)合密度最大化：

p(u,y|β,G,R)∝

(8)

產(chǎn)生亨德森的混合模型方程式(MME)：

(9)

1.2 回歸方程

i=i:N，j=i:ni

(10)

式中，yij為與第i次地震的第j次記錄相對應(yīng)的目標(biāo)變量，xij為固定效應(yīng)協(xié)變量矢量，β為固定效應(yīng)矢量，sij為包含與目標(biāo)變量沒有特定關(guān)系的預(yù)測變量的矢量，zij為隨機(jī)效應(yīng)協(xié)變量向量，ui為隨機(jī)效應(yīng)向量，eij為與yij相關(guān)的隨機(jī)誤差項(xiàng)。向量β包含常數(shù)項(xiàng)作為其第1項(xiàng)，向量xij的第1項(xiàng)為1。假定向量ui和ei=[ei1,ei2,…,ein]T遵循正態(tài)分布，分別具有均值零和協(xié)方差矩陣Gi和Ri，向量ui和ei不相關(guān)，則式(10)可用矩陣形式寫為：

y=Xβ+Ωα+Zu+e

(11)

R-1(y-Xβ-Zu-Ωα)+uTG-1u+αTA-1α]

(12)

設(shè)定l(β,u,α)關(guān)于β、u和α的偏導(dǎo)數(shù)為零，則得出以下方程組：

(13)

在前面的推導(dǎo)中，假設(shè)協(xié)方差矩陣G、A和R是已知的，并且已選擇具有相關(guān)參數(shù)的適當(dāng)核函數(shù)。但在大多數(shù)實(shí)際情況下，協(xié)方差矩陣和核函數(shù)必須由用戶確定，可以使用多種方法來估計(jì)這3個協(xié)方差矩陣，例如最小化廣義交叉驗(yàn)證得分法或最大化殘差最大似然(REML)法。一旦確定了G、A、R和Ω，就可利用式(13)確定固定效應(yīng)、隨機(jī)效應(yīng)和拉格朗日乘數(shù)的最佳值，并將其代入式(11)。模型對(x(i),s(i),z(i))描述的第i個觀測點(diǎn)的估計(jì)為：

E[y(i)|x(i),s(i),z(i)]=

(14)

未觀測點(diǎn)(x*,s*)的預(yù)測方程為：

(15)

1.3 逐點(diǎn)置信度和預(yù)測間隔

由于對條件期望的核估計(jì)可能存在偏差，因此必須注意以正確置信區(qū)間為中心。假設(shè)觀測值可以寫成均值函數(shù)f(x)和標(biāo)準(zhǔn)偏差函數(shù)σ(x)：

yi=f(xi)+σ(xi)εi,i=1,2,…,n

(16)

其中，E[εi|xi]=0,var[εi|xi]=1，以x為條件的期望和方差為E[y|x]=f(x)、var[y|x]=σ2(x)。f(x)的100(1-a)%置信區(qū)間對應(yīng)于邊界qa：

(17)

其中，P(·)表示概率。

可以證明式(17)為最小二乘支持向量機(jī)器模型，式(6)為線性平滑器。也就是說，存在一個向量L(x)=[l1(x),l2(x),…,ln(x)]T。

(18)

(19)

(20)

(21)

可以使用估算器近似為：

(22)

(23)

式中，Σ為由訓(xùn)練點(diǎn)處的預(yù)測方差組成的對角矩陣。使用式(22)和式(23)得到條件均值的近似及偏差校正的100(1-a)%點(diǎn)狀置信區(qū)間為：

(24)

(25)

z(i)GZT)(INn-XTUXTV-1)

(26)

Ω*AΩV-1(INn-XTUXTV-1)

(27)

式中，Ω*=[K(s*,s1),K(s*,s2),…,K(s*,sNn]。

以上算法在MATLAB 9.0軟件上測試通過，算法流程見圖1。

圖1 算法流程Fig.1 Algorithm flowchart

2 仿真實(shí)驗(yàn)

為驗(yàn)證本文模型算法的可行性和準(zhǔn)確性，選用地震地球物理觀測數(shù)據(jù)在此算法基礎(chǔ)上開展仿真實(shí)驗(yàn)，并與傳統(tǒng)的經(jīng)典神經(jīng)網(wǎng)絡(luò)算法實(shí)例作比較。首先選取水溫分鐘觀測值訓(xùn)練該神經(jīng)網(wǎng)絡(luò)模型，用訓(xùn)練數(shù)據(jù)進(jìn)化200次，訓(xùn)練過程中神經(jīng)網(wǎng)絡(luò)模型預(yù)測誤差變化趨勢見圖2。

圖2 網(wǎng)絡(luò)進(jìn)化過程Fig.2 Network evolution

當(dāng)訓(xùn)練次數(shù)達(dá)到200次時，預(yù)測誤差小于2.37×10-3，達(dá)到訓(xùn)練要求的目標(biāo)誤差0.01。隨機(jī)選取2020-10-11湖北省荊門地震臺4測點(diǎn)的水溫儀實(shí)測時序數(shù)據(jù)作為測試數(shù)據(jù)，圖3為BP神經(jīng)網(wǎng)絡(luò)預(yù)測、RBF神經(jīng)網(wǎng)絡(luò)預(yù)測、模糊神經(jīng)網(wǎng)絡(luò)預(yù)測及線性神經(jīng)網(wǎng)絡(luò)預(yù)測結(jié)果對比，圖4為基于核混合效應(yīng)回歸模型預(yù)測。

圖3 幾種神經(jīng)網(wǎng)絡(luò)預(yù)測算法的預(yù)測結(jié)果(荊門)Fig.3 Prediction results of several neural network prediction algorithms(Jingmen)

圖4 基于核混合效應(yīng)回歸模型的預(yù)測結(jié)果(荊門)Fig.4 Prediction results of regression model prediction based on nuclear mixed effects(Jingmen)

由圖3和圖4可知，本文回歸模型預(yù)測效果較好，其預(yù)測曲線的擬合度明顯高于其他算法。BP神經(jīng)網(wǎng)絡(luò)預(yù)測、RBF神經(jīng)網(wǎng)絡(luò)預(yù)測、模糊神經(jīng)網(wǎng)絡(luò)預(yù)測、線性神經(jīng)網(wǎng)絡(luò)預(yù)測及基于核混合效應(yīng)回歸模型預(yù)測的相對誤差分別為0.92%、0.83%、0.78%、0.54%、0.05%，本文算法模型的相對誤差曲線見圖5，其他幾種算法的相對誤差曲線類似，不再贅述。

圖5 基于核混合效應(yīng)回歸模型預(yù)測的相對誤差(荊門)Fig.5 Predicting relative error based on the regression model of nuclear mixed effects(Jingmen)

同理，選取2020-09-24湖北省房縣三海村地震臺1測點(diǎn)水位儀實(shí)測時序數(shù)據(jù)，結(jié)合上述5種神經(jīng)網(wǎng)絡(luò)算法進(jìn)行仿真實(shí)驗(yàn)，結(jié)果見圖6～7。

圖6 幾種神經(jīng)網(wǎng)絡(luò)預(yù)測算法的預(yù)測結(jié)果(房縣)Fig.6 Prediction results of several neural network prediction algorithms(Fangxian)

圖7 基于核混合效應(yīng)回歸模型的預(yù)測結(jié)果(房縣)Fig.7 Prediction results of regression model prediction based on nuclear mixed effects(Fangxian)

經(jīng)計(jì)算，對于水位數(shù)據(jù)預(yù)測而言， BP神經(jīng)網(wǎng)絡(luò)預(yù)測、RBF神經(jīng)網(wǎng)絡(luò)預(yù)測、模糊神經(jīng)網(wǎng)絡(luò)預(yù)測、線性神經(jīng)網(wǎng)絡(luò)預(yù)測及基于核混合效應(yīng)回歸模型預(yù)測的相對誤差分別為2.26%、1.78%、1.59%、1.34%、0.48%，其中基于核混合效應(yīng)回歸模型預(yù)測的相對誤差曲線見圖8。對于其他地球物理觀測手段(形變、重力、電磁等)產(chǎn)出的觀測數(shù)據(jù)資料，分析方法相同。

圖8 基于核混合效應(yīng)回歸模型預(yù)測的相對誤差(房縣)Fig.8 Predicting relative error based on the regression model of nuclear mixed effects(Fangxian)

3 結(jié) 論

本文提出一種基于核混合效應(yīng)回歸模型的數(shù)據(jù)預(yù)測新方法，并將該方法用于地震地球物理觀測數(shù)據(jù)的預(yù)測。對不同神經(jīng)網(wǎng)絡(luò)算法進(jìn)行仿真實(shí)驗(yàn)對比，結(jié)果表明，本文算法能很好地預(yù)測未來曲線趨勢，且相對誤差較低。結(jié)論如下：

1)本文方法使用最小二乘支持向量機(jī)來擴(kuò)展線性混合模型方程，該擴(kuò)展為復(fù)雜關(guān)系的建模提供了靈活的工具，同時考慮到觀測值之間的相關(guān)性，提出的方法在計(jì)算上是有效的，其將優(yōu)化問題簡化為線性方程組。

2)當(dāng)前常用的參數(shù)最小二乘法是普通最小二乘模型的擴(kuò)展，該模型旨在用于線性模型。對于線性模型，已知最小二乘估計(jì)量是無偏的，并且在所有線性無偏估計(jì)量中方差最低，但像LS-SVM估計(jì)器這樣的有偏估計(jì)器可能會比最小二乘估計(jì)器獲得更低的均方誤差。

3)通過仿真研究，假設(shè)線性數(shù)據(jù)生成函數(shù)，并將本文方法與參數(shù)最小二乘法的精度進(jìn)行比較，仿真結(jié)果表明，使用本文模型可以更好地處理涉及相關(guān)預(yù)測變量的問題。關(guān)于LS-SVM的許多結(jié)論都適用于非線性模型，這是因?yàn)長S-SVM使用非線性核處理非線性問題，并且一旦選擇了核，LS-SVM模型就相當(dāng)于分配給每個點(diǎn)的權(quán)重保持線性。

4)本文模型算法還可應(yīng)用于地震地球物理其他測項(xiàng)(地磁、形變、重力等)的數(shù)據(jù)預(yù)測，應(yīng)用前景廣泛，同時也可為較復(fù)雜的深度學(xué)習(xí)類算法框架模型的構(gòu)建提供實(shí)踐基礎(chǔ)。