圖像信息熵約束的淺地層層界劃分方法

趙建虎，馮 杰，施 鳳, 張紅梅, 何林幫

(1.武漢大學 測繪學院，武漢 430079；2.武漢大學 動力與機械學院，武漢 430072)

圖像信息熵約束的淺地層層界劃分方法

趙建虎1，馮 杰1，施 鳳1, 張紅梅2, 何林幫1

(1.武漢大學 測繪學院，武漢 430079；2.武漢大學 動力與機械學院，武漢 430072)

為實現(xiàn)快速、精確、自動化、智能化的海底淺地層層界提取，克服傳統(tǒng)淺地層層界在復雜海洋環(huán)境下提取時的低效、模糊、主觀性等缺點，提出一種基于圖像信息熵約束的淺地層層界劃分方法．首先，將淺剖圖像分割為不同區(qū)塊；然后，在不同區(qū)塊計算信息熵，并結(jié)合鉆孔數(shù)據(jù)，建立信息熵與顯著性參數(shù)關系模型；最后，據(jù)此模型對整個淺剖圖像進行層界劃分．研究表明,該方法克服了現(xiàn)有方法的不足，實現(xiàn)了淺地層剖面層界的自適應、準確劃分，試驗中取得了與鉆孔層界深度、厚度同量級的精度．由此可知采用圖像信息熵約束進行層界提取，可以實現(xiàn)淺地層層界提取的自動化與智能化．

淺地層剖面圖像；淺地層層界及提??；二維熵；中誤差系數(shù)；自適應層界提取

借助淺剖儀的換能器發(fā)射的聲波在海底不同介質(zhì)面發(fā)生的反射信號，可以獲得淺剖圖像[1-8]．但受復雜海底環(huán)境、多樣性底質(zhì)和聲波在層界間多次反射造成的多次波干擾等影響，清晰的淺剖圖像難以得到，基于一般的圖像處理技術很難實現(xiàn)層界的準確、連續(xù)劃分[9]．為此，本文結(jié)合淺剖Ping(一個完整的聲波發(fā)射—接收過程)回波強度時序特點、多Ping構(gòu)成的圖像的復雜程度以及層界變化的連續(xù)性，提出了一種顧及圖像信息熵[10-12]的淺地層層界劃分方法，以期消除諸因素影響，實現(xiàn)層界的自適應、準確劃分．

1 淺地層剖面圖像生成及其預處理

1.1 圖像生成

淺地層剖面儀(淺剖)測量中，換能器接收到的時序回波強度數(shù)據(jù)常以二進制的SEG-Y格式文件存儲，為獲得淺地層數(shù)字圖像，需首先對其解碼[13-14]．由于接收信號常以相對發(fā)射聲波振幅形式記錄，需借助Hilbert變換才能獲得接收信號的實際振幅或強度[15]．

為提高數(shù)據(jù)質(zhì)量，對Hilbert變換獲得的聲強進行截取有效聲強、剔除無效數(shù)據(jù)、聲強集中部分拉伸等處理．結(jié)合走航過程中每Ping的平面位置、測量中每個回波的深度，在地理框架下即可得到淺地層剖面圖像．

1.2 圖像預處理

受海洋環(huán)境、聲波傳播機制等影響，淺剖換能器接收到的回波數(shù)據(jù)中存在Ping丟失、噪聲和多次波等干擾，導致圖像質(zhì)量偏低，為此需對原始圖像開展如下預處理.

1)Ping修復. 基于底質(zhì)結(jié)構(gòu)變化的連續(xù)性，根據(jù)丟失Ping前后Ping中的有效Ping數(shù)據(jù)，采用均值內(nèi)插或基于底質(zhì)漸變原則的多項式模型內(nèi)插，實現(xiàn)丟失Ping修復．

2)噪聲抑制. 中值濾波具有不損失特征邊緣、對噪聲消除徹底等特點，可用于抑制淺剖圖像中的噪聲，提高圖像質(zhì)量．

3)多次波壓制. 利用多次波的周期性特征，通過設計預測誤差濾波器，基于預測反褶積法削弱多次波影響[16]．

2 層界劃分基本原理

2.1 基于Ping回波時序強度變化特點的層界峰谷劃分方法

根據(jù)水聲學原理，Ping聲信號在進入不同介質(zhì)的層界面時會產(chǎn)生強反射，在淺地層圖像中則表現(xiàn)為灰度突變；而在層界間，由于底質(zhì)近似，回波強度接近，圖像中則表現(xiàn)為灰度均勻變化．根據(jù)以上特點，在每Ping形成的列灰度序列中定義一個由一定數(shù)量采樣點的灰度值組合形成的統(tǒng)計窗口，并給出如下5個層界劃分原則.

1)灰度陡變點為對應統(tǒng)計窗口內(nèi)的峰值點.

2)灰度陡變點有足夠的顯著性．對于第m個檢驗點，定義檢驗模型為

G(m)-μ>kσ.

式中:G(m)為第m個檢驗點的灰度值，統(tǒng)計窗口的范圍為(m-10，m+10);k為中誤差系數(shù);μ、σ分別為統(tǒng)計窗口內(nèi)灰度的均值、標準差.

3)考慮聲能隨傳播距離衰減，后一個峰值點灰度值應小于前一個.

另外由于淺地層剖面測量時受到來自海洋、儀器、多次波等各類噪聲的影響，使得單Ping灰度曲線變化劇烈，存在眾多噪聲點或偽層界．為削弱這些影響，達到提取真實層界的目的，本文采用NPing均值濾波方法，即認為在NPing變化范圍內(nèi)層界是連續(xù)的、結(jié)構(gòu)近似一致．在參與平均的Ping數(shù)選擇時，應遵循如下2個原則.

4)涵蓋的范圍不能太大，避免出現(xiàn)層界細部信息丟失.

5)涵蓋的范圍不能太小，否則不能達到噪聲壓制效果.

基于以上5個原則，對每Ping開展層界劃分,將相鄰Ping相同層界連接，獲得整個圖像中的層界及其分布.

2.2 圖像二維信息熵

以上層界峰谷劃分法中，中誤差系數(shù)k及參與平均的Ping數(shù)N的確定問題需要解決．k和N均與圖像中灰度變化的復雜程度、噪聲、層界圖像的顯著程度等密切相關．為此引入圖像二維信息熵，以期實現(xiàn)上述2個參數(shù)的自適應確定．

2.3N、k與V關系建立

獲得二維熵V后，可以根據(jù)鉆孔取芯數(shù)據(jù)提供的真實層界為參考，結(jié)合對應層界的圖像二維熵V，以及基于不同N和k下層界的劃分結(jié)果，構(gòu)建彼此間的關系模型.

1)N、k和V參數(shù)集確定．為方便計算，首先將整個圖像根據(jù)鉆孔位置劃分成若干區(qū)塊，區(qū)塊大小為S，并提取S對應圖像的二維熵V；然后在一個區(qū)塊內(nèi)，不斷改變N和k，借助層界峰谷劃分法開展層界劃分；以鉆孔數(shù)據(jù)為參考，獲得劃分層界與實際層界一致時的N、k及V；在不同區(qū)塊開展上述工作，最終形成N、k和V構(gòu)成的集合．{N，k，V}集合的尋找過程如圖1所示．

圖1 k，V，N確定流程

2)利用{N，k，V}集合繪制三維圖，根據(jù)N和V在三維圖中的變化曲線，發(fā)現(xiàn)變化拐點，結(jié)合最大信息熵原則，最終確定N．

3)根據(jù)輸出的{N，k，V}序列，分析k隨N的變化及成因．

考慮k、V變化特點，采用高斯擬合函數(shù)建立k-V關系模型為

(1)

式中，ai、bi、ci(i=1,2,…,8)分別為模型系數(shù)．與傳統(tǒng)的多項式擬合相比，高斯擬合采用高斯函數(shù)系，具有收斂快、擬合準確等特點．根據(jù)實際擬合情況，本文采用式(1)所示8項高斯疊加函數(shù)作為k-V關系模型．

基于式(1)給出的k-V關系模型及N-V關系，在不同的子區(qū)，根據(jù)對應的圖像二維熵可得相應的N及k值，據(jù)此借助層界峰谷劃分法實現(xiàn)層界劃分．

3 層界劃分流程

根據(jù)以上原理，基于圖像二維信息熵約束的淺地層層界劃分流程如下,如圖2所示.

圖2 層界劃分流程

4 實驗與分析

為驗證上述方法的正確性，在煙臺某水域借助C-Boom淺地層剖面儀開展了1個航次的淺地層剖面測量和鉆孔取芯工作．測量中，C-Boom的采樣頻率為4 Hz，船速為4～5節(jié)，Ping間隔約為0.5 m．走航斷面平均設有7個鉆孔點，經(jīng)實驗室分析得到所在位置的淺地層層界．

對淺地層剖面觀測數(shù)據(jù)分別開展SEG-Y格式原始文件解碼和轉(zhuǎn)換，并開展Ping修復、消噪和多次波壓制等圖像預處理．預處理后的圖像如圖3所示．

選擇4個代表性鉆孔，在對應位置劃分區(qū)塊，根據(jù)各區(qū)塊的V、N及k構(gòu)建{N，k，V}集合，并繪制如圖4所示圖形．可以看出：1)整體上N越大V越小．分析認為隨著N變大，濾波所得圖像細節(jié)信息損失，導致V減小. 在N為35～40區(qū)間出現(xiàn)不一致現(xiàn)象．該現(xiàn)象主要是由于當N超過圖像區(qū)塊大小S時，方法將S放大進而造成V增大; 2)V越大k越大．主要由于V增大，層界提取時圖像中的信息量就越大，只有選取較大的k才能界定峰點．

圖3 預處理后的淺地層剖面

圖4 {k, V, N}集合三維分布

以上僅定性分析了三者間關系，下面定量確定彼此間關系．考慮直接建立N-k-V關系模型比較復雜，為簡化問題，首先確定N，然后構(gòu)建k-V模型．

從圖4可看出，N≤35時N增大V減?。籒在 36～37區(qū)間時N增大V增大；N>37時N增大V減小．因此，兩者關系在35～40間存在兩個拐點：35和37．考慮層界處圖像灰度會產(chǎn)生突變，信息熵會較大，取N=37，此處V最大，同時還兼顧了N與V之間的兩個關系．

為驗證上述結(jié)論，讓k=2.5，變化N從2到60，對不同N下所得平均Ping灰度序列，借助峰谷法提取層界,如圖5所示．從圖 5可看出：N越小，提取出的層界越不準確；隨著N增大，層界變得越來越接近鉆孔層界位置；但隨N進一步增大，層界位置逐漸變得模糊，最終偏離實際．以上關系變換的拐點即在N=37，與圖4給出的結(jié)論一致．

N確定后，提取出{N，k，V}集中的{k，V}集，將之繪制到圖6中，并利用這些數(shù)據(jù)構(gòu)建k-V關系模型，模型參數(shù)和模型曲線如表1和圖6所示．可以看出，模型內(nèi)符合精度較高，RMSE達到0.006 8；擬合曲線與實際非常吻合，表明確定的模型參數(shù)真實可靠．

計算其他3個未參與建模鉆孔位置對應區(qū)塊的V，依據(jù)上述模型得到各區(qū)塊的N及k．同樣借助圖1所示流程，得到這3個區(qū)塊對應的N0和k0，并將V對應的k0繪制到圖7中．從圖7可以看出，基于模型得到k-V變化曲線與實際k0取得了一致，同時N-V關系得到的N與N0相同，進一步表明基于上述方法得到的N-V關系以及k-V關系模型適用于整個圖像．

圖5 N與層界劃分準確性關系曲線

圖6 k-V關系曲線

參數(shù)類型1階2階3階4階5階6階7階8階a-5.930.677.282.05003701.001.52b5.815.655.82-369.204.564.614.785.03c0.080.050.102442.000.020.060.030

注：模型精度REMS=0.006 8.

得到可靠的N-V關系及k-V關系模型后，便可進行層界劃分．首先進行分塊，S為平均N的2倍；然后計算不同區(qū)塊的V，N和k；采用逐Ping遍歷法，對每一Ping對應的灰度序列基于峰谷法開展層界提取，完成所有區(qū)塊層界劃分；最后，將相鄰層界連接，獲得整個層界．比較圖8和圖3可看出，相對原始圖像，本文方法劃分出的層界真實、清晰．

圖7 k-V關系模型與k0的比較

為定量評估本文方法的正確性及精度，以鉆孔給出的層界為參考，將上述劃分層界與之比較(如圖9所示)，劃分層界位置及層厚度偏差統(tǒng)計結(jié)果如表2所示．從圖9可以看出，二者吻合較好；表2中的統(tǒng)計結(jié)果表明，二者的層界深度偏差最大絕對值小于20 cm，厚度偏差最大絕對值小于30 cm，與鉆孔數(shù)據(jù)的給出的層界精度(分米級)處同量級，表明本文給出的方法不但實現(xiàn)了層界的自適應劃分，層界劃分的精度也滿足了工程需要．

圖8 劃分出的層界分布圖

圖9 用于檢驗層界提取的鉆孔位置

Tab.2 Biases of layer’ position and thickness determined by layer extraction relative to the drilling data

鉆孔層界深度偏差/m層厚度偏差/m層1層2層3層1-2層2-3ZK2-0.100.100.07-0.200.02ZK30.10-0.20—0.30—ZK70.03-0.150.150.18-0.30

5 結(jié) 論

1)淺剖圖像的預處理(包括修復、去噪、壓制多次波等)，可提高圖像信噪比，得到層界較清晰的淺剖圖像．該環(huán)節(jié)得到的圖像質(zhì)量直接影響后續(xù)淺剖層界自動提取的效果.

2)根據(jù)淺剖層界位置處，Ping灰度值突變原理提取層界，實現(xiàn)了淺剖層界提取自動化．但該方法需要一個準確的顯著性系數(shù)k來界定灰度突變點，以期實現(xiàn)淺剖層界的準確提取.

3)采用二維圖像信息熵作為約束，即根據(jù)不同區(qū)塊的不同信息熵值設置顯著性系數(shù)k，提高了淺剖層界提取的準確性，減少了圖像細節(jié)信息的損失，避免了強噪聲或偽層界造成的干擾，實現(xiàn)了智能化、自動化提取真實層界的目的．

[1] 李斌,楊文達,張異彪,等. 海底管道的淺地層剖面圖上反射特征與判讀方法[J]. 海洋測繪, 2010,30(5):56-58. DOI:10.3969/j.issn.1671-3044.2010.05.016.

LI Bin, YANG Wenda, ZHANG Yibiao, et al. Characteristics of the seabed pipeline’s reflection in the sub-bottom profiling and interpretation[J]. Hydrographic Surveying and Charting, 2010, 30(5):56-58. DOI:10.3969/j.issn.1671-3044.2010.05.016.

[2] HONSHO C, URA T, ASADA A, et al. High-resolution acoustic mapping to understand the ore deposit in the Bayonnaise knoll caldera, Izu-Ogasawara arc[J]. Journal of Geophysical Research: Solid Earth, 2015, 120(4):2070-2092. DOI: 10.1002/2014JB011569.

[3] 王方旗. 淺地層剖面儀的應用及資料解譯研究[D]. 青島：國家海洋局第一海洋研究所,2010.

WANG Fangqi. The research of the application and data interpretation of sub-bottom profiler[D]. Qingdao: The First Institute of Oceanography, SOA, 2010.

[4] MACLEAN B, BLASCO S, BENNETT R, et al. New marine evidence for a Late Wisconsinan ice stream in Amundsen Gulf, Arctic Canada[J]. Quaternary Science Reviews, 2015, 114: 149-166. DOI: 10.1016/j.quascirev.2015.02.003.

[5] PLETS R M K, DIX J K, ADAMS J R, et al. The use of a high-resolution 3D Chirp sub-bottom profiler for the reconstruction of the shallow water archaeological site of the Grace Dieu (1439), River Hamble, UK[J]. Journal of Archaeological Science, 2009, 36(2):408-418.DOI: 10.1016/j.jas.2008.09.026.

[6] CARLIN J A, DELLAPENNA T M, STROM K, et al. The influence of a salt wedge intrusion on fluvial suspended sediment and the implications for sediment transport to the adjacent coastal ocean: a study of the lower Brazos River TX, USA[J]. Marine Geology, 2015, 359: 134-147.DOI: 10.1016/j.margeo.2014.11.001.

[7] BAETEN N J, LABERG J S, VANNESTE M, et al. Origin of shallow submarine mass movements and their glide planes-Sedimentological and geotechnical analyses from the continental slope off northern Norway[J]. Journal of Geophysical Research: Earth Surface, 2014,119(11): 2335-2360.DOI: 10.1002/2013JF003068.

[8] HAO Yanjun, MCINTOSH K, MAGNANI MB. Long-lived deformation in the southern Mississippi Embayment revealed by high-resolution seismic reflection and sub-bottom profiler data[J]. Tectonis,2015,34(3): 555-570. DOI: 10.1002/2014TC003750.

[9] 張恒超. 幾種簡單實用的多次波去除方法及其應用[J]. 內(nèi)蒙古石油化工，2011(21)：21-24．DOI: 10.3969/j.issn.1006-7981.2011.21.010.

ZHANG Hengchao. Several practical methods of multiple waves suppression and their applications[J]. Inner mongulia petrochemical industry, 2011(21)：21-24．DOI: 10.3969/j.issn.1006-7981.2011.21.010.

[10]華長發(fā),范建平,高傳善,等. 基于二維熵閾值的圖像分割及其快速算法[J]. 模式識別與人工智能, 2000,13(1)：42-46. DOI: 10.3969/j.issn.1003-6059.2000.01.009.

HUA Changfa，F(xiàn)AN Jianpin，GAO Chuanshan, et al. Image segmentation based on 2d entropic thresholding and its fast algorithm[J]. Pattern Recognition and Artificial Intelligence, 2000,13(1)：42-46．DOI: 10.3969/j.issn.1003-6059.2000.01.009.

[11]吳澤鵬,郭玲玲,朱明超,等. 結(jié)合圖像信息熵和特征點的圖像配準方法[J]. 紅外與激光工程, 2013,42(10)：2846-2852. DOI: 10.3969/j.issn.1007-2276.2013.10.046.

WU Zepeng, GUO Lingling, ZHU Mingchao, et al. Improved image registration using feature points combined with image entropy[J]. Infrared and Laser Engineering, 2013,42(10)：2846-2852. DOI: 10.3969/j.issn.1007-2276.2013.10.046.

[12]GABARDA S, CRISTOBAL G. The generalized Renyi image entropy as a noise indicator-art[C]// Proceedings of the SPIE 6603, Noise and Fluctuations in Photonics, Quantum Optics, and Communications. Florence: SPIE,2007. DOI: 10.1117/12.725086.

[13]王增波,李雁鴻,趙劍,等. SEG-Y地震數(shù)據(jù)格式解析及轉(zhuǎn)換方法[J]. 物探裝備, 2012,22(3):177-182.DOI: 10.3969/j.issn.1671-0657.2012.03.010.

WANG Zengbo，LI Yanhong，ZHAO Jian, et al. Analytical method and conversion method for SEG-Y data[J]. Equipment for Geophysical Prospecting, 2012, 22(3): 177-182. DOI: 10.3969/j.issn.1671-0657.2012.03.010.

[14]肖梅,劉國華,李慶春. 在VC++環(huán)境下讀取地震勘探SEG-Y格式數(shù)據(jù)及其應用[J]. 中國科技信息, 2005,(9):17-17, 29. DOI: 10.3969/j.issn.1001-8972.2005.09.012.

[15]李枝文,宋偉,肖柏勛,等. Hilbert變換與小波變換在探地雷達資料處理中的應用[J]. 工程地球物理學報, 2012,9(4):428-432.DOI: 10.3969/j.issn.1672-7940.2012.04.011.

LI Zhiwen, SONG Wei, XIAO Boxun, et al. Application of hilbert transform and wavelet transform to data processing of ground penetrating radar[J]. Chinese Journal of Engineering Geophysics, 2012,9(4):428-432.DOI: 10.3969/j.issn.1672-7940.2012.04.011.

[16]張軍華,繆彥舒,鄭旭剛,等. 預測反褶積去多次波幾個理論問題探討[J]. 物探化探計算技術,2009,31(1)：6-10. DOI: 10.3969/j.issn.1001-1749.2009.01.003.

ZHANG Junhua, MIAO Yanshu, ZHENG Xugang, et al. Discussion of several theoretical questions to remove seismic multiples using predictive deconvolution[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2009,31(1)：6-10. DOI: 10.3969/j.issn.1001-1749.2009.01.003.

(編輯 張 紅)

Demarcation of the sub-bottom layers based on image information entropy constraint

ZHAO Jianhu1, FENG Jie1, SHI Feng1, ZHANG Hongmei2, HE Linbang1

(1.School of Geodesy and Geomatics, Wuhan University,Wuhan 430079, China; 2.Shool of Power and Mechanical Engineering, Wuhan University, Wuhan 430072, China)

To address the issue of the sub-bottom profile layer extraction in complex circumstance, this paper proposes a new demarcating method based on constraint of image information entropy. Firstly, the image of sub-bottom is divided into different blocks; then, the information entropy in each block is calculated and a relation model of information entropy and significant parameters are established according to drilling data; finally, the whole sub-bottom profiling is demarcated according to the model. It is revealed that this method has overcome the shortcomings of existing methods, realized the self-adapting and exacted demarcation of sub-bottom layers. The experiment has gained the same accuracy as the depth and thickness of layers got by drilling data.

sub-bottom profiling; sub-bottom layer and its extraction; 2-D entropy; coefficient of mean square error; self-adaption extraction

10.11918/j.issn.0367-6234.201504135

2015-04-30

國家自然科學基金(41376109，41176068，41576107)

趙建虎(1970—)，男，教授，博士生導師

馮 杰，1017587956@qq.com

P229

A

0367-6234(2017)08-0165-06