李 穩(wěn),吳 晨,竇海岳,呂子強

(1.中國地震局地球物理勘探中心,鄭州 450002;

2.山東省地震局,濟南 250014)

地震波視出射角求取軟件的C++實現(xiàn)

李 穩(wěn)1,吳 晨2,竇海岳2,呂子強2

(1.中國地震局地球物理勘探中心,鄭州 450002;

2.山東省地震局,濟南 250014)

在分析數(shù)字地震波形記錄數(shù)據(jù)格式和信號處理流程的基礎上,選用C++作為編程語言開發(fā)了專用的地震波視出射角求取軟件。該軟件利用了C++語言的優(yōu)勢,具有高性能數(shù)值計算能力和高級交互式界面設計。軟件中使用的波形資料處理流程對地震信號處理研究有參考價值。具體算法均在處理大量實際資料的過程中得到了嚴格檢驗,高效可靠。

地震波視出射角;聯(lián)合定位;地震波理論與應用;C++

0 引言

通過充分挖掘地震波所攜帶的豐富信息,增加獨立的觀測參數(shù)以減輕地震定位反問題的多解性是地震科學研究中的重要課題。地震波視出射角參數(shù)是一個獨立的觀測參數(shù)。理論計算和實際地震資料處理結果表明,利用地震波視出射角參數(shù)與到時數(shù)據(jù)進行聯(lián)合定位可有效提高地震定位對震源深度的分辨,從而提高地震定位精度[1]。

本文在分析數(shù)字地震波形數(shù)據(jù)格式和處理流程的基礎上,選用C++作為編程語言,在Visual C++6.0和MFC 6.0環(huán)境下編寫了ProcessDSR可視化交互分析軟件。該軟件專用于求取數(shù)字地震記錄中的地震波視出射角參數(shù)。

1 地震波視出射角的概念

人們在地球表面觀測和接收地震波。地震波入射射線與地面的夾角是真出射角,常用e表示。地面位移矢量與地面的夾角稱為視出射角,常用ˉe表示[2]。視出射角的計算公式為公式(1)。視出射角與真出射角的關系式為公式(2)。圖1是視出射角與真出射角示意圖。

圖1 視出射角與真出射角示意圖

式中,AZ為地動位移的垂直分量,ANS為N S分量,AEW為EW分量。VP為縱波速度,VS為橫波速度,其比值稱為介質的泊松比。

2 數(shù)字地震波形記錄數(shù)據(jù)格式

數(shù)字地震波形記錄以特定的數(shù)據(jù)格式存儲,常見的有 SEED、SAC、EVT等。甘肅省地震臺網(wǎng)提供的波形記錄采用EVT格式。

EVT格式包含事件頭段結構、臺網(wǎng)參數(shù)塊結構、臺站參數(shù)結構、通道參數(shù)結構、索引塊結構和數(shù)據(jù)塊結構幾個部分。地震波原始數(shù)據(jù)樣本經(jīng)STEIM2壓縮算法壓縮后保存在數(shù)據(jù)塊結構中[3]。

3 數(shù)字地震波處理流程

臺網(wǎng)記錄到的地震波原始信號包含著一些干擾因素,在求視出射角參數(shù)之前需對其進行扣除儀器響應、濾波等常規(guī)處理。本文參考前人的研究成果[4],設計了規(guī)范有效的數(shù)字波形資料處理流程,見圖2。

圖2 求取地震波視出射角參數(shù)的流程圖

4 軟件實現(xiàn)

C++是靜態(tài)類型語言的典型代表。在高級系統(tǒng)程序設計(如用戶界面設計、圖像和聲音的編排等)和科學計算等領域有明顯的優(yōu)勢[5]。地震波視出射角參數(shù)求取軟件既需實現(xiàn)復雜、高效的數(shù)字信號處理算法,又要提供高性能的用戶界面設計。經(jīng)過分析討論,我們選用C++作為編程語言進行軟件開發(fā)。在Visual C++6.0和MFC 6.0環(huán)境下編寫、編譯,在win32平臺上進行系統(tǒng)測試。

4.1 數(shù)據(jù)讀取

在EVT數(shù)據(jù)格式中存在由某些字段值決定內存分配的情況,因此在定義類時需用到一級和二級指針。為了安全,必須禁止類對象之間的按值傳遞[6]。另外,為避免內存泄露,在使用指針前將其初始化,在不再使用時釋放其所指向的內存并將其復位為0,同樣非常重要。

設計了CDSR類用來存儲EVT波形數(shù)據(jù)。該類具有基本的數(shù)據(jù)讀入、提取功能,并通過聲明私有拷貝構造函數(shù)的方法禁止按值傳遞。

ReadDSR()函數(shù)的實現(xiàn)利用了 EVT文件將二進制數(shù)據(jù)封裝在文件頭、臺站參數(shù)、波形數(shù)據(jù)塊3個部分的特點。首先建立好嵌套結構;利用fread語句讀取文件頭字段信息;根據(jù)其中提供的臺站總數(shù)計算臺站參數(shù)部分所占內存,然后將臺站參數(shù)部分以數(shù)據(jù)塊的模式讀入;再根據(jù)臺站參數(shù)部分提供的采樣率、記錄長度、數(shù)采字長等信息分配內存,讀取波形數(shù)據(jù)塊。

將數(shù)據(jù)讀入CDSR類的實例之后,可調用DistillData()函數(shù)在地震波初至附近提取4096個采樣點(通常為81.92s)的計算數(shù)據(jù)。該取值既能保證截取的信號具有足夠高的頻率分辨力,又適宜于進行快速傅氏變換和小波變換。

4.2 扣除儀器響應

為了消除波形記錄中各頻段靈敏度的影響,必須扣除儀器響應[7-8]。首先根據(jù)頻率響應參數(shù)進行頻點靈敏度插值(利用友元函數(shù)MyChaZhi()進行計算);然后對提取的波形數(shù)據(jù)做傅氏變換,將所得復數(shù)序列的幅值除以頻點靈敏度插值結果;最后進行反傅氏變換,所得復數(shù)序列的實部即為地震波數(shù)據(jù)扣除儀器響應后的結果。

實現(xiàn)傅氏變換時我們選用了時間抽選奇偶分解傅氏變換算法,代碼如下:

4.3 濾波

長周期低頻信號和高頻干擾信號對地震波有效信號有著嚴重的影響,可設計數(shù)字濾波器將它們?yōu)V除。

Ⅰ 采用窗函數(shù)法設計線性相位的 FIR帶通濾波器[9]。

(1)根據(jù)實際問題確定技術指標

選取初始的濾波器通帶下邊緣截止頻率為2.5Hz,通帶上邊緣截止頻率為20.0Hz,過渡帶寬為1.25Hz,這些參數(shù)以后可在交互分析界面中進行調整。要求阻帶衰減不小于-60dB。

(2)用因果離散時間系統(tǒng)逼近能滿足上述技術指標的濾波器

用有限長的單位沖擊響應 h(n)逼近無限長的理想單位沖擊響應hd(n)的有效方法是用一個有限長的窗函數(shù)序列 w(n)來截取 hd(n)?？紤]到幅度要求指標(阻帶衰減不小于-60dB),我們選用了Blackman窗作為窗函數(shù),其阻帶衰減為-74dB。最終計算出的 h(n)的表達式為:

式中,ωc1=0.075π,ωc2=0.825π,它們是兩個理想低通濾波器的截止頻率;N=221,是單位沖擊響應 h(n)的長度;α=110,是為滿足線性相位響應所必需的移位值。

(3)系統(tǒng)實現(xiàn)

得到了以單位沖擊響應表示的系統(tǒng)表達式以后,濾波器既可以用硬件來實現(xiàn),還可以在計算機上通過編程用軟件實現(xiàn)。

實現(xiàn)濾波器的算法很多,本軟件采用了由Grant R.Griffin 發(fā)布的fir—split算法(文件fir—algs—1-0.c[10]),該算法運算速度快,且避免了使用環(huán)形緩存。

Ⅱ 小波分析濾波

考慮到地震波信號頻率成分隨時間變化的特點,進一步采用了小波分析濾波。小波分析濾波可通過在二進小波分析(常用Mallat二進小波算法)中適當?shù)剡x擇尺度因子來實現(xiàn)[11-13]。

4.4 傾斜校正與水平基線校正

進行傾斜校正和水平基線校正的目的是去掉記錄中的直流分量。傾斜校正可通過最小二乘法線性回歸實現(xiàn)。最小二乘法線性回歸的計算公式和實現(xiàn)代碼如下:

設回歸直線為:y=a+bx用最小二乘法擬合回歸直線的解為:

水平基線校正即先求得地震記錄序列的平均值,再將該值從整個序列中減去。水平基線校正應在傾斜校正之后進行。其實現(xiàn)代碼如下:

4.5 數(shù)值積分

數(shù)字地震記錄是離散的數(shù)據(jù)集,不適合用牛頓-萊布尼茨公式積分,可用梯形公式(1階的牛頓-柯特斯公式)等數(shù)值積分法積分[14]。其計算公式和實現(xiàn)代碼如下:

4.6 二次濾波

一次積分運算在頻率域中相當于頻譜除以圓頻率[4],所以在積分運算中,原始信號的低頻成分被放大了,甚至將有用的信號淹沒。因此,在積分運算之后還應再次進行濾波處理。

4.7 處理結果屏幕顯示與交互分析

在地震波視出射角求取軟件中除保證數(shù)據(jù)處理功能的穩(wěn)定高效之外,提供一個高性能的交互分析界面也非常重要。

本論文參考“EDSP-IAS”軟件界面設置,應用文檔/視圖結構和非模態(tài)對話框實現(xiàn)交互分析功能。在軟件設計過程中嚴格遵循面向對象編程的思想,軟件的修改、維護非常方便。該軟件提供的具體功能包括:交互選擇計算數(shù)據(jù)、通過下拉組合框設置濾波器參數(shù)、通過單選按鈕選擇顯示速度記錄還是位移記錄、通過復選框選擇是否在屏幕上疊加原始記錄、通過拖動滑動條設置放大倍數(shù)等。計算結果均在屏幕上即時顯示。

5 實際應用

利用該軟件分析處理了2002年12月玉門地震序列的數(shù)字地震波形記錄,共求得247條地震波的視出射角參數(shù)。應用地震波到時與視出射角聯(lián)合定位方法得到74個地震事件的定位結果。將聯(lián)合定位結果與目前工作中常用的雙差法定位結果進行了對比研究。

參照走時殘差的定義,定義視出射角殘差為理論值與觀測值的差值,見公式(3)。則平均視出射角殘差的計算公式為公式(4),視出射角殘差的均方差的計算公式為公式(5)。

表1顯示了不同方法定位結果的殘差分析結果。可見,增加對視出射角參數(shù)利用的聯(lián)合定位結果,具有更小的平均走時殘差和走時殘差均方差。這說明,聯(lián)合定位結果不但整體上更加準確,而且更加穩(wěn)定。

表1 玉門地震序列地震定位殘差分析結果表

圖3、圖4分別是玉門地震序列近垂直于斷層走向的深度剖面圖和斷裂、三維震中分布圖?？梢娐?lián)合定位結果在震源深度方面更加收斂,顯示出了明顯的傾角,能更有效地反映發(fā)震構造的深部信息。據(jù)此勾畫出的斷層性質與前人的地質構造研究成果[15-16]相吻合。

圖3 玉門地震序列深度剖面圖

圖4 玉門地震序列斷裂、三維震源位置分布圖

通過殘差分析和與地質研究成果的對照分析,可以看出增加對視出射角參數(shù)的利用后所得的定位結果更加準確、可靠。

6 結論與討論

地震定位受諸多因素的影響。在現(xiàn)有條件下,通過增加對地震波視出射角(尤其是震中附近臺站監(jiān)測到的地震波視出射角)的利用,獲得更準確的地震深度信息,對研究發(fā)震斷層的性質及其空間展布有實用價值。我國已在全國范圍內建立了數(shù)字地震觀測臺網(wǎng),地震波到時與視出射角聯(lián)合定位方法具有廣泛適用性。

C++語言具有結構規(guī)范、便于調試和類型安全的特點,在數(shù)值/科學計算、通用應用程序設計等領域具有廣泛的應用前景。地震波分析處理軟件既需具備高性能數(shù)值計算能力又要提供高級交互式界面設計,宜選用C++語言開發(fā)實現(xiàn)。

在頻率域中求取視出射角參數(shù),以提高計算效率和準確度是可行的研究方向,相關研究已經(jīng)取得了初步成果[17]。另外,完善該軟件使其可以直接處理多種不同格式(如SEED、SAC等)的地震數(shù)據(jù),也是一個切實的問題。

[1] 李穩(wěn),張元生,何斌.地震波到時與視出射角聯(lián)合定位方法研究[J].西北地震學報,2009,31(3):207-210.

[2] 傅淑芳,劉寶誠.地震學教程[M].北京:地震出版社,1991:83-85.

[3] 馬中華.臺站數(shù)據(jù)共享技術研究[D].甘肅:中國地震局蘭州地震研究所,2008.

[4] 陳運泰,吳忠良,王培德,等.數(shù)字地震學[M].北京:地震出版社,2000:55-65.

[5] Bjarne Stroustrup.C++語言的設計和演化[M].裘宗燕譯.北京:機械工業(yè)出版社,2002:150-154.

[6] Bruce Eckel.C++編程思想(第二版)[M].劉宗田,袁兆山,潘秋菱,等譯.北京:機械工業(yè)出版社,2002:256-257.

[7] 李媛媛,吳東.山西數(shù)字遙測地震臺網(wǎng)十五勘選子臺臺址地動噪聲分析[J].華北地震科學,2004,22(1):60-62.

[8] 楊晶瓊,楊周勝.云南地震數(shù)字遙測臺網(wǎng)子臺噪聲分析[J].地震研究,2005,28(1):86-89.

[9] 陳玉東.數(shù)字信號處理[M].北京:地質出版社,2005:146-148.

[10] Grant R and Griffin.Digital Signal Processing Software[CP/OL].2008-11-28.www.dspguru.com/sw/lib/fir—algs—1-0.c.

[11] 連海寧,謝禮立.三種變換在強震記錄分析方面的研究[J].華北地震科學,2007,25(3):28-33.

[12] 劉澤民.小波分析在數(shù)字地震資料中的應用[J].地震地磁觀測與研究,2004,1(25):88-91.

[13] 曾憲偉.利用小波包變換識別地震和爆破[D].甘肅:中國地震局蘭州地震研究所,2008.

[14] Pallab Ghosh.數(shù)值方法(C++描述)[M].北京:清華大學出版社,2008:362-371.

[15] 何文貴,鄭文俊,趙廣堃,等.2002年12月14日甘肅玉門5.9級地震的發(fā)震構造研究[J].地震地質,2004,26(4):688-697.

[16] 陳柏林,劉建生,張永雙,等.玉門斷裂全新世活動特征及其與玉門地震的關系[J].地質論評,2005,51(2):138-142.

[17] 何斌,張元生,李穩(wěn).計算地震初至波視出射角方法[J].西北地震學報,2010,32(1):11-15.

The C++Realization of Seismic Apparent Emergence Angle Calculation Softw are

LI Wen1,WU Chen2,DOU Hai-yue2,LV Zi-qiang2
(1.The Geophysical Exploration Center of CEA,Zhengzhu 450002,China;2.Earthquake Administration of Shandong Province,Jinan 250014,China)

On the basis of analysis on digital seismic wave data format and signal processing,using the C++language,we developed software for seismic apparent emergence angle calculating.Utilized the advantages of C++,the software has high-performance numerical computing ability and advanced interface design.The seismic wave processing design used in the software is valuable for the study of seismic signal processing.The computerized algorithms had stood a severe proof in the process of processing arithmetic real data.They are efficient and reliable.

apparent emergence angle;multiparameter seismic location;theory and applications of seismic waves;C++

P315.69

A

1003-1375(2010)04-0009-06

2010-04-01

國家自然科學基金(40874029)資助;

李穩(wěn)(1983-),男(漢族),河南開封人,助理工程師,主要從事地震波理論與應用研究.E-mail:liwen—forwork@163.com