史鵬飛+張翔+豐繼林+王茂發(fā)+韓定良

摘 要：隨著地震觀測(cè)技術(shù)與手段的不斷更新，數(shù)字化地震觀測(cè)逐漸取代傳統(tǒng)的模擬地震觀測(cè)。為了保存歷史數(shù)據(jù)，如何有效提取模擬地震記錄中的珍貴信息是地震工作者的研究重點(diǎn)。該項(xiàng)目提出了基于OpenCV的模擬地震記錄矢量化算法，使用C#開發(fā)相應(yīng)軟件。該軟件利用有震波形區(qū)域提取方式可以快速、準(zhǔn)確地將模擬地震記錄矢量化。

關(guān)鍵詞：模擬地震記錄；數(shù)字化；OpenCV；有震波形區(qū)域提??；波形追蹤

0 引言

世界上第一張地震圖要追溯到1889年，英國(guó)人米爾恩和尤因利用安置在德國(guó)波茨坦的現(xiàn)代地震儀記錄下了發(fā)生在日本的一次地震。此后70年的時(shí)間里，人類觀測(cè)的地震信息都是以模擬地震記錄的方式存在。模擬地震記錄蘊(yùn)含著大量寶貴的地震信息，具有不可估量的價(jià)值。直到19世紀(jì)70年代，美國(guó)率先建成了由十個(gè)數(shù)字化臺(tái)站組成的數(shù)字化臺(tái)網(wǎng)。隨后數(shù)字地震記錄技術(shù)逐漸取代模擬地震記錄技術(shù)，一直發(fā)展沿用至今。這造成了地震數(shù)據(jù)應(yīng)用研究中，模擬地震信息與數(shù)字地震信息之間有一條明顯的數(shù)據(jù)斷裂。如今大多數(shù)模擬地震記錄存放在倉(cāng)庫(kù)里經(jīng)受著自然的侵蝕，蘊(yùn)藏著寶貴信息的模擬地震記錄無法被計(jì)算機(jī)處理。若能有效提取這些模擬地震信息，將會(huì)推動(dòng)地震研究工作快速發(fā)展。

1 研究現(xiàn)狀

目前為止，對(duì)于模擬地震記錄的研究主要分為兩個(gè)階段：存儲(chǔ)與矢量化。

存儲(chǔ)是將模擬地震記錄中的信息由易損壞的紙質(zhì)信息轉(zhuǎn)化為方便存儲(chǔ)與備份的數(shù)字圖片信息。2012年開始，哈佛大學(xué)投入大量精力進(jìn)行相關(guān)模擬地震記錄電子存檔工作，該項(xiàng)計(jì)劃被命名為哈佛大學(xué)地震波形歸檔計(jì)劃。2015年，河北地震局蔣宏毅等人研發(fā)了一套模擬地震記錄數(shù)字化存儲(chǔ)管理系統(tǒng)，實(shí)現(xiàn)了模擬地震記錄的數(shù)字化存儲(chǔ)與查詢。但數(shù)字存儲(chǔ)只是將紙質(zhì)模擬地震記錄掃描為數(shù)字圖片，沒有進(jìn)行矢量化，計(jì)算機(jī)仍無法識(shí)別和處理模擬地震記錄，也無法挖掘記錄中的寶貴信息‘5]。

許多科研工作者進(jìn)行模擬地震記錄矢量化的研究：Teves-Costa（ 1999）使用CADcore和AutoCAD開展模擬地震記錄的矢量化工作；Pintore（2005）基于人工神經(jīng)網(wǎng)絡(luò)提出了Teseo方法，可以手動(dòng)和自動(dòng)矢量化模擬地震記錄；徐濤等（2014）利用Matlab GUI開發(fā)了一個(gè)交互的模擬地震矢量化程序；王茂發(fā)等（2016）提出了CSF算法，并開發(fā)了相應(yīng)的程序。雖然矢量化研究工作取得了顯著成績(jī)，但許多矢量化工作只是將一小部分含有復(fù)雜波形的圖片截取出來單獨(dú)進(jìn)行矢量化研究，脫離了整張圖片使得后期波形拼接與時(shí)間拼接變得難以實(shí)現(xiàn)。 鑒于跨平臺(tái)計(jì)算機(jī)視覺庫(kù)OpenCV在人臉識(shí)別、車牌識(shí)別、動(dòng)態(tài)手勢(shì)檢測(cè)等方面的準(zhǔn)確性與高效性，因此，本文基于OpenCV設(shè)計(jì)了新的模擬地震記錄矢量化算法，并使用程序設(shè)計(jì)語(yǔ)言C#開發(fā)軟件。其中包含分塊提取波動(dòng)區(qū)域、自動(dòng)提取復(fù)雜波等特有的功能模塊，能夠快速、準(zhǔn)確地處理整張模擬地震記錄。

2 軟件矢量化過程

2.1 流程

本軟件實(shí)現(xiàn)模擬地震記錄矢量化的過程共分為8個(gè)步驟（圖1所示）：加載地震圖紙、圖紙二值化、搜索起點(diǎn)、有震波形區(qū)域提取、掃描波動(dòng)區(qū)域、判斷波形類型、拼接和反演。下面給出每一個(gè)步驟的具體描述。

2.2加載地震圖紙

將整張圖紙加載至矢量化軟件中是模擬地震記錄矢量化最基本的工作。早期的模擬監(jiān)測(cè)是將一段時(shí)間（一般為1天）內(nèi)3個(gè)方向（東西、南北、上下）的地震波動(dòng)記錄在一張圖紙上。觀察整張圖（如圖2）紙可知，波形密集，包含的信息量巨大，該項(xiàng)目使用C#常用控件PictureBox實(shí)現(xiàn)了軟件的圖片加載功能，可以加載、處理整張模擬地震記錄。

2.3 圖紙二值化

許多因素如光線、圖紙的新舊程度等，會(huì)對(duì)模擬地震記錄掃描而成的圖片產(chǎn)生影響，所以必須對(duì)圖紙進(jìn)行二值化處理，它是在圖片預(yù)處理階段減少圖片噪聲非常有效的方式。一張灰度的光柵圖片類似一個(gè)二維的數(shù)組，每一個(gè)元素在數(shù)組中的值在0-255之間。利用OpenCV的Threshold函數(shù)對(duì)圖片進(jìn)行二值化處理：設(shè)置一個(gè)閾值（目前實(shí)驗(yàn)獲得經(jīng)驗(yàn)數(shù)值為155），將大于閾值的像素點(diǎn)置為255，小于閾值的置為0，如圖3所示：

2.4搜索起點(diǎn)

為了確保數(shù)據(jù)的一致性，軟件在提取模擬地震記錄中的信息時(shí)，始終在同一個(gè)坐標(biāo)系下處理。以圖片的左上角為原點(diǎn)，圖片的上邊緣為x軸正值方向，圖片的左邊緣為y軸正值方向。取原點(diǎn)附近的一個(gè)點(diǎn)c（x，0）為起點(diǎn)，沿著垂直方向向下搜索，每次步長(zhǎng)為一個(gè)像素點(diǎn)。搜素算法具體描述如下：

（I）水平方向

水平方向x值保持不變（公式1），為了避免圖片左側(cè)空白邊緣對(duì)搜索結(jié)果的影響，一般x取值稍大于0（實(shí)驗(yàn)取值為10）。 X_i=x（i=l，2…m）

（1）

（2）垂直方向

垂直方向y值由0不斷變大，每次以一個(gè)步長(zhǎng)即一個(gè)像素向下搜索，并判定像素值的變化。當(dāng)遇到相鄰兩個(gè)像素點(diǎn)值跳變時(shí)記錄這個(gè)像素點(diǎn)坐標(biāo)到數(shù)組中（x，Y_n1），然后繼續(xù)向下搜索，再次遇到跳變時(shí)繼續(xù)記錄（x，Y_n2）…。

則第一條波的起點(diǎn)坐標(biāo)為公式（2）：

直到搜索y值到圖紙的下邊緣，數(shù)組中記錄的點(diǎn)（x，Y_n1），（x，Y_n2），……，（x，_nm）。（點(diǎn)的編號(hào)為1，2…m）這樣每一條波形起點(diǎn)的信息都可以確定。第h（l≤h≤m/2）條波的起點(diǎn)坐標(biāo)信息為（公式3）：

2.5有震波形區(qū)域提取

模擬地震記錄中最重要的數(shù)據(jù)是記錄地震發(fā)生時(shí)的復(fù)雜波，矢量化工作的難點(diǎn)是矢量化這震波形區(qū)域的波形。通過觀察發(fā)現(xiàn)這些復(fù)雜波只是整張圖紙很小的一個(gè)局部，所以采取有震波形區(qū)域提取的方式：先找到復(fù)雜波所在的第一個(gè)區(qū)域，在這個(gè)波動(dòng)區(qū)域內(nèi)進(jìn)行矢量化操作，接著尋找復(fù)雜波所在的第二個(gè)區(qū)域進(jìn)行矢量化操作，直至完成所有波動(dòng)區(qū)域的波形矢量化（圖4）。所有的矢量化工作基于圖紙的局部卻沒有脫離整體圖紙，確保了波形信息與時(shí)間信息拼接操作的可行性。如果在整張模擬地震圖紙上進(jìn)行矢量化復(fù)雜波的操作，這樣既增加難度，也降低效率和準(zhǔn)確率。endprint

2.6掃描波動(dòng)區(qū)域

2.6.1 判斷波形類別

首先，使用搜索起點(diǎn)算法將波動(dòng)區(qū)域的所有波形起點(diǎn)找到。然后，按照Y值由小變大的順序以每條波的起點(diǎn)水平掃描這條波。掃描算法具體描述如下：

（1）選取在此區(qū)域內(nèi)第一條波的起點(diǎn)（X_c1，y_c2）。

（2）水平方向，X_C1不變。垂直方向，以一個(gè)像素點(diǎn)為步長(zhǎng)，先沿著v軸正方向搜索至（ X_C1，Y_C1+△y） （Ay最大為波形寬度的平均值，若在此區(qū)域內(nèi)出現(xiàn)像素跳變點(diǎn)則記錄其坐標(biāo)（ X_c1.y_D1），若不出現(xiàn)跳變點(diǎn)則記錄點(diǎn)（ X_C1，y_D1）=（X_C1，Y_C1+△y）；后沿著y軸負(fù)方向搜索至（X_C1，Y_C1-△y），若在此區(qū)域內(nèi)出現(xiàn)像素跳變點(diǎn)則記錄其坐標(biāo)（X_C1，y_D2），若不出現(xiàn)跳變點(diǎn)則（X_C1，y_D2）=（X_C1，Y_C1-△y）。然后記錄在（X_C1，Y_C1±△y）范圍內(nèi)中點(diǎn)坐標(biāo)同時(shí)記錄該點(diǎn)的像素值。

（3）水平方向：X_C2=X_C1+1。

（4）重復(fù)（2-3）操作直至波動(dòng)區(qū)域邊界。

使用掃描算法掃描平滑波時(shí)，由于Ay的限制，可以基本消除與復(fù)雜波交叉的影響，追蹤到平滑波的中線，所以記錄下來平滑波中線的像素基本都為黑點(diǎn)。但由于復(fù)雜波是有振幅的，而△_y值最大為線寬，所以掃描復(fù)雜波時(shí)許多中點(diǎn)的像素點(diǎn)為白色（圖5）。因此，當(dāng)掃描一條波形記錄下的白色像素點(diǎn)多于其他波形時(shí)，則判定該波為復(fù)雜波同時(shí)進(jìn)行標(biāo)記。

2.6.2 追蹤平滑波

判斷波形類型后，根據(jù)標(biāo)記信息進(jìn)行平滑波的跟蹤。利用改進(jìn)的掃描算法追蹤平滑波：掃描時(shí)，若在區(qū)域（y+△y，y-△y）}H現(xiàn)像素跳變點(diǎn)，記錄此范圍內(nèi)的中點(diǎn)并將此范圍的所有黑色像素點(diǎn)置為白色像素點(diǎn)，反之則不變。這樣既保證記錄平滑坡中線坐標(biāo)的正確性，也不會(huì)破壞復(fù)雜波的波形，追蹤效果如圖5所示：

2.6.3 去除噪聲

平滑波追蹤完成后，復(fù)雜波追蹤變得更加簡(jiǎn)單。但儀器記錄和保存時(shí)各種因素可能會(huì)導(dǎo)致圖片存在干擾信息，為了避免這些噪聲對(duì)追蹤復(fù)雜波的影響，必須再次進(jìn)行去噪處理。首先，利用OpenCV對(duì)圖片進(jìn)行SmoothMedian處理，消除椒煙噪聲與斑點(diǎn)噪聲。其次，利用FindContours函數(shù)提取所有圖形的輪廓，記錄每一個(gè)輪廓的id，計(jì)算每一個(gè)輪廓的周長(zhǎng)。最后，通過比較保留周長(zhǎng)最長(zhǎng)的輪廓，其余全部去除——這樣就保存了復(fù)雜波信息而去除了噪聲，圖6是保留的復(fù)雜波輪廓。

2.6.4追蹤復(fù)雜波

如同上面一樣，可以計(jì)算出復(fù)雜波起始處的中點(diǎn)坐標(biāo)，然后以此點(diǎn)為起點(diǎn)追蹤復(fù)雜波。在追蹤過程中，值得注意的是，在波峰和波谷處需要記錄此處的v值上邊緣值與下邊緣值，這樣可以正確反映地震的能量釋放情況。追蹤復(fù)雜波的算法與追蹤平滑波的算法相似，不同之處在于：沒有△y的限制，當(dāng)中間點(diǎn)的y值大于相鄰兩點(diǎn)的v值時(shí)（波谷，v軸正方向向下），記錄其下邊緣的v值；當(dāng)中間點(diǎn)的v值小于相鄰兩點(diǎn)的v值日寸（波峰），記錄其上邊緣的v值；搜素到區(qū)域邊界，圖7中紅線為波形的中線，藍(lán)點(diǎn)為提取的關(guān)鍵點(diǎn)。

2.7 拼接和反演

在前期工作中已經(jīng)實(shí)現(xiàn)了拼接與反演，這里只做簡(jiǎn)單描述。每條波形上的點(diǎn)都以坐標(biāo)的形式存儲(chǔ)在數(shù)組中并且標(biāo)有唯一的id，首先根據(jù)點(diǎn)的坐標(biāo)信息進(jìn)行波形的拼接，然后根據(jù)圖紙上標(biāo)注的時(shí)間信息進(jìn)行時(shí)間的拼接，最后將波形信息反演，如圖8所示：

3 結(jié)束語(yǔ)

該項(xiàng)目設(shè)計(jì)了基于OpenCV的模擬地震記錄矢量化軟件，分8個(gè)步驟矢量化模擬地震記錄的柵格圖片，通過給出的測(cè)試結(jié)果可以看出，反演出來的波形較好的符合原始圖紙中的波形。下一步，將進(jìn)一步提升算法效率和改善軟件，期望開發(fā)出自動(dòng)校驗(yàn)?zāi)K以確保矢量化的質(zhì)量。endprint