基于OpenCV的模擬地震記錄矢量化軟件

史鵬飛+張翔+豐繼林+王茂發+韓定良

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

關鍵詞：模擬地震記錄；數字化；OpenCV；有震波形區域提取；波形追蹤

0 引言

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

1 研究現狀

目前為止，對于模擬地震記錄的研究主要分為兩個階段：存儲與矢量化。

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

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

2 軟件矢量化過程

2.1 流程

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

2.2加載地震圖紙

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

2.3 圖紙二值化

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

2.4搜索起點

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

（I）水平方向

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

（1）

（2）垂直方向

垂直方向y值由0不斷變大，每次以一個步長即一個像素向下搜索，并判定像素值的變化。當遇到相鄰兩個像素點值跳變時記錄這個像素點坐標到數組中（x，Y_n1），然后繼續向下搜索，再次遇到跳變時繼續記錄（x，Y_n2）…。

則第一條波的起點坐標為公式（2）：

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

2.5有震波形區域提取

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

2.6掃描波動區域

2.6.1 判斷波形類別

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

（1）選取在此區域內第一條波的起點（X_c1，y_c2）。

（2）水平方向，X_C1不變。垂直方向，以一個像素點為步長，先沿著v軸正方向搜索至（ X_C1，Y_C1+△y） （Ay最大為波形寬度的平均值，若在此區域內出現像素跳變點則記錄其坐標（ X_c1.y_D1），若不出現跳變點則記錄點（ X_C1，y_D1）=（X_C1，Y_C1+△y）；后沿著y軸負方向搜索至（X_C1，Y_C1-△y），若在此區域內出現像素跳變點則記錄其坐標（X_C1，y_D2），若不出現跳變點則（X_C1，y_D2）=（X_C1，Y_C1-△y）。然后記錄在（X_C1，Y_C1±△y）范圍內中點坐標同時記錄該點的像素值。

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

（4）重復（2-3）操作直至波動區域邊界。

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

2.6.2 追蹤平滑波

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

2.6.3 去除噪聲

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

2.6.4追蹤復雜波

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

2.7 拼接和反演

在前期工作中已經實現了拼接與反演，這里只做簡單描述。每條波形上的點都以坐標的形式存儲在數組中并且標有唯一的id，首先根據點的坐標信息進行波形的拼接，然后根據圖紙上標注的時間信息進行時間的拼接，最后將波形信息反演，如圖8所示：

3 結束語

該項目設計了基于OpenCV的模擬地震記錄矢量化軟件，分8個步驟矢量化模擬地震記錄的柵格圖片，通過給出的測試結果可以看出，反演出來的波形較好的符合原始圖紙中的波形。下一步，將進一步提升算法效率和改善軟件，期望開發出自動校驗模塊以確保矢量化的質量。endprint