本發(fā)明涉及一種地震前兆觀測數(shù)據(jù)異常檢測方法，具體涉及一種基于主成分分析的鉆孔應(yīng)變數(shù)據(jù)異常提取方法。

背景技術(shù)：

地球科學(xué)是一門以觀測和測量資料為基礎(chǔ)的學(xué)科，對觀測系統(tǒng)的不斷完善是發(fā)展地球科學(xué)、減輕災(zāi)害的必由之路。地球表面變形和地殼內(nèi)部的構(gòu)造運動及其產(chǎn)生的各種地質(zhì)災(zāi)害都與地殼應(yīng)力作用密切相關(guān)，地殼應(yīng)力狀態(tài)的變化是導(dǎo)致斷裂、褶皺乃至發(fā)生地震的最直接動因。鉆孔應(yīng)變觀測通過對地層內(nèi)部應(yīng)變狀態(tài)依時間連續(xù)變化的精細(xì)觀測，發(fā)現(xiàn)和掌握地震應(yīng)變前兆的(長)中期-短期-臨震以及震后調(diào)整的時空分布與發(fā)展變化規(guī)律，是重要的地震前兆觀測手段。鉆孔應(yīng)變觀測憑借精度高、抗干擾能力強、多測項觀測、儀器易于維護(hù)等特點成為了僅次于gps觀測的前兆手段。

主成分分析方法理論上可以將原始信號中相對微弱的信號從較強干擾背景中識別出來，解決了鉆孔應(yīng)變觀測中地殼變化微弱信號淹沒在較強干擾背景中的問題。

cn104390733a公開了一種地應(yīng)力大小和方向的確定方法，在測量應(yīng)力的環(huán)氧樹脂空心包體上沿圓周方向設(shè)置三個應(yīng)變花，相鄰應(yīng)變花的夾角為120°，每個應(yīng)變花上有四個應(yīng)變片，每個應(yīng)變花上的應(yīng)變片依次旋轉(zhuǎn)45°設(shè)置；利用圓周方向上12個不同方向應(yīng)變值，能夠準(zhǔn)確有效的計算出地應(yīng)力和應(yīng)力方向。

cn202452947u公開了一種四分量鉆孔應(yīng)變儀觀測系統(tǒng)，包括井上設(shè)置的數(shù)據(jù)存儲網(wǎng)絡(luò)傳輸設(shè)備、井下設(shè)置的主四分量應(yīng)變儀及連接井上井下的信號線纜，還包括井下設(shè)置的輔助四分量應(yīng)變儀。通過在觀測系統(tǒng)內(nèi)設(shè)置一主、一輔兩個四分量應(yīng)變儀，增加了系統(tǒng)的可靠性，避免了因一個裝置損壞就使得整個系統(tǒng)癱瘓的情形。

邱澤華等(用小波-超限率分析提取寧陜臺汶川地震體應(yīng)變異常，2012)利用小波分解和改進(jìn)的超限率分析方法對姑咱臺站四分量鉆孔應(yīng)變數(shù)據(jù)進(jìn)行研究，將細(xì)微的地震異常變化精細(xì)的提取出來。池順良等(2013年廬山ms7.0地震的震前及臨震應(yīng)變異常，2013)利用廬山地震前后姑咱臺站四分量鉆孔應(yīng)變觀測資料，對地震同震信號進(jìn)行分析。但到目前為止，尚未見運用主成分分析的方法對鉆孔應(yīng)變數(shù)據(jù)進(jìn)行異常提取的報道。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的就在于針對上述現(xiàn)有技術(shù)的不足，提供一種基于主成分分析的鉆孔應(yīng)變數(shù)據(jù)異常提取方法。

本發(fā)明的目的是通過以下技術(shù)方案實現(xiàn)的：

本發(fā)明利用主成分分析對鉆孔應(yīng)變數(shù)據(jù)進(jìn)行異常提取，首先是將同一臺站的鉆孔應(yīng)變數(shù)據(jù)序列進(jìn)行應(yīng)變換算，對換算后的數(shù)據(jù)進(jìn)行預(yù)處理；將預(yù)處理后的鉆孔應(yīng)變數(shù)據(jù)構(gòu)造成一個矩陣；并對每一天的矩陣進(jìn)行主成分分析，以獲取每個矩陣的特征值和特征向量；將得到的特征值與計算出來的特征向量角度與地震事件相對應(yīng)，以得到異常的檢測結(jié)果。通過本發(fā)明能夠有效的利用主成分分析的方法對鉆孔應(yīng)變數(shù)據(jù)進(jìn)行分析，根據(jù)鉆孔應(yīng)變各測項的相關(guān)性，對可能的地震前兆異常進(jìn)行提取。

一種基于主成分分析的鉆孔應(yīng)變數(shù)據(jù)異常提取方法，包括以下步驟：

a、錄入鉆孔應(yīng)變數(shù)據(jù)，并進(jìn)行數(shù)據(jù)有效性的驗證，“是”，進(jìn)行下一步；

b、對鉆孔應(yīng)變數(shù)據(jù)進(jìn)行應(yīng)變換算；

c、數(shù)據(jù)預(yù)處理；

d、制作樣本數(shù)據(jù)；

e、求樣本數(shù)據(jù)的協(xié)方差矩陣；

f、求出協(xié)方差矩陣的特征值和特征向量；

g、求出特征值向量的角度；

h、輸出特征值變化曲線；

i、輸出特征向量角度變化曲線及統(tǒng)計直方圖。

步驟a，所述的錄入鉆孔應(yīng)變數(shù)據(jù)是選取一個地區(qū)不同臺站的四分量鉆孔應(yīng)變數(shù)據(jù)，制作成按照分鐘值采樣的時間序列。數(shù)據(jù)有效性的驗證是按照不同分量，數(shù)據(jù)可表示為s1,s2,s3,s4；根據(jù)由帶孔平板彈性理論模型導(dǎo)出的圓孔徑向變形與區(qū)域應(yīng)力的關(guān)系式，對鉆孔應(yīng)變數(shù)據(jù)進(jìn)行自洽分析，其關(guān)系式為：

s1+s3＝k(s2+s4)，(1)

選取自洽系數(shù)k在0.9以上的數(shù)據(jù)為有效數(shù)據(jù)。

步驟b，所述的對鉆孔應(yīng)變數(shù)據(jù)進(jìn)行應(yīng)變換算是根據(jù)公式(2)將四分量鉆孔應(yīng)變觀測數(shù)據(jù)換算成兩個剪應(yīng)變s13、s24和一個面應(yīng)變sa，

步驟c，所述的數(shù)據(jù)預(yù)處理是對數(shù)據(jù)進(jìn)行調(diào)和分析，用來去除數(shù)據(jù)的周期項；其調(diào)和分析函數(shù)表達(dá)式為：

其中，a0為時間序列的直流分量，m為諧波的次數(shù)，系數(shù)am,bm是權(quán)重因子，表示各次諧波對總序列的貢獻(xiàn)；原始數(shù)據(jù)減去擬合后的周期數(shù)據(jù)就是我們想要的預(yù)處理后數(shù)據(jù)。

步驟d，所述的制作樣本數(shù)據(jù)是將預(yù)處理后臺站每天的數(shù)據(jù)按時間序列表示為z1＝[z1(1),z1(2),...,z1(1440)],…,zn＝[zn(1),zn(2),...,zn(1440)]，可以得到樣本矩陣y＝[z1,z2,z3,...,zn]^t。樣本矩陣y的表達(dá)式為：

步驟e，樣本數(shù)據(jù)y(n×1440)其協(xié)方差矩陣cy(n×n)的元素γpq由公式(5)計算得到，

其中，和分別為第i行第p和第q分鐘數(shù)據(jù)；和分別是n行數(shù)據(jù)的第p和第q分鐘數(shù)據(jù)的均值。

步驟f，求出協(xié)方差矩陣的特征值是對矩陣y的協(xié)方差矩陣cy進(jìn)行特征值分解，

cy＝r∧r^t，(6)

式中，∧(n×n)為大小降序排列的特征值對角矩陣；r(n×1)為與特征值對角矩陣相對應(yīng)的特征向量矩陣。特征值表示為λ1,λ2,...,λn(λ1＞λ2＞...＞λn)。

步驟g，所述求出特征向量角度是設(shè)特征向量為r＝[v1v2...vn]^t，由圖2所示，特征向量角度為：

與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果在于：本發(fā)明基于主成分分析的鉆孔應(yīng)變數(shù)據(jù)異常提取方法，利用主成分分析中的特征值和特征向量角度分別將地殼的微弱變化表征出來；實現(xiàn)了在有較強背景干擾的情況下對鉆孔應(yīng)變數(shù)據(jù)異常的精確提取。

附圖說明

圖1為基于主成分分析的鉆孔應(yīng)變數(shù)據(jù)異常提取方法流程圖；

圖2為特征向量角度示意圖；

圖3為姑咱地震前兆監(jiān)測臺站與蘆山地震震中位置示意圖；

圖4為第一主成分特征值變化曲線；

圖5為第一主成分特征向量角度變化示意圖；

圖6a為全部時間段的特征向量角度三維直方圖；

圖6b為無震期間特征向量角度三維直方圖；

圖6c為震前5個月地震前兆異常數(shù)據(jù)的特征向量角度三維直方圖；

圖6d為地震前后異常數(shù)據(jù)的特征向量角度三維直方圖。

具體實施方式

下面結(jié)合實施實例對本發(fā)明進(jìn)行詳細(xì)說明。

針對蘆山地震，以四川地區(qū)姑咱地震前兆監(jiān)測臺站的鉆孔應(yīng)變分鐘值數(shù)據(jù)為例。姑咱臺站與蘆山地震位置如圖3所示。該數(shù)據(jù)于2012年1月1日至2013年12月31日使用yry四分量鉆孔應(yīng)變儀測得。

a、錄入姑咱臺站2012年1月1日至2013年12月31日鉆孔應(yīng)變分鐘值時間序列，按照北南分量、東西分量、北東分量、北西分量的順序?qū)?shù)據(jù)分別記為s1、s2、s3、s4；根據(jù)由帶孔平板彈性理論模型導(dǎo)出的圓孔徑向變形與區(qū)域應(yīng)力的關(guān)系式，對鉆孔應(yīng)變數(shù)據(jù)進(jìn)行自洽分析，其關(guān)系式為：

s1+s3＝k(s2+s4)，(1)

計算出的k值，由k≥0.9可知姑咱臺站2012年1月1日至2013年12月31日鉆孔應(yīng)變分鐘值數(shù)據(jù)有效。

b、對鉆孔應(yīng)變數(shù)據(jù)進(jìn)行應(yīng)變換算是根據(jù)公式(2)將四分量鉆孔應(yīng)變觀測數(shù)據(jù)換算成兩個剪應(yīng)變s13、s24和一個面應(yīng)變sa；

c、數(shù)據(jù)預(yù)處理，即對換算后的數(shù)據(jù)進(jìn)行調(diào)和分析，去掉應(yīng)變數(shù)據(jù)的周期響應(yīng)；調(diào)和分析函數(shù)表達(dá)式為：

其中a0為時間序列的直流分量，m為諧波的次數(shù)，系數(shù)am,bm是權(quán)重因子，表示各次諧波對總序列的貢獻(xiàn)。原始數(shù)據(jù)減去擬合后的周期數(shù)據(jù)就是我們想要的預(yù)處理后數(shù)據(jù)。

d、將預(yù)處理后的鉆孔應(yīng)變每天的數(shù)據(jù)表示為

zn＝[zn(1),zn(2),...,zn(1440)],n＝1,2,3，

可以得到樣本矩陣y＝[z1,z2,z3]^t。樣本矩陣y的表達(dá)式為：

e、求樣本數(shù)據(jù)y的協(xié)方差矩陣，樣本數(shù)據(jù)y(3×1440)其協(xié)方差矩陣cy(3×3)的元素γpq由公式(5)計算得到，

其中，和分別為第i行數(shù)據(jù)的第p和第q分鐘數(shù)據(jù)；和分別是第n行數(shù)據(jù)的第p和第q分鐘數(shù)據(jù)的均值。

f、對矩陣y的協(xié)方差矩陣cy進(jìn)行特征值分解，

cy＝r∧r^t，(6)

式中，∧(3×3)為大小降序排列的特征值對角矩陣；r為與特征值對角矩陣相對應(yīng)的特征向量矩陣，特征值表示為λ1,λ2,λ3(λ1＞λ2＞λ3)。

g、求出特征向量角度；以第一主成分特征向量為例設(shè)r1＝[v1,v2,v3]，根據(jù)圖2所示，由公式(7)求出

第二、第三主成分特征向量同理可求。

h、得到特征值變化曲線，如圖4所示，圖中，實線為第一主成分特征值，---虛線為蘆山地震時刻，矩形框中的為提取到的地震相關(guān)異常，可以看出特征值很好的表征了地震相關(guān)異常。

i、得到特征向量角度變化圖及統(tǒng)計直方圖，如圖5，圖6a-圖6d所示。圖5中空心圓點為特征向量角度，---虛線為蘆山地震時刻，矩形框中的為提取到的地震相關(guān)異常。圖6a-圖6d為第一主成分特征向量角度變化三維直方圖，從圖6d可以看出，在異常部分的特征向量角度有明顯的聚集現(xiàn)象，很好的提取到了地震前兆異常。