CN112529841A - 多波束水柱数据中海底气体羽状流处理方法、系统及应用 - Google Patents

多波束水柱数据中海底气体羽状流处理方法、系统及应用 Download PDF

Info

Publication number
CN112529841A
CN112529841A CN202011278968.2A CN202011278968A CN112529841A CN 112529841 A CN112529841 A CN 112529841A CN 202011278968 A CN202011278968 A CN 202011278968A CN 112529841 A CN112529841 A CN 112529841A
Authority
CN
China
Prior art keywords
data
point
wci
water
ping
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011278968.2A
Other languages
English (en)
Other versions
CN112529841B (zh
Inventor
丁咚
秦浩森
李广雪
任昕
王蕊
韩慧慧
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ocean University of China
Original Assignee
Ocean University of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ocean University of China filed Critical Ocean University of China
Priority to CN202011278968.2A priority Critical patent/CN112529841B/zh
Publication of CN112529841A publication Critical patent/CN112529841A/zh
Application granted granted Critical
Publication of CN112529841B publication Critical patent/CN112529841B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/20—Image enhancement or restoration using local operators
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/70—Determining position or orientation of objects or cameras
    • G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/28—Quantising the image, e.g. histogram thresholding for discrimination between background and foreground patterns
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20024—Filtering details
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • Multimedia (AREA)
  • Quality & Reliability (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明属于海洋测绘、海洋探测技术领域，公开了一种多波束水柱数据中海底气体羽状流提取方法、系统及应用，气体或气‑液混合的水体会比正常的水体环境产生更大的声阻抗，为利用水体的散射信息来识别气体羽状流提供了可能性。针对多波束原始数据进行解析提取，使用提取到的船位坐标及艏向、测量船姿态信息等信息计算获取各类多波束WCI并依据其空间参考进行三维重建，构建基于水体声学散射信息的三维点云模型。利用图像学方法实现WCI降噪、滤波、二值化以从整体水体点云中初步分离羽状流部分数据；将Kmeans、局部异常因子数据挖掘算法引入水体信息处理，在尽可能保证羽状流保留的基础上进一步去噪，提升羽状流的提取效果。

Description

多波束水柱数据中海底气体羽状流处理方法、系统及应用
技术领域
本发明属于海洋测绘、海洋探测技术领域，尤其涉及一种多波束水柱(WaterColumn，WC)数据中海底气体羽状流处理方法、系统及应用。
背景技术
目前，海底气体渗漏是在全球海洋环境中广泛分布的一种自然现象，渗漏出的气体(如CH4、CO2等) 可能影响全球的气候变化。同时，海底气体渗漏也是天然气水合物成藏区识别的重要标志之一。海底气体的渗漏现象，在全球多个海域，尤其是在陆架边缘区普遍存在，往往伴随着一些特殊的地质活动现象，如热液、冷泉的分布、天然气水合物成藏、海底火山的发育等。海底气体渗漏区的圈定和海底气体羽状流的快速识别，有益于高效率的获取海底气体渗漏点分布，对气体渗漏区地形地貌和底质分布研究、渗漏气体通量估算、资源勘查以及气候环境变化研究具有重要意义；此外，基于多波束探测的水柱数据也可以用于海底输气管道的渗漏检测。多波束测深系统是在单波束回声测深仪的基础上发展起来的，与单波束测深仪相比，有测量范围大、测量速度快、精度和效率高的特点。随着现代声呐技术的发展，多波束水体信息记录的巨大磁盘空间需求和数字化的水体信息记录问题得以解决，新一代的多波束测深系统除了可以通过声波确定水深外，也可以记录声波传播过程中途经的水体和海底的声散射信息。
对于海底气体渗漏区的探测可以采用多种调查技术手段实现，如卫星遥感等针对海面的探测，其探测效率高、范围大，对于浅水区或海水表层的渗漏表征现象探测具有优势，但由于电磁波穿透性弱，无法探测水深较大的海底渗漏现象；多道地震声反射数据成像可以结合地层分布信息进行海底渗漏的识别，但由于其为保证声波穿透性，发射声波波长大、分辨率低、覆盖面窄，无法完全反馈水体环境信息，难以探测气体渗漏点；CTD观测及水体取样、ROV探测及深海摄像、海底观测站原位监测等探测手段对工作量、观测成本高，事先需要了解测区内气体渗漏概况，在前期准备上都有较高的要求。
通过上述分析，现有技术存在的问题及缺陷为：
(1)现有技术对于浅水区或海水表层的渗漏表征现象探测具有优势，但无法探测水深较大的海底渗漏现象；
(2)现有技术多道地震声反射数据成像可以结合地层分布信息进行海底渗漏的识别，但难以准确探测气体渗漏点位置；
(3)现有技术CTD观测及水体取样、ROV探测及深海摄像、海底观测站原位监测等探测手段对工作量、观测成本高，事先需要了解测区内气体渗漏概况，在前期准备上有较高的要求且探测成功率低。
解决以上问题及缺陷的难度为：
(1)多波束具有探测效率高、分辨率高、覆盖宽的优点，且可有效记录水体中羽状流目标，但从多波束原始数据中解析、计算获取测区完整水体点云信息、多波束WCI波束阵列图、扇面图(相对坐标和投影坐标)、航向剖面图为难点之一。
(2)多波束WC数据记录了异常目标，但其中存在成因广、范围大的噪声给提取气体羽状流带来极大阻碍。如何在尽可能高的自动化程度基础上同时将不同来源、分布的噪声去除，将羽状流剥离出来为难点之一。
(3)在去噪提取羽状流的过程中，选用何种方法以保证羽状流提取完整性，并良好去除噪声为难点之一。
解决以上问题及缺陷的意义为：
(1)给出了从多波束原始数据中解析计算获取测区完整水体点云信息、多波束WCI波束阵列图、扇面图(相对坐标和投影坐标)、航向剖面图的方法。
(2)采用图像学方法初步降噪分离羽状流水体点云，数据挖掘进一步去噪提取羽状流的方法，参数可变动范围大，自动化程度较高。
(3)可快速获取气体渗漏点位置，了解测区羽状流喷溢状态，提高探测效率，削减探测成本。
发明内容
针对现有技术存在的问题，本发明提供了一种多波束水柱数据中海底气体羽状流处理方法、系统及应用。
本发明是这样实现的，一种多波束水柱数据中海底气体羽状流处理方法，所述多波束水柱数据中海底气体羽状流处理方法包括：
解析提取多波束原始数据信息，根据数据记录信息完成水体、水深采样点归位，计算回波强度值真实的地理位置，获取测区完整水体点云信息、水深地形图、各类多波束WCI；通过对称作差、图像滤波实现多波束水柱影像提取气体羽状流前处理工作；若WCI无第三方声呐干扰，图像二值化后即可获取初步分离气体羽状流的水体点云；若WCI有第三方声呐干扰，则先选取N Ping不包含羽状流的二值化WCI，采用开操作设置像素阈值滤除小连通分量的方式将第三方声呐干扰提取出来，并对N Ping提取结果求并生成模板，滤除全部二值化WCI每Ping中与模板相交部分，以获取初步分离气体羽状流后的水体点云；引入 Kmeans、局部异常因子方法进一步迭代去噪提取气体羽状流体；绘制各类多波束WCI、气体羽状流点云图、多波束水深地形图。
进一步，所述多波束水柱数据中海底气体羽状流处理方法还包括：
第一步，根据数据存储方式解析提取原始数据，转换为可用格式数据；
第二步，转换后的数据显示多波束记录了水深、底质、水体回波强度、采样频率、声速及波束入射角信息，利用提取到的信息，对每一Ping数据进行采样点归位计算，得到每一水体、水深采样点实际位置信息，测区完整水体点云信息、多波束WCI波束阵列图、扇面图(相对坐标和投影坐标)、航向剖面图；
第三步，船舶、换能器自身产生的噪声仅在水体浅表层数十次采样序列中存在，直接将浅表层水体一定数量等时采样序列去除，并以中央波束为分界线，将两侧WC数据对称作差初步压制WCI中噪声；
第四步，多波束WCI滤波：对滤除表层噪声并对称作差法初步去噪的WCI采用图像滤波方法，平滑增强气体羽状流目标，完成数据前处理；
第五步，基于大津算法(OSTU)逐一对全部Ping多波束WCI进行二值化处理，并设置系数w供用户设置，以改善二值化结果；
第六步，若多波束WCI有第三方声呐干扰，则先选取N Ping不包含羽状流的二值化WCI，采用开操作设置像素阈值P滤除小连通分量的方式提取第三方声呐干扰，并对N Ping提取结果求并生成模板，滤除全部二值化WCI每Ping中与模板相交部分，以获取初步分离气体羽状流后的水体点云；若无则第五步处理结果即为初步分离气体羽状流的水体点云；
第七步，将每Ping初步分离气体羽状流的水体点云在以换能器为原点的相对坐标空间叠加，首次尝试以点云水深值为特征值，通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤ r则不去除；
第八步，采用LOF方法迭代去除离散分布噪点；
第九步，再次尝试以点云水深值为特征值，通过Kmeans将点云分为2类，若类间聚类中心比值≥r 去除上层分布噪点，≤r则不去除，完成羽状流提取；
第十步，利用提取结果和解析的数据，实现气体羽状流点云图、多波束各类WCI、地形图的成图。
进一步，所述第一步原始数据的解析方法包括：根据文件格式，编写程序，解析提取原始数据，保存为可用的格式；
所述第二步的水深及水体采样点计算包括：
根据原始数据中水深地形及POS数据完成水深数据投影坐标计算，公式：
Figure RE-GDA0002867399030000031
[THETA，RHO]＝cart2pol(AcrosstrackY，AlongtrackX)
Figure RE-GDA0002867399030000032
S_East＝S_Eastsonar+P_prcEast
S_North＝S_Northsonar+P_prcNorth
式中：P_prc为每Ping数据获取时投影坐标(P_prcEast、P_prcNorth)、航向(P_prcHead)、子午线收敛角(P_prcCon)；posdown、posup为与每Ping水深数据记录时间最邻近的两组POS点；AcrosstrackY， AlongtrackX为以换能器为原点采样点横向、纵向相对坐标，m；S_Eastsonar，S_Northsonar为水深采样点到换能器距离，m。编程计算各地形采样点投影位置即可获取水深地形图；
结合导航数据包、水柱数据包内存储数据，计算水体后向散射强度采样点以换能器为原点的相对坐标，公式：
PBS_SR＝(PBS_idx+PBS_sRSNum)×c/(2·f)
PBS_H＝-PBS_SR×cos(PBS_BPA)
PBS_X＝-PBS_SR×sin(PBS_BPA)
获取水体后向散射强度采样点相对坐标后，结合每Ping数据换能器投影坐标，计算获取水体后向散射强度采样点投影坐标，公式：
Figure RE-GDA0002867399030000041
PBS_East＝P_prcEast+PBS_X×cos(theta)
PBS_North＝P_prcNorth+PBS_X×sin(theta)
式中：PBS_SR为水体采样点到换能器距离，m；BS_idx为每柱接收波束回波采样点编号；PBS_sRSNum 为每柱接收波束起始采样点编号；c为声速，m；f为采频率；BS_BPA为波束发射角；PBS_X，PBS_H为水体后向散射强度采样点沿船横向位置及水深，m；sonarHeadingOffset为换能器安装角度，°；theta为波束方位角，°；PBS_East、PBS_North为水体采样点投影坐标，m。编程计算各水体采样点相对坐标和投影坐标后，获取测区完整水体点云信息、多波束WCI波束阵列图、扇面图(相对坐标和投影坐标)、航向剖面图。
进一步，所述第三步噪声初步压制包括：
针对测区内噪声以中央波束为轴对称分布特征，对多波束WC数据单Ping内呈镜面对称位置采样点的水体后向反射强度进行Pearson相关性分析，并求和取均值可见单Ping内镜面对称位置水体后向散射强度采样点呈中等程度正相关，可以中央波束为分界线，将两侧WC数据对称作差初步压制WCI中噪声；
进一步，所述第四步的WCI滤波包括：
通过图像滤波方法对多波束WCI数据进行滤波，平滑增强气体羽状流目标，完成数据前处理；
进一步，所述第五步获取多波束WCI二值化包括：
以中央波束为轴将每Ping WCI图像分成两侧，采用OSTU算法初步确定两侧二值化阈值，并设置系数w作为此阈值的系数，以改善二值化结果，公式：
Figure RE-GDA0002867399030000042
式中：mG为WCI平均灰度；概率P1(k)为灰度阈值为k时被分到小于k一类的概率；mk为灰度级k的累加均值；δB 2为类间方差，将每Ping单侧多波束WCI中包含的灰度级依次代入上式，当δB 2取最大时的 k*即为此算法下计算获取阈值，再乘以预先设置的w，即可得到二值化阈值。分别计算每Ping WCI两侧阈值完成二值化后拼接，即可完成多波束WCI二值化。
进一步，所述第六步多波束WCI第三方干扰去除包括：
若多波束WCI有第三方声呐干扰，则选取N Ping不包含羽状流的二值化WCI图像执行开操作，即对图像腐蚀，然后膨胀的处理方式，图像腐蚀、膨胀公式分别为：
Figure RE-GDA0002867399030000043
Figure RE-GDA0002867399030000044
式中：Z为某Ping WCI；A为Z中包含的连通分量；E为腐蚀后Z中包含的连通分量；SE为结构体；开操作处理后WCI中指示第三方声呐干扰的连通分区像素数量较多，设置合理阈值P滤除像素数量小的连通分区。然后对N Ping处理结果求并完成第三方声呐干扰模板提取，逐一滤除全部Ping二值化WCI中与模板相交部分，以获取初步分离气体羽状流后的水体点云，若无第三方声呐干扰，则通过多波束WCI 二值化结果即可获取初步分离气体羽状流的水体点云。
进一步，所述第七步和第九步Kmeans去噪包括：
在LOF迭代前后分别尝试以点云水深值为特征值，采用无监督聚类方法K-means将点云分成2类，设置最上层聚类中心与下层簇中心比值大于r为限。当比值大于r时，去除聚类后水体上层由噪声干扰产生的数据点；
进一步所述第八步LOF迭代去噪包括：
将每Ping初步分离气体羽状流的水体点云在以换能器为原点的扇形空间坐标(Depth-Acrosstrack空间， D-A空间)下叠加，剩余下层噪点是在获取初步分离气体羽状流的水体点云时没有完全消除的旁瓣、环境噪声、第三方声呐等类型噪点，往往为离散分布噪点。利用下层噪点离散且叠加在D-A空间后，单位区域数量远小于气体羽状流数据点特征，通过从点云中随机选取一定数量点的方法使噪点在抽稀后的D-A空间中表现为离群点，再采用LOF迭代方法去除该类噪点。所述LOF迭代去噪具体方法如下：
①点云数据输入，将每Ping OSTU二值化(或OSTU二值化并去除第三方声呐干扰)后获取的初步分离羽状流点云叠加到Depth-Acrosstrack空间，构成初始输入数据；
②设置参数，设置参数邻域k、LOF阈值thresh、迭代次数iter、随机抽取样本点数量num；
③随机抽取样本点集合O，从点云数据中随机抽取num个样本点，形成样本点集合O＝{o1，o2， o3，…on}，其中oi是集合O中第i个样本点；
④计算第k距离，计算集合O每个样本点oi的第k距离distk(oi)，即距oi第k远的样本点的距离；
⑤计算k邻域内样本点真实距离，计算集合O每个样本点oi的k邻域内(Nk(oi))点oi′到点oi的距离 dist(oi，oi′)，即点oi第k距离以内的所有点，包括第k距离在内的点到oi的欧式距离；
⑥计算可达距离，计算集合O每个样本点oi的k邻域内(Nk(oi))点oi′到点oi的第k可达距离，其定义如下：
reachdistk(oi←oi′)＝max{distk(oi′)，dist(oi，oi′)}；
即点oi′到点oi的第k可达距离，看出其值至少是oi′的第k距离，或者为oi、oi′间的真实距离；
⑦计算局部可达密度，计算集合O中每个样本点oi的局部可达密度lrd，其定义如下：
Figure RE-GDA0002867399030000051
即点oi的第k邻域内点oi′到oi平均可达距离的倒数。其意在量化oi与其k邻域内点之间密度差异，若lrdk越高则越可能属于同一簇点，越低则越可能为离群点；
⑧计算LOF，根据计算所得lrd，计算集合O中每个样本点oi的局部离群因子LOF，其定义如下：
Figure RE-GDA0002867399030000052
表示点oi的邻域点Nk(oi)的局部可达密度与点oi的局部可达密度之比的平均数，若LOF值越接近于1，两者的密度越接近；若LOF值小于1，则oi密度高于其邻域点；若LOF值大于1，则值越大说明点oi密度小于邻域点，越可能是离群点，即更可能为噪点；
⑨滤除LOF异常高点，根据预先设置好的LOF阈值thresh，滤除大于thresh的点云数据，若迭代次数小于iter则转至③，否则顺序执行；
⑩有禁忌搜索随机抽取计算
Figure RE-GDA0002867399030000061
滤除空间内未去除噪点，在首次LOF迭代去噪完成后，采用有禁忌搜索随机选取集合O对全部点云重新完整迭代计算LOF滤除首次迭代后依旧未完全去除的噪点，完善提取结果后结束。
本发明的另一目的在于提供一种实施所述多波束水柱数据中海底气体羽状流处理方法的多波束水柱数据中海底气体羽状流处理系统，所述多波束水柱数据中海底气体羽状流处理系统包括：
数据解析模块，用于不同公司不同型号不同存储模式多波束原始数据的解析提取，将其转换为可用格式数据；
水体及水深采样点归位模块，用于根据原始数据提取到的多波束水深、底质、水体回波强度、采样频率、声速及波束入射角等信息，对每一Ping数据进行采样点归位计算，得到每一水体、水深采样点实际位置信息，进而获取测区水深地形图、水体点云信息、多波束WCI波束阵列图、扇面图(投影和相对坐标)、航向剖面图；
噪声初步压制模块，用于去除船舶、换能器自身在水体浅表层数十次采样序列中产生的噪声；以中央波束为分界线，将两侧WC数据对称作差初步压制WCI中噪声；
多波束WCI滤波模块，用于对滤除表层噪声并对称作差法初步去噪的WCI采用图像滤波方法，平滑增强气体羽状流目标，完成数据前处理；
第三方声呐干扰去除模块，用于多波束WCI有第三方声呐干扰时，先选取N Ping不包含羽状流的二值化WCI，采用开操作设置像素阈值P滤除小连通分量的方式提取第三方声呐干扰，并对N Ping提取结果求并生成模板。逐一滤除全部Ping中与模板相交部分，以获取初步分离气体羽状流后的水体点云；若无则多波束WCI二值化处理后即可获取初步分离气体羽状流的水体点云；
数据挖掘去噪模块：用于去除初步分离气体羽状流水体点云中的剩余噪点，LOF迭代去噪前，以点云水深值为特征值，尝试通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤r不去除。然后，采用LOF方法迭代去除剩余离散分布噪点，再次尝试通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤r不去除，完成羽状流提取；
绘图模块，用于利用提取结果和解析的数据，实现气体羽状流点云图、各类多波束WCI、地形图的成图。
本发明的另一目的在于提供一种用于海底输气管道的渗漏检测方法，所述用于海底输气管道的渗漏检测方法包括以下步骤：
解析提取多波束原始数据信息，根据数据记录信息完成水体、水深采样点归位，计算回波强度值真实的地理位置，获取测区完整水体点云信息、水深地形图、各类多波束WCI；通过对称作差、图像滤波实现多波束水柱影像提取油气渗漏前处理工作；若WCI无第三方声呐干扰，图像二值化后即可获取初步分离油气渗漏的水体点云；若WCI有第三方声呐干扰，则先选取N Ping不包含油气渗漏的二值化WCI，采用开操作设置像素阈值滤除小连通分量的方式将第三方声呐干扰提取出来，并对N Ping提取结果求并生成模板，逐一滤除全部Ping二值化WCI中与模板相交部分，以获取初步分离油气渗漏后的水体点云；引入 Kmeans、局部异常因子方法进一步迭代去噪提取油气渗漏点；绘制各类多波束WCI、油气渗漏点云图、多波束水深地形图。
结合上述的所有技术方案，本发明所具备的优点及积极效果为：气体或气-液混合的水体会比正常的水体环境产生更大的声阻抗，这为利用水体的散射信息来识别气体羽状流提供了可能性。针对多波束原始数据进行解析提取，使用提取到的船位坐标及艏向、测量船姿态信息、声速输入信息、激发Ping编号及时间、波束编号、水体散射强度、采样时间、采样距离、波束角等信息，计算获取各类多波束WCI并依据其空间参考进行三维重建，构建基于水体声学散射信息的三维点云模型。利用图像学方法实现WCI降噪、滤波、二值化以从整体水体点云中初步分离羽状流部分数据点；将Kmeans、局部异常因子数据挖掘算法引入水体信息处理，在尽可能保留羽状流的基础上进一步去噪，提升羽状流的提取效果，进而提出了一种基于水体点云数据挖掘提取多波束声呐水柱数据中海底气体羽状流的处理方法。
本发明的目的在于提供从多波束WCI中提取气体羽状流的方法，本发明首次对通过利用多波束原始数据中获取的水体点云，结合图像学方法实现WCI降噪、滤波、二值化以从整体水体点云中初步分离羽状流部分数据，并采用数据挖掘方法进一步去噪完成气体羽状流，且与噪声点良好分离，实现羽状流高效提取，提高实际探测效率。此方法可用于海底气体渗漏探测及海底管线油气渗漏检测。
气体或气-液混合的水体会比正常的水体环境产生更大的声阻抗，这为利用水体的散射信息来识别气体羽状流提供了可能性。针对多波束原始数据进行解析提取，使用提取到的船位坐标及艏向、测量船姿态信息、声速输入信息、激发Ping编号及时间、波束编号、水体散射强度、采样时间、采样距离、波束角等信息计算获取各类多波束WCI并依据其空间参考进行三维重建，构建基于水体声学散射信息的三维点云模型。利用图像学方法实现WCI降噪、滤波、二值化以从整体水体点云中初步分离羽状流部分数据；将 Kmeans、局部异常因子(Local Outlier Factor，LOF)数据挖掘算法引入水体信息处理，在尽可能保证羽状流保留的基础上进一步去噪，提升羽状流的提取效果，进而提出了一种基于水体点云数据挖掘提取多波束声呐水柱数据中海底气体羽状流处理方法。本发明有益于高效获取海底气体渗漏点位置，对气体渗漏区地形地貌和底质分布研究、渗漏气体通量估算、资源勘查以及气候环境变化研究、海底输气管道渗漏检测以及水体中生物或障碍物或其他悬浮物体探测与探测等方面具有重要意义。
附图说明
为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。
图1是本发明实施例提供的多波束水柱数据中海底气体羽状流处理方法流程图。
图2是本发明实施例提供的多波束水柱数据中海底气体羽状流处理系统的结构示意图；
图2：1、数据解析模块；2、水体及水深采样点归位模块；3、噪声初步压制模块；4、多波束WCI 滤波模块；5、多波束WCI二值化模块；6、多波束WCI第三方干扰去除模块；7、数据挖掘去噪模块；8、绘图模块。
图3是本发明实施例提供的测区原始数据航向剖面图及示例WCI图像
图4本发明实施例提供的相关系数分析结果示意图。
图5是本发明实施例提供的多波束水深地形图。
图6是本发明实施例提供的各类多波束WCI。
图7是本发明实施例提供的对称作差法降噪对比图。
图8是本发明实施例提供的巴特沃斯低通滤波效果图。
图9是本发明实施例提供的多波束WCI二值化效果图。
图10是本发明实施例提供的第三方声呐干扰模板。
图11是本发明实施例提供的初步分离气体羽状流点云图。
图12是本发明实施例提供的初步分离羽状流数据点D-A空间叠加图。
图13是本发明实施例提供的Kmeans上层噪点去除效果图。
图14是本发明实施例提供的LOF迭代去噪流程图。
图15是本发明实施例提供的LOF迭代过程中某次LOF计算结果(圆圈大小代表LOF值)。
图16本发明实施例提供的最终气体羽状流提取结果。
图17是本发明验证例提供的测区原始数据航向剖面图及示例WCI图像。
图18是本发明验证例提供的多波束水深地形图。
图19是本发明验证例提供的相关系数分析结果示意图。
图20是本发明验证例提供的对称作差法降噪对比图。
图21是本发明验证例提供的巴特沃斯低通滤波效果图。
图22是本发明验证例提供的多波束WCI二值化效果图。
图23是本发明验证例提供的初步分离气体羽状流点云图。
图24是本发明验证例提供的Kmeans上层噪点去除效果图。
图25是本发明验证例提供的初步分离羽状流点数据点D-A空间叠加图。
图26是本发明验证例提供的LOF迭代过程中某次LOF计算结果(圆圈大小代表LOF值)。
图27是本发明验证例提供的最终气体羽状流提取结果。
图28是本发明提供的处理流程框架图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例(图3)，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
针对现有技术存在的问题，本发明提供了一种多波束水柱数据中海底气体羽状流处理方法、系统及应用，下面结合附图对本发明作详细的描述。
如图1所示，本发明提供的多波束水柱数据中海底气体羽状流处理方法包括以下步骤：
S101：根据数据存储方式解析提取原始数据，转换为可用格式数据；
S102：转换后的数据显示多波束记录了水深、底质和水体回波强度、采样频率、声速及波束入射角等信息，利用这些提取到的信息，对每一Ping数据进行采样点归位计算，得到每一水体、水深采样点实际位置信息，进而获取测区水深地形图、完整水体点云信息、各类多波束WCI；
S103：去除船舶、换能器自身在水体浅表层数十次采样序列中产生的噪声，并以中央波束为分界线，将两侧呈镜面对称WC数据对称作差初步压制WCI中噪声；
S104：对滤除表层噪声并对称作差初步去噪的WCI采用图像滤波方法，平滑增强气体羽状流目标，完成数据前处理；
S105：基于大津算法(OSTU)逐一对全部Ping多波束WCI进行二值化处理，并设置系数w供用户设置，以改善二值化结果；以改善二值化结果；
S106：若多波束WCI有第三方声呐干扰，则先选取N Ping不包含羽状流的二值化WCI，采用开操作设置像素阈值P滤除小连通分量的方式提取第三方声呐干扰，并对N Ping提取结果求并生成模板，逐一滤除全部Ping二值化WCI中与模板相交部分，以获取初步分离气体羽状流后的水体点云；若无则第五步处理结果即可获取初步分离气体羽状流的水体点云；
S107:将每Ping初步分离气体羽状流的数据点在以换能器为原点的相对坐标空间叠加，首次尝试以点云水深值为特征值，通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤r则不去除；
S108：采用LOF方法迭代去除剩余离散分布噪点；
S109：再次尝试以点云水深值为特征值，通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤r则不去除，完成羽状流提取；
S110：利用提取结果和解析的数据，可实现气体羽状流点云图、多波束各类WCI、地形图的成图。
本发明提供的多波束水柱数据中海底气体羽状流处理方法业内的普通技术人员还可以采用其他的步骤实施，图1的本发明提供的多波束水柱数据中海底气体羽状流处理方法仅仅是一个具体实施例而已。
如图2所示，本发明提供的多波束水柱数据中海底气体羽状流处理系统包括：
数据解析模块1，用于不同公司不同型号不同存储模式多波束原始数据的解析提取，将其转换为可用格式数据；
水体及水深采样点归位模块2，用于根据原始数据提取到的水深、底质、水体回波强度、采样频率、声速及波束入射角等信息，对每一Ping数据进行采样点归位计算，得到每一水体、水深采样点实际位置信息，进而获取测区完整水体点云信息、水深地信息、各类多波束WCI；
噪声初步压制模块3，用于去除船舶、换能器自身在水体浅表层数十次采样序列中产生的噪声，并以中央波束为分界线，将两侧呈镜面对称WC数据对称作差初步压制WCI中噪声；
多波束WCI滤波模块4，用于对滤除表层噪声并对称作差法初步去噪的WCI采用图像滤波方法，平滑增强气体羽状流目标，完成数据前处理；；
多波束WCI二值化模块5，基于大津算法(OSTU)逐一对全部Ping多波束WCI进行二值化处理，并设置系数w供用户设置，以改善二值化结果；
多波束WCI第三方干扰去除模块6：用于多波束WCI有第三方声呐干扰时，先选取NPing不包含羽状流的二值化WCI，采用开操作设置像素阈值P滤除小连通分量的方式提取第三方声呐干扰，并对N Ping 提取结果求并生成模板。逐一滤除全部二值化WCI中与模板相交部分，以获取初步分离气体羽状流后的水体点云；若无则多波束WCI二值化处理后即可获取初步分离气体羽状流的水体点云；
数据挖掘去噪模块7：用于去除水体点云数据中的剩余噪点，LOF迭代去噪前后，以点云水深值为特征值，尝试通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤r不去除。然后，采用LOF方法迭代去除剩余离散分布噪点，完成羽状流提取；
绘图模块8，用于利用提取结果和解析的数据，可实现气体羽状流点云图、多波束各类WCI、地形图的成图。
下面结合附图对本发明的技术方案作进一步的描述。
本发明多波束测深仪采用发射、接收指向性正交的两组换能器阵向海底发射声波信号并接收回波信号，通过接收指向性与发射指向性叠加，获得一系列垂直于航向分布的窄波束。而在通过束控实现换能器基阵指向性的同时，除主瓣外还会产生一系列旁瓣，对换能器到海底最短距离半径即最小倾斜距离(MSR) 内WCI产生严重干扰。因此，MSR内的WCI受旁瓣干扰小质量高，MSR以外WCI受旁瓣干扰大，一般难以反馈水体真实信息，相关学者选取MSR以内WCI分析提取水体中异常目标，本方法也采用同样选择。此外，多波束系统将海洋风浪、浮游动植物、悬浮物等产生的海洋环境噪声，船舶、第三方声呐等产生的电子噪声也以异常值的形式加以记录，对多波束WCI造成污染，给利用WCI目标探测提取产生不利影响。因此，本文通过对多波束水柱影像中存在的噪声及气体羽状流研究发现：以中央波束为轴镜面对称分布采样点的水体后向散射强度呈正相关(图4)，且羽状流反射强度高于对称位置水体声强特点，可利用这一特征实现多波束WCI初步降噪；WCI降噪、滤波、二值化初步分离羽状流(若WCI有第三方声呐干扰，则还需滤除此干扰)的水体点云，其下层噪点为离散分布噪点，上层多为广泛分的大范围噪点，且与下层羽状流数据分层明显，可结合Kmeans和局部异常因子数据挖掘方法进一步去除。基于上述特征，提出一种基于水体点云数据挖掘提取多波束声呐水柱数据中海底气体羽状流处理方法，具体包括以下步骤：
(1)数据解析：不同公司不同型号多波束原始数据都有其相应的存储模式，需根据其数据存储方式解析提取原始数据，转换为可用格式数据；
(2)水体及水深采样点归位：转换后的数据显示多波束记录了水深、底质和水体回波强度、采样频率、声速及波束入射角等信息，利用这些提取到的信息，对每一Ping数据进行采样点归位计算，得到每一水体、水深采样点实际位置信息，进而获取测区水深地形图(图5)、完整水体点云图、航向剖面图、多波束WCI波束阵列图(Time-Angle空间，T-A空间，本方法处理过程中所选取的WCI，且由于换能器到海底最短距离半径外WCI受旁瓣影响大无法使用，选取半径内WCI波束阵列图)、以换能器为原点的扇面图 (Depth-Acrosstrack空间，D-A空间)及投影坐标下多波束WCI(Depth-Projection空间，D-P空间)(图6)；
(3)噪声初步压制：去除船舶、换能器自身在水体浅表层数十次采样序列中产生的噪声，并以中央波束为分界线，将两侧呈镜面对称WC数据对称作差以初步压制WCI中噪声(图7)；
(4)多波束WCI滤波：对滤除表层噪声并对称作差法初步去噪的WCI采用图像滤波方法，平滑增强气体羽状流目标，完成数据前处理工作(图8)。
(5)获取二值化分割后的水体点云数据：基于大津算法(OSTU)逐一对全部Ping多波束WCI进行二值化处理，并设置参数w供用户设置，以改善二值化结果(图9)。
(6)多波束WCI第三方干扰去除：用于多波束WCI有第三方声呐干扰时，先选取NPing不包含羽状流的二值化WCI，采用开操作设置像素阈值P滤除小连通分量的方式提取第三方声呐干扰，并对N Ping 提取结果求并生成模板(图10)。逐一滤除全部Ping二值化WCI中与模板相交部分，以获取初步分离气体羽状流后的水体点云；若无则多波束WCI二值化处理后即可获取初步分离气体羽状流的水体点云(图 11)。
(7)Kmeans首次尝试去除上层噪点：将每Ping初步分离气体羽状流的水体点云在以换能器为原点的相对坐标空间叠加(D-A空间，图12)，首次尝试以点云水深值为特征值，通过Kmeans将点云分为2类，若类间聚类中心比值≥≥r去除上层分布噪点，≤r则不去除(图13)；
(8)LOF迭代去除离散分布噪点：利用D-A空间剩余噪点为离散分布且单位区域噪点数量远少于羽状流特点，采用从点云中随机选取一定数量点的方法使噪点在抽稀后表现为离群点以满足LOF方法要求 (图15)，再通过LOF方法迭代去除离散分布噪点(图16)，LOF迭代去噪具体流程如图14。。
(9)Kmeans再次尝试去除上层噪点：再次尝试以点云水深值为特征值，通过Kmeans将点云分为2 类，若类间聚类中心比值≥≥r去除上层分布噪点，≤r则不去除；
(10)绘图：利用提取结果和解析的数据，可实现气体羽状流点云图、多波束各类WCI、地形图的成图。
本发明采用如图28所示的处理方法：
(1)解析原始数据
多波束测深系统保存的各类型数据包均由导航、时间、姿态、水深、水体、POS等信息构成，并以二进制存储方式对系统输入数据进行存储，需将原始数据解析转化为可用格式数据。本方法以挪威Kongsberg 公司EM302型深水多波束系统为例，进行了具体实施。该系统可形成288个波束，系统发射波束最大扇面开角180°，形成波束宽度1°×1°。EM型多波束系统的数据存储格式包括ALL格式和WCD格式，其中 ALL文件包含测深、定位、日期、时间、姿态传感器、罗经等数据包，WCD可只用于存储水体相关信息，也可包含ALL文件中信息。首先，需根据文件格式，编写程序，解析提取原始数据信息，保存为可用的格式。
(2)水深及水体采样点计算
根据原始数据中水深地形及POS数据可完成水深数据投影坐标计算，公式如下：
Figure RE-GDA0002867399030000111
[THETA，RHO]＝cart2pol(AcrosstrackY，AlongtrackX)
Figure RE-GDA0002867399030000112
S_East＝S_Eastsonar+P_prcEast
S_North＝S_Northsonar+P_prcNorth
式中：P_prc为每Ping数据获取时投影坐标(P_prcEast、P_prcNorth)、航向(P_prcHead)、子午线收敛角(P_prcCon)；posdown、posup为与每Ping水深数据记录时间最邻近的两组POS点；AcrosstrackY， AlongtrackX为以换能器为原点采样点横向、纵向相对坐标，m；S_Eastsonar，S_Northsonar为水深采样点到换能器距离，m。编程计算各地形采样点投影位置即可获取水深地形图(图5)。
结合导航数据包、水柱数据包内存储数据，可计算水体后向散射强度采样点以换能器为原点的相对坐标，公式如下：
PBS_SR＝(PBS_idx+PBS_sRSNum)×c/(2·f)
PBS_H＝-PBS_SR×cos(PBS_BPA)
PBS_X＝-PBS_SR×sin(PBS_BPA)
获取水体后向散射强度采样点相对坐标后，结合每Ping数据换能器投影坐标，可计算获取水体后向散射强度采样点投影坐标，公式如下：
Figure RE-GDA0002867399030000121
PBS_East＝P_prcEast+PBS_X×cos(theta)
PBS_North＝P_prcNorth+PBS_X×sin(theta)
式中：PBS_SR为水体采样点到换能器距离，m；BS_idx为每柱接收波束回波采样点编号；PBS_sRSNum 为每柱接收波束起始采样点编号；c为声速，m；f为采频率；BS_BPA为波束发射角；PBS_X，PBS_H为水体后向散射强度采样点沿船横向位置及水深，m；sonarHeadingOffset为换能器安装角度，°；theta为波束方位角，°；PBS_East、PBS_North为水体采样点投影坐标，m。编程计算各水体采样点相对坐标和投影坐标后，即可获取测区完整点云数据、各类多波束WCI(图6)。
(3)噪声初步压制
针对测区内噪声以中央波束为轴对称分布特征，对多波束WC数据单Ping内呈镜面对称位置水体后向反射强度采样点进行Pearson相关性分析(图4)，并求和取均值可见单Ping内镜面对称位置水体后向散射强度采样点呈中等程度正相关。以中央波束为分界线，将两侧呈镜面对称WC数据对称作差初步压制 WCI中噪声(图7)。
(4)WCI滤波
通过图像滤波方法对多波束WCI数据进行预处理，进一步降低噪声影响，并对WCI内气体羽状流部分目标增强(图8，示例为巴特沃斯滤波结果)。
(5)多波束WCI二值化
以中央波束为轴将每Ping WCI图像分成两侧，采用OSTU算法初步确定两侧二值化阈值，并设置系数w作为此阈值的系数，以改善二值化结果，公式：
Figure RE-GDA0002867399030000122
式中：mG为WCI平均灰度；概率P1(k)为灰度阈值为k时被分到小于k一类的概率；mk为灰度级k的累加均值；δB 2为类间方差。从式中不难看出，类间方差越大，则可分性越好。因此，将每Ping单侧多波束WCI中包含的灰度级依次代入上式，当δB 2取最大时的k*即为此算法下计算获取阈值，再乘以预先设置的系数w，即可依此灰度值完成图像二值化。通过上述方法分别计算每Ping WCI两侧阈值完成二值化，并按需求设置合理阈值P滤除像素数量小的连通分区后拼接，完成多波束WCI二值化(图9)。
(6)多波束WCI第三方干扰去除
若多波束WCI有第三方声呐干扰，则选取N Ping不包含羽状流的二值化WCI图像执行开操作，即对图像腐蚀，然后膨胀的处理方式。其中，图像腐蚀、膨胀公式分别为：
Figure RE-GDA0002867399030000131
Figure RE-GDA0002867399030000132
式中：Z为某Ping WCI；A为Z中包含的连通分量；E为腐蚀后Z中包含的连通分量；SE为结构体；开操作处理后WCI中指示第三方声呐干扰的连通分区像素数量较多，设置合理阈值P滤除像素数量小的连通分区。然后对N Ping处理结果求并完成第三方声呐干扰模板提取(图10)，逐一滤除全部Ping二值化WCI中与模板相交部分，以获取初步分离气体羽状流后的水体点云，若无第三方声呐干扰，则通过多波束WCI二值化结果即可获取初步分离气体羽状流的水体点云(图11)。
(7)数据挖掘去噪
WCI中噪点数据主要分为上下两部分：上层噪点数据可能来源海洋风浪、浮游动植物、悬浮物等产生的海洋环境噪声，该类噪声在不同型号多波束、不同水域采集的WCI中均有存在，且大区域分布，可采用经典无监督聚类方法K-means以点云水深值作为特征值，将点云数据分成2类。上下层数据往往分割明显，但仍可能与噪点接近，故设置最上层簇中心与底层簇中心比值大于r为限，去除聚类后水体最上层由噪声干扰产生的数据点(图13)；下层噪点主要是在OSTU二值化分割开操作后没有完全消除的旁瓣、环境噪声、第三方声呐等产生的噪点，往往为离散分布噪点，且单位区域数量远小于气体羽状流部分。从图 12点云数据叠加图可见，残余噪点区域内点间距离与气体羽状流区域内点间距离相近，无法满足LOF算法异常数据点为离群点要求，直接使用LOF算法无法滤除点云中残余噪点。但利用残余噪点为离散分布噪点且单位区域数量远小于气体羽状流数据点特征，从点云中随机选取一定数量点的方法使噪点在抽稀后的 Depth-Acrosstrack空间中表现为离群点，然后采用LOF迭代去除噪点(图15)，LOF迭代去噪具体方法如下(图14)：
①点云数据输入，将OSTU二值化(或OSTU二值化并去除第三方声呐干扰)后获取的初步分离羽状流点云点云投影到Depth-Acrosstrack空间，构成初始输入数据。
②设置参数，设置参数邻域k、LOF阈值thresh、迭代次数iter、随机抽取样本点数量num。
③随机抽取样本点集合O，从点云数据中随机抽取num个样本点，形成样本点集合O＝{o1，o2， o3，…on}，其中oi是集合O中第i个样本点。
④计算第k距离，计算集合O每个样本点oi的第k距离distk(oi)(k-distance)，即距oi第k远(不包括 oi)的样本点的距离。
⑤计算k邻域内样本点真实距离，计算集合O每个样本点oi的k邻域内(Nk(oi))点oi′到点oi的距离 dist(oi，oi′)，即点oi第k距离以内的所有点，包括第k距离在内的点到oi的欧式距离。
⑥计算可达距离，计算集合O每个样本点oi的k邻域内(Nk(oi))点oi′到点oi的第k可达距离，其定义如下：
reachdistk(oi←oi′)＝max{distk(oi′)，dist(oi，oi′)}；
即点oi′到点oi的第k可达距离，可以看出其值至少是oi′的第k距离，或者为oi、oi′间的真实距离。
⑦计算局部可达密度，计算集合O中每个样本点oi的局部可达密度(Locareachability density，lrd)，其定义如下：
Figure RE-GDA0002867399030000141
即点oi的第k邻域内点oi′到oi平均可达距离的倒数。其意在量化oi与其k邻域内点之间密度差异，若lrdk越高则越可能属于同一簇点，越低则越可能为离群点。
⑧计算LOF，根据(7)中计算所得lrd，计算集合O中每个样本点oi的局部离群因子(LOF)，其定义如下：
Figure RE-GDA0002867399030000142
表示点oi的邻域点Nk(oi)的局部可达密度与点oi的局部可达密度之比的平均数，若LOF值越接近于1，两者的密度越接近；若LOF值小于1，则oi密度高于其邻域点；若LOF值大于1，则值越大说明点oi密度小于邻域点，越可能是离群点，即更可能为噪点。
⑨滤除LOF异常高点，根据预先设置好的LOF阈值thresh，滤除大于thresh的点云数据，若迭代次数小于iter则转至第(3)步，否则顺序执行。
⑩有禁忌搜索随机抽取计算
Figure RE-GDA0002867399030000143
滤除空间内未去除噪点，由于上述LOF迭代去噪过程中采用随机抽取点云数据的方式将噪点转化为离群点，伴随着迭代过程的进行，噪点被选取得可能性越来越小，导致迭代去噪处理过程中有可能有部分噪点没有被选取未被去出。因此，在首次LOF迭代去噪完成后，采用有禁忌搜索随机选取集合O对全部点云重新完整迭代计算LOF滤除首次迭代后依旧未完全去除的噪点，完善提取结果后结束。通过Kmeans与LOF两种数据挖掘方法相结合，即可去除点与数据中绝大部分噪点 (本发明示例Kmeans满足2类聚类中心比值大于r要求，若无法满足用户可自主判断去除上层噪点)，完成气体羽状流提取(图16)。
(8)绘图，利用提取结果和解析的数据，可实现气体羽状流点云图、多波束各类WCI、水深地形图的成图。
证明部分(具体实施例/实验/仿真/能够证明本发明创造性的正面实验数据等)
为证明本方法的可行性，本文选取另一地区另一组示例作为验证例加以证明，原始水柱数据概况见图 17，测区水深地形见图18。对该组数据进行按图28流程所示提取WCI中气体羽状流。从图19对每Ping 水柱数据对称位置声散射强度相关性进行分析，可见在没有第三方声呐干扰的情况下，相关系数提高且成正相关；由于该组数据无第三方声呐干扰的影响，无需提取模板去除该类噪声。在WCI预处理完成之后 (图20-21，验证例滤波截至频率为100)，可直接设置系数w完成OSTU二值化(图22)，获取初步分离气体羽状流的点云数据(图23)；在LOF迭代去噪前后尝试Kmeans去除上层噪点，若类间聚类中心比值≥r去除上层分布噪点(图24)，≤r不去除；将每Ping二值化分离结果叠加到D-A空间内(25)，LOF 迭代去除剩余噪点(图26)，即可完成羽状流提取(图27)，从最终提取结果可见方法的可行性。
应当注意，本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的普通技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现，例如在诸如磁盘、CD或DVD-ROM的载体介质、诸如只读存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、诸如逻辑芯片、晶体管等的半导体、或者诸如现场可编程门阵列、可编程逻辑设备等的可编程硬件设备的硬件电路实现，也可以用由各种类型的处理器执行的软件实现，也可以由上述硬件电路和软件的结合例如固件来实现。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。

Claims (10)

1.一种多波束水柱数据中海底气体羽状流处理方法，其特征在于，所述多波束水柱数据中海底气体羽状流处理方法包括：解析提取多波束原始数据信息，根据数据记录信息完成水体、水深采样点归位，计算回波强度值真实的地理位置，以获取测区完整水体点云信息、水深地形图、各类多波束WCI；通过对称作差、图像滤波实现多波束水柱影像提取气体羽状流前处理工作；若WCI无第三方声呐干扰，图像二值化后获取初步分离气体羽状流的水体点云数据；若WCI有第三方声呐干扰，则先选取N Ping不包含羽状流的二值化WCI，采用开操作设置像素阈值滤除小连通分量的方式将第三方声呐干扰提取出来，并对N Ping提取结果求并生成模板，滤除每Ping二值化WCI中与模板相交部分，以获取初步分离气体羽状流后的水体点云数据；引入Kmeans、局部异常因子方法进一步迭代去噪提取气体羽状流；绘制各类多波束WCI、气体羽状流点云图、多波束水深地形图。
2.如权利要求1所述的多波束水柱数据中海底气体羽状流处理方法，其特征在于，所述多波束水柱数据中海底气体羽状流处理方法还包括：
第一步，根据数据存储方式解析提取原始数据，转换为可用格式数据；
第二步，转换后的数据显示多波束记录了水深、底质和水体回波强度、采样频率、声速及波束入射角等信息，利用提取到的信息，对每一Ping数据进行采样点归位计算，得到每一水体、水深采样点实际位置信息，进而获取测区水深地形图、完整水体点云信息、多波束WCI波束阵列图、扇面图、航向剖面图；扇面图包括投影和先对坐标；
第三步，船舶、换能器自身产生的噪声仅在水体浅表层数十次采样序列中存在，直接将浅表层水体一定数量等时采样序列去除，并以中央波束为分界线，将两侧呈镜面对称WC数据对称作差初步压制WCI中噪声；
第四步，对滤除表层噪声并对称作差法初步去噪的WCI采用图像滤波方法，平滑增强气体羽状流目标，完成数据前处理；
第五步，基于大津算法OSTU逐一对全部Ping多波束WCI进行二值化处理，并设置系数w供用户设置；
第六步，若多波束WCI有第三方声呐干扰，则先选取N Ping不包含羽状流的二值化WCI，采用开操作设置像素阈值P滤除小连通分量的方式提取第三方声呐干扰，并对N Ping提取结果求并生成模板，滤除全部二值化WCI每Ping中与模板相交部分，以获取初步分离气体羽状流后的水体点云；若无则第五步处理结果即为初步分离气体羽状流的水体点云数据；
第七步，将每Ping初步分离气体羽状流的水体点云在以换能器为原点的相对坐标空间叠加，首次尝试以点云水深值为特征值，通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤r则不去除；
第八步，采用LOF方法迭代去除剩余区域噪点；
第九步，再次尝试以点云水深值为特征值，通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤r则不去除，完成羽状流提取；
第十步，利用提取结果和解析的数据，实现气体羽状流点云图、多波束各类WCI、地形图的成图。
3.如权利要求2所述的多波束声呐水柱数据中海底气体羽状流处理方法，其特征在于，所述第二步的水深及水体采样点计算包括：
根据原始数据中水深地形及POS数据完成水深数据投影坐标计算，公式：
Figure RE-FDA0002867399020000021
[THETA，RHO]＝cart2pol(AcrosstrackY，AlongtrackX)
Figure RE-FDA0002867399020000022
Figure RE-FDA0002867399020000023
Figure RE-FDA0002867399020000024
式中：P_prc为每Ping数据获取时投影坐标(P_prcEast、P_prcNorth)、航向(P_prcHead)、子午线收敛角(P_prcCon)；posdown、posup为与每Ping水深数据记录时间最邻近的两组POS点；AcrosstrackY，AlongtrackX为以换能器为原点采样点横向、纵向相对坐标，m；S_Eastsonar，S_Northsonar为水深采样点到换能器距离，m。编程计算各地形采样点投影位置即可获取测区水深地形图；
结合导航数据包、水柱数据包内存储数据，计算水体后向散射强度采样点以换能器为原点的相对坐标，公式：
PBS_SR＝(PBS_idx+PBS_sRSNum)×c/(2·f)
PBS_H＝-PBS_SR×cos(PBS_BPA)
PBS_X＝-PBS_SR×sin(PBS_BPA)；
获取水体后向散射强度采样点相对坐标后，结合每Ping数据换能器投影坐标，计算获取水体后向散射强度采样点投影坐标，公式：
Figure RE-FDA0002867399020000031
Figure RE-FDA0002867399020000032
Figure RE-FDA0002867399020000033
式中：PBS_SR为水体采样点到换能器距离，m；BS_idx为每柱接收波束回波采样点编号；PBS_sRSNum为每柱接收波束起始采样点编号；c为声速，m；f为采频率；BS_BPA为波束发射角；PBS_X，PBS_H为水体后向散射强度采样点沿船横向位置及水深，m；sonarHeadingOffset为换能器安装角度，°；theta为波束方位角，°；PBS_East、PBS_North为水体采样点投影坐标，m。编程计算各水体采样点相对坐标和投影坐标后，获取测区完整水体点云信息、多波束WCI波束阵列图、扇面图、航向剖面图。
4.如权利要求2所述的多波束水柱数据中海底气体羽状流处理方法，其特征在于，所述第三步噪声初步压制包括：
针对测区内噪声以中央波束为轴对称分布特征，对多波束WC数据单Ping内呈镜面对称位置采样点水体后向反射强度进行Pearson相关性分析，并求和取均值可见单Ping内镜面对称位置水体后向散射强度采样点呈中等程度正相关，以中央波束为分界线，将两侧呈镜面对称WC数据对称作差初步压制WCI中噪声；
5.如权利要求2所述的多波束声呐水柱数据中海底气体羽状流处理方法，其特征在于，所述第四步WCI滤波包括：通过图像滤波方法对多波束WCI数据进行滤波，平滑增强气体羽状流目标，完成数据前处理。
6.如权利要求2所述的多波束水柱数据中海底气体羽状流处理方法，其特征在于，所述第五步多波束WCI二值化包括：
以中央波束为轴将每Ping WCI图像分成两侧，采用OSTU算法初步确定两侧二值化阈值，并设置系数w作为此阈值的系数，以改善二值化结果，公式：
Figure RE-FDA0002867399020000041
式中：mG为WCI平均灰度；概率P1(k)为灰度阈值为k时被分到小于k一类的概率；mk为灰度级k的累加均值；δB 2为类间方差，将每Ping单侧多波束WCI中包含的灰度级依次代入上式，当δB 2取最大时的k*即为此算法下计算获取的阈值，再乘以预先设置的w，即可得到单侧WCI二值化阈值。分别计算每Ping WCI两侧阈值完成二值化后拼接，即可完成多波束WCI二值化。
7.如权利要求2所述的多波束水柱数据中海底气体羽状流处理方法，其特征在于，所述第六步多波束WCI第三方干扰去除包括：
若多波束WCI有第三方声呐干扰，则选取N Ping不包含羽状流的二值化WCI图像执行开操作，即对图像腐蚀，然后膨胀的处理方式，图像腐蚀、膨胀公式分别为：
Figure RE-FDA0002867399020000042
Figure RE-FDA0002867399020000043
式中：Z为某Ping WCI；A为Z中包含的连通分量；E为腐蚀后Z中包含的连通分量；SE为结构体；开操作处理后WCI中指示第三方声呐干扰的连通分区像素数量较多，设置合理阈值P滤除像素数量小的连通分区。然后对N Ping处理结果求并完成第三方声呐干扰模板提取，滤除全部二值化WCI中与模板相交部分，以获取初步分离气体羽状流后的水体点云，若无第三方声呐干扰，则通过多波束WCI二值化结果即可获取初步分离气体羽状流的水体点云；
所述第七步和第九步Kmeans去噪包括：
在LOF迭代前后分别尝试以点云水深值为特征值，采用无监督聚类方法K-means将点云分成2类，设置最上层聚类中心与下层簇中心比值大于r为限。当比值大于r时，去除聚类后水体最上层由噪声干扰产生的数据点；
8.如权利要求2所述的多波束水柱数据中海底气体羽状流处理方法，其特征在于，所述第八步LOF迭代去噪包括：
将每Ping初步分离气体羽状流的水体点云在以换能器为原点的扇形空间坐标D-A空间下叠加，剩余下层噪点是在获取初步分离气体羽状流的水体点云时没有完全消除的旁瓣、环境噪声、第三方声呐等类型噪点，往往为离散分布噪点。利用下层噪点离散且叠加在D-A空间后，单位区域数量远小于气体羽状流数据点特征，通过从点云中随机选取一定数量点的方法使噪点在抽稀后的D-A空间中表现为离群点，再采用LOF迭代方法去除该类噪点。所述LOF迭代去噪具体方法如下：
①点云数据输入，将每Ping OSTU二值化或OSTU二值化并去除第三方声呐干扰后获取的初步分离羽状流点云到D-A空间，构成初始输入数据；
②设置参数，设置参数邻域k、LOF阈值thresh、迭代次数iter、随机抽取样本点数量num；
③随机抽取样本点集合O，从点云中随机抽取num个样本点，形成样本点集合O＝{o1，o2，o3，…on}，其中oi是集合O中第i个样本点；
④计算第k距离，计算集合O每个样本点oi的第k距离distk(oi)，即距oi第k远的样本点的距离；
⑤计算k邻域内样本点真实距离，计算集合O每个样本点oi的k邻域内(Nk(oi))点oi′到点oi的距离dist(oi，oi′)，即点oi第k距离以内的所有点，包括第k距离在内的点到oi的欧式距离；
⑥计算可达距离，计算集合O每个样本点oi的k邻域内(Nk(oi))点oi′到点oi的第k可达距离，其定义如下：
reachdistk(oi←oi′)＝max{distk(oi′)，dist(oi，oi′)}；
即点oi′到点oi的第k可达距离，看出其值至少是oi′的第k距离，或者为oi、oi′间的真实距离；
⑦计算局部可达密度，计算集合O中每个样本点oi的局部可达密度lrd，其定义如下：
Figure RE-FDA0002867399020000061
即点oi的第k邻域内点oi′到oi平均可达距离的倒数，其意在量化oi与其k邻域内点之间密度差异，若lrdk越高则越可能属于同一簇点，越低则越可能为离群点；
⑧计算LOF，根据计算所得lrd，计算集合O中每个样本点oi的局部离群因子LOF，其定义如下：
Figure RE-FDA0002867399020000062
表示点oi的邻域点Nk(oi)的局部可达密度与点oi的局部可达密度之比的平均数，若LOF值越接近于1，两者的密度越接近；若LOF值小于1，则oi密度高于其邻域点；若LOF值大于1，则值越大说明点oi密度小于邻域点，越可能是离群点，即更可能为噪点；
⑨滤除LOF异常高点，根据预先设置好的LOF阈值thresh，滤除大于thresh的点云数据，若迭代次数小于iter则转至③，否则顺序执行；
⑩有禁忌搜索随机抽取计算LOFk(oi)滤除空间内未去除噪点，在首次LOF迭代去噪完成后，采用有禁忌搜索随机选取集合O对全部点云重新完整迭代计算LOF滤除首次迭代后依旧未完全去除的噪点，完善提取结果。
9.一种实施权利要求1～8任意一项所述多波束水柱数据中海底气体羽状流处理方法的多波束水柱数据中海底气体羽状流处理系统，其特征在于，所述多波束水柱数据中海底气体羽状流处理系统包括：
数据解析模块，用于不同公司不同型号不同存储模式多波束原始数据的解析提取，将其转换为可用格式数据；
水体及水深采样点归位模块，用于根据原始数据提取到的多波束水深、底质、水体回波强度、采样频率、声速及波束入射角等信息，对每一Ping数据进行采样点归位计算，得到每一水体、水深采样点实际位置信息，进而获取测区水深地形图、水体点云信息、多波束WCI波束阵列图、扇面图(投影和相对坐标)、航向剖面图；
噪声初步压制模块，用于去除船舶、换能器自身在水体浅表层数十次采样序列中产生的噪声；以中央波束为分界线，将两侧WC数据对称作差以初步压制WCI中噪声；
多波束WCI滤波模块，用于通过图像滤波方法对多波束WCI数据进行滤波，平滑增强气体羽状流目标，完成数据前处理；
多波束WCI二值化模块，用于基于大津算法OSTU逐一对全部Ping多波束WCI进行二值化处理，并设置系数w供用户设置，以改善二值化结果；
多波束WCI第三方干扰去除模块：用于多波束WCI有第三方声呐干扰时，先选取N Ping不包含羽状流的二值化WCI，采用开操作，并设置像素阈值P滤除小连通分量的方式提取第三方声呐干扰，并对N Ping提取结果求并生成模板。滤除全部二值化WCI每Ping中与模板相交部分，以获取初步分离气体羽状流后的水体点云数据；若无则多波束WCI二值化处理后即可获取初步分离气体羽状流的水体点云；
数据挖掘去噪模块：用于去除初步分离气体羽状流水体点云中的剩余噪点，LOF迭代去噪前，以点云水深值为特征值，尝试通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤r不去除。然后，采用LOF方法迭代去除剩余离散分布噪点，再次尝试通过Kmeans将点云分为2类，若类间聚类中心比值≥r去除上层分布噪点，≤r不去除，完成气体羽状流提取；
绘图模块，用于利用提取结果和解析的数据，实现气体羽状流点云图、多波束各类WCI、地形图的成图。
10.一种用于海底输气管道的渗漏提取方法，其特征在于，所述用于海底输气管道的渗漏检测方法包括以下步骤：解析提取多波束原始数据信息，根据数据记录信息完成水体、水深采样点归位，计算回波强度值真实的地理位置，获取测区完整水体点云信息、水深地形图、各类多波束WCI；通过对称作差、图像滤波实现多波束水柱影像提取油气渗漏前处理工作；若WCI无第三方声呐干扰，图像二值化后即可获取初步分离油气渗漏的水体点云；若WCI有第三方声呐干扰，则先选取N Ping不包含油气渗漏的二值化WCI，采用开操作设置像素阈值滤除小连通分量的方式将第三方声呐干扰提取出来，并对N Ping提取结果求并生成模板，逐一滤除全部Ping二值化WCI中与模板相交部分，以获取初步分离油气渗漏后的水体点云；引入Kmeans、局部异常因子方法进一步迭代去噪提取油气渗漏点；绘制各类多波束WCI、油气渗漏点云图、多波束水深地形图。
CN202011278968.2A 2020-11-16 2020-11-16 多波束水柱数据中海底气体羽状流处理方法、系统及应用 Active CN112529841B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011278968.2A CN112529841B (zh) 2020-11-16 2020-11-16 多波束水柱数据中海底气体羽状流处理方法、系统及应用

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011278968.2A CN112529841B (zh) 2020-11-16 2020-11-16 多波束水柱数据中海底气体羽状流处理方法、系统及应用

Publications (2)

Publication Number Publication Date
CN112529841A true CN112529841A (zh) 2021-03-19
CN112529841B CN112529841B (zh) 2023-01-31

Family

ID=74980901

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011278968.2A Active CN112529841B (zh) 2020-11-16 2020-11-16 多波束水柱数据中海底气体羽状流处理方法、系统及应用

Country Status (1)

Country Link
CN (1) CN112529841B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113175974A (zh) * 2021-03-22 2021-07-27 中交华南勘察测绘科技有限公司 智能水深测量的成图方法、装置及系统
CN115880189A (zh) * 2023-02-22 2023-03-31 山东科技大学 一种海底地形多波束点云滤波方法
CN117538881A (zh) * 2024-01-10 2024-02-09 海底鹰深海科技股份有限公司 一种声呐水体成像波束形成方法、系统、设备及介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106097315A (zh) * 2016-06-03 2016-11-09 河海大学常州校区 一种基于声呐图像的水下构筑物裂缝提取方法
CN109035224A (zh) * 2018-07-11 2018-12-18 哈尔滨工程大学 一种基于多波束点云的海底管道检测与三维重建方法
US20190114747A1 (en) * 2016-04-07 2019-04-18 Carmel Haifa University Economic Corporation Ltd. Image dehazing and restoration
CN110675410A (zh) * 2019-09-25 2020-01-10 江苏海洋大学 基于选择性搜索算法的侧扫声呐沉船目标非监督探测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20190114747A1 (en) * 2016-04-07 2019-04-18 Carmel Haifa University Economic Corporation Ltd. Image dehazing and restoration
CN106097315A (zh) * 2016-06-03 2016-11-09 河海大学常州校区 一种基于声呐图像的水下构筑物裂缝提取方法
CN109035224A (zh) * 2018-07-11 2018-12-18 哈尔滨工程大学 一种基于多波束点云的海底管道检测与三维重建方法
CN110675410A (zh) * 2019-09-25 2020-01-10 江苏海洋大学 基于选择性搜索算法的侧扫声呐沉船目标非监督探测方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
CHAO XU ET AL。: "Optical Flow-Based Detection of Gas Leaks from Pipelines Using Multibeam Water Column Images", 《REMOTE SENSING》 *
GIACOMO MONTEREALE-GAVAZZI ET AL.: "Seafloor change detection using multibeam echosounder backscatter:case study on the Belgian part of the North Sea", 《MAR GEOPHYS RES》 *
JIANHU ZHAO ET AL。: "Automatic Detection and Segmentation on Gas Plumes from Multibeam Water Column Images", 《REMOTE SENSING》 *
蒲定 等: "水中气体目标的多波束声呐成像与检测算法", 《应用科技》 *
龙睿捷 等: "多波束水柱影像中三维羽流数据提取方法研究", 《江西科学》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113175974A (zh) * 2021-03-22 2021-07-27 中交华南勘察测绘科技有限公司 智能水深测量的成图方法、装置及系统
CN115880189A (zh) * 2023-02-22 2023-03-31 山东科技大学 一种海底地形多波束点云滤波方法
CN115880189B (zh) * 2023-02-22 2023-05-30 山东科技大学 一种海底地形多波束点云滤波方法
CN117538881A (zh) * 2024-01-10 2024-02-09 海底鹰深海科技股份有限公司 一种声呐水体成像波束形成方法、系统、设备及介质
CN117538881B (zh) * 2024-01-10 2024-05-07 海底鹰深海科技股份有限公司 一种声呐水体成像波束形成方法、系统、设备及介质

Also Published As

Publication number Publication date
CN112529841B (zh) 2023-01-31

Similar Documents

Publication Publication Date Title
CN112529841B (zh) 多波束水柱数据中海底气体羽状流处理方法、系统及应用
Calder et al. Automatic processing of high‐rate, high‐density multibeam echosounder data
Luo et al. Sediment classification of small-size seabed acoustic images using convolutional neural networks
Singh et al. Microbathymetric mapping from underwater vehicles in the deep ocean
Fox et al. Detection of changes in ridge‐crest morphology using repeated multibeam sonar surveys
CN105787886A (zh) 一种基于多波束图像声纳的实时图像处理方法
Zhang et al. Submarine pipeline tracking technology based on AUVs with forward looking sonar
Villar et al. Evaluation of an efficient approach for target tracking from acoustic imagery for the perception system of an autonomous underwater vehicle
CN112539886B (zh) 一种海底气体羽状流提取方法及应用
Mitchell Processing and analysis of Simrad multibeam sonar data
Vàsquez Tuning the CARIS implementation of CUBE for Patagonian Waters
CN115422981A (zh) 面向单频机载激光测深数据的水陆分类方法、系统及应用
Li et al. Automatic detection of pipelines from sub-bottom profiler sonar images
Villar et al. A framework for acoustic segmentation using order statistic-constant false alarm rate in two dimensions from sidescan sonar data
Foresti et al. A voting-based approach for fast object recognition in underwater acoustic images
CN108460773B (zh) 一种基于偏移场水平集的声纳图像分割方法
CN113703050B (zh) 一种深海地震垂直缆二次定位方法
CN115690035A (zh) 基于声纳图形的水下结构缺陷检测方法、装置、设备和介质
Veloso-Alarcón et al. Quantitatively monitoring bubble-flow at a seep site offshore Oregon: Field trials and methodological advances for parallel optical and hydroacoustical measurements
Wang et al. Seafloor classification based on deep-sea multibeam data—Application to the southwest Indian Ridge at 50.47° E
Zhang et al. Reconstruction of large complex sand-wave bathymetry with adaptive partitioning combining satellite imagery and sparse multi-beam data
Rizaev et al. The effect of sea state on ship wake detectability in simulated SAR imagery
Li et al. Surface extraction and segmentation from 3-D underwater sub-bottom point clouds using enhancement filtering and global energy optimization
Casalbore Multibeam Echosounder
Lingsch et al. Acoustic imagery using a multibeam bathymetric system

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant