CN109101996A - 一种多类型探测传感器综合观测的海底热液探测方法 - Google Patents
一种多类型探测传感器综合观测的海底热液探测方法 Download PDFInfo
- Publication number
- CN109101996A CN109101996A CN201810735835.XA CN201810735835A CN109101996A CN 109101996 A CN109101996 A CN 109101996A CN 201810735835 A CN201810735835 A CN 201810735835A CN 109101996 A CN109101996 A CN 109101996A
- Authority
- CN
- China
- Prior art keywords
- data
- value
- time
- detection
- target
- 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
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/25—Fusion techniques
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
- G01D21/00—Measuring or testing not otherwise provided for
- G01D21/02—Measuring two or more variables by means not covered by a single other subclass
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Computer Hardware Design (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Testing Or Calibration Of Command Recording Devices (AREA)
Abstract
本发明公开了一种多类型探测传感器综合观测的海底热液探测方法。专人负责某一项探测参数的后处理工作，数据融合性较差。本发明步骤：定位数据准备；匹配浊度传感器探测资料的位置信息；浊度传感器探测资料的处理；匹配化学传感器探测资料的位置信息；化学传感器探测资料的处理；匹配甲烷传感器探测资料的位置信息；甲烷传感器探测资料的处理；成图。本发明可对多种热液调查过程中采集的多类型外业数据进行融合并快速处理，提高了数据处理的效率，操作方便，费用相对较低。
Description
技术领域
本发明属于海底探测技术领域，具体涉及一种以海底热液流体为主要目标的探测采集方法及后处理技术，可同时对温度、浊度、氧化还原电位、甲烷等异常信息进行观测和判断。
背景技术
海底热液活动喷出的热液流体是现代海底热液探测的核心要素，也是大洋中脊向海洋输送物质和能量的关键介质。热液流体，从海底热液喷口喷出并与周围海水相遇后迅速混合，改变了周围海水的温度、盐度、浊度等属性值。海底热液成因的多金属硫化物矿床作为一种海底资源，目前已成为各国科学家研究的热点。
当前对海底热液活动研究的主要方法是采用探测海水温度、浊度、氧化还原电位、硫化氢、pH、甲烷、铁锰金属元素异常等手段，来实现对热液流体的定位和分布研究，从而锁定热液喷口。近年来我国在热液调查上的投入力度越来越大，发现了一系列的海底热液区。但在现场海洋调查中，涉及海洋学科领域多，海洋数据获取手段迥异，直接观测的海洋技术中与热液调查密切相关的有：温盐深仪CTD、浊度传感器、化学传感器、甲烷传感器；这些传感器有各自不同的资料采集系统，记录方法也不同，数据的格式也不统一，在现场调查中需要专人单独负责某一项探测参数的后处理工作，由此探测参数的处理方法也会因人而异。因此在综合判断热液流体的定位与分布过程中，数据融合性较差，往往会导致现场工作效率低下，因此亟需提出一种多类型探测传感器综合观测的海底热液探测方法。
发明内容
本发明的目的在于针对现有多参数探测数据后处理方法依赖多人负责、方法不统一、数据融合性较差、工作效率低下的问题，提出一种温度、浊度、氧化还原电位、硫化氢、pH、甲烷等异常信息综合处理的数据处理方法，以实现现场探测工作的高效性、准确性。
本发明包括下列步骤：
步骤1：定位数据准备工作
(1.1)整理GPS定位数据
GPS定位数据，采样周期为1秒，有时间、经度、纬度三列数据。其中原始的GPS定位数据中，时间序列是字符串格式，表现形式为“20170930140500”，即“YYYYMMDDHHMMSS”，分别表示“年份”、“月份”、“日期”、“小时”、“分钟”和“秒”。时间序列是多参数探测数据进行融合的唯一标准，因此在数据准备阶段，通过MATLAB软件的datenum函数将字符串的时间序列转换成具体的日期值。根据实际测线的探测时间，选取该时间范围内的GPS定位数据存入矩阵ship(i1，j1)中，i1＝1，...，nxl，j1＝1，...，nyl，其中nxl为矩阵总行数，nyl为矩阵总列数，i1、j1、nx1、ny1均为自然数。在海上作业现场，往往只需要对刚结束调查的单独一条测线数据进行处理和分析，因此按测线逐条处理即可；而在实验室后期处理中，可将调查航次中其中一个航段的GPS数据进行汇总，统一处理，方便探测数据的位置匹配工作。
(1.2)整理超短基线定位系统(USBL)数据
在海底热液流体探测工作(包括综合摄像拖体下放和回收的过程，整个过程中，综合探测传感器持续探测数据)中，以温度、浊度、氧化还原电位、甲烷、硫化氢、pH参数为主要指标，分别在距离综合摄像拖体0m(即加载在综合摄像拖体本体上)、5m、30m、80m、L1、L2的高度加挂综合探测传感器(包括浊度传感器、化学传感器和甲烷传感器)，L1在200-300m中取值，L2在300-500m中取值。当综合摄像拖体近底作业时，超短基线定位系统开始记录综合摄像拖体的位置信息。超短基线定位系统的采样周期不固定，通常在9-11秒之间浮动，有时间、经度、纬度三列数据，与GPS数据的处理过程一样，将超短基线定位系统的时间序列通过MATLAB软件的datenum函数转换成具体的日期值，并将超短基线定位系统数据存入矩阵target(i2，j2)中，其中i2＝1，...，nx2，j2取值为1、2或3分别对应时间、经度、纬度，target(i2，j2)为第i2行第j2列的超短基线定位系统信息，nx2为target矩阵总行数，i2、j2、nx2均为自然数。
(1.3)利用GPS数据将超短基线定位系统数据加密处理
将GPS数据与超短基线定位系统数据进行加密处理，根据实际探测的测线有无超短基线定位系统数据的情况，分两种情况进行加密处理。
(1.3.1)若实际探测测线有超短基线定位系统数据，则进行以下操作：
a)若在超短基线定位系统记录时间以外的时间，一般都是综合摄像拖体垂直下放和回收的过程，此时认为综合摄像拖体位于科考船的正下方，将船的GPS信息等同综合摄像拖体的位置信息。若在超短基线定位系统记录到综合摄像拖体位置信息的时间范围内，则计算测线中前100组超短基线定位系统数据的采样周期，取最小的采样周期作为超短基线定位系统的采样周期t。
b)遍历time1(i1)＝ship(i1，1)，计算每个time1(i1)与所有i2对应的超短基线定位系统时间time2(i2)＝target(i2，1)的差值。若差值在超短基线定位系统的采样周期t范围内，即time1(i1)-time2(i2)＜t，则计算time1(i1)时刻对应的母船位置与time2(i2)时刻对应的综合摄像拖体位置之间的距离差，存入矩阵[Lon(i2)，Lat(i2)]，距离差包含经度差Lon(i2)和纬度差Lat(i2)，其中Lon(i2)＝ship(i1，2)-target(i2，2)，Lat(i2)＝ship(i1，3)-target(i2，3)。
c)若在[d1，d2]之间，则将该时刻的Lon(i2)和Lat(i2)存入距离差矩阵LonLat(K，j3)，其中，d1的取值一般为0.005°，d2的取值一般为0.02°，j3取1或2，LonLat(K，1)＝Lon(i2)，LonLat(K，2)＝Lat(i2)，K初值为0，K用来计算time1(i1)时间点符合要求的time2(i2)时间点个数，当time1(i1)时间点对应的在[d1，d2]之间时K加1。
d)当K≥7时，对符合要求的最后7个距离差取平均值，经度差平均值和纬度差平均值分别用mean1和mean2表示，并通过平均值来校正GPS定位信息，形成加密后的超短基线定位系统位置信息，并存入矩阵TARGET(k，l)，其中k＝1，...，n，l取1、2或3分别对应时间、通过平均值校正后的经度、通过平均值校正后的纬度，TARGET(k，l)为第k行第l列的加密后的超短基线定位系统信息，n为TARGET矩阵总行数，l为TARGET矩阵总列数，k、l、n均为自然数；TARGET(k，1)＝ship(i1，1)，TARGET(k，2)＝ship(i1，2)-meanl，TARGET(k，3)＝ship(i1，3)-mean2。
若K＜7，则选取符合要求的最后一个距离差作为校正值，形成加密后的超短基线定位系统位置信息，并存入矩阵TARGET(k，l)，则此时TARGET(k，1)＝ship(i1，1)，TARGET(k，2)＝ship(i1，2)-LonLat(K，1)，TARGET(k，3)＝ship(i1，3)-LonLat(K，2)。
(1.3.2)若实际探测的测线无超短基线定位系统数据，则利用船、综合摄像拖体以及船下方海底三个位置形成的直角三角形推算出综合摄像拖体位置。在综合摄像拖体作业过程中，为保证整个作业的安全，需要将拖缆保持位于船体的竖直对称面上，因此我们认为综合摄像拖体在作业过程中无横向漂移，将作业现场记录的缆长与水深代入计算，利用直角三角形原理估算综合摄像拖体的位置，认为其等效于超短基线定位系统数据，因此将计算结果存入矩阵TARGET(k，l)。
步骤2：匹配浊度传感器探测资料的位置信息
(2.1)读取浊度传感器探测所得的信息MAPR(x，y)(包含探测时间、深度、浊度和温度四列数据)和加密后的超短基线定位系统信息TARGET(k，l)。x＝1，...，s，s为MAPR(x，y)的总行数，y取1、2、3或4分别表示时间、深度、浊度和温度信息的浊度传感器探测数据。默认MAPR矩阵数据的第一列均为探测时间，将其转换为具体的日期值；读取MAPR(x，y)矩阵数据的前100组数据，将数据记录间隔的最小值设为浊度传感器探测数据的采样周期c1。
(2.2)遍历加密后的超短基线定位系统时间序列TARGET(k，1)，计算每个TARGET(k，1)与所有浊度传感器探测数据MAPR(x，1)的时间差值，若时间差的绝对值小于探测数据的采样周期cl，且TARGET(k，1)时刻对应的k在c1/2+1至n-c1/2范围内，则将TARGET(k，1)时刻前半个采样周期取整后的时刻至TARGET(k，1)时刻后半个采样周期取整后的时刻时间范围内的位置信息取平均，作为该时刻的位置信息Loc，将位置信息与浊度传感器探测数据一并存入矩阵MAPR1(x1，y1)，其中，x1＝1，...，s，y1＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、浊度和温度信息。
(2.3)若在同一条测线上有n1(n1≥2)个层位的浊度传感器分别测得探测数据，则取第q(2≤q≤n1)个层位的传感器循环步骤(2.1)和(2.2)，且将循环后的MAPR1(x1，y1)导出并保存至第q个数据文档中。
步骤3：浊度传感器探测资料的处理
(3.1)噪声点处理
选取深度大于1000米的温度和浊度探测值，分别依次进行滤波圆滑、梯度判断和趋势线提取。
滤波圆滑中，滤波窗口大小m设置为12(因为在实际近底拖曳探测过程中，假定小范围(30m)内羽状流的温度和浊度变化较小，船速一般控制在1节左右，因此滤波窗口大小m设置为12比较合适)，滤波处理过程具体为：通过MATLAB软件的nanstd函数计算每一个滤波窗口的标准偏差std，通过MATLAB软件的nanmean函数计算每一个滤波窗口的平均值var；若某一点探测值与对应滤波窗口的平均值var之差是标准偏差std的三倍以上，则认定为噪声点，用平均值var替换该点数值，将计算结果存入矩阵MAPR2(x2，y2)中，x2＝1，...，s，y2＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、滤波圆滑处理后的浊度和滤波圆滑处理后的温度信息。
梯度判断具体如下：根据矩阵MAPR2(x2，y2)数据计算纬度方向(假定测线呈南北走向)和深度方向的变化梯度，分别记为G1和G2，即：
G1(x2)＝(NTU(x2)-NTU(x2-1))/(Lat2(x2)-Lat2(x2-1))；
G2(x2)＝(NTU(x2)-NTU(x2-1))/(Depth(x2)-Depth(x2-1))；
其中，进行温度探测值梯度判断时，NTU(x2)＝MAPR2(x2，6)；进行浊度探测值梯度判断时，NTU(x2)＝MAPR2(x2，5)；
Lat2(x2)＝MAPR2(x2，3)；Depth(x2)＝MAPR2(x2，4)。
根据变化梯度的整体分布情况设定阈值t1和t2，对变化梯度进行判断。若G1(x2)＜t1且G2(x2)＜t2，则保留NTU(x2)原始数据；若G1(x2)＞t1或G2(x2)＞t2，则视为噪声点去除，得到梯度判断后的探测数据矩阵MAPR3(x3，y3)，x3＝1，...，s，y3＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、梯度判断后的浊度和梯度判断后的温度信息。
趋势线提取具体如下：对MAPR1(x1，y1)数据每5个点取一个平均值，连接各个平均值得到最大数据密度包络线，对最大数据密度包络线进行线性插值，获得与MAPR1(x1，y1)数据个数相同的探测值趋势线；将探测值趋势线与梯度处理后的数据进行对比，选取在探测值趋势线上下浮动幅度为M的数据点作为该测线的浊度值或温度值，得到趋势线处理后的数据矩阵MAPR4(x4，y4)，x4＝1，...，s，y4＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、趋势线提取后的浊度和趋势线提取后的温度信息；当探测值为浊度时，M≤0.02NTU；当探测值为温度时，M≤0.05℃。
(3.2)背景海水处理
不同海区海水的成分和性质存在差别，即使在不受热液羽状流影响的深度，所有的浊度传感器仍然记录到一个较小的浊度值，这个值即为当地海水的浊度背景值nb(一般取深度1000-1200m处海水的浊度平均值作为该区域的浊度背景值)。将浊度传感器探测得到的热液羽状流浊度数据扣除海水本身的浊度背景值，得到背景海水处理后的数据矩阵MAPR5(x5，y5)，x5＝1，...，s，y5＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、背景海水处理后的浊度和趋势线提取后的温度信息。MAPR5(x5，y5)真实反映热液异常情况。
(3.3)若在同一条测线上有n1(n1≥2)个层位的浊度传感器分别测得探测数据，则取第q(2≤q≤n1)个层位的传感器循环步骤(3.1)和(3.2)，且将循环后的MAPR5(x5，y5)导出并保存至第q个数据文档中；由于浊度传感器本身性能的问题，即使在假定海水浊度值为定值的海水环境中(深度1000-1200m)作业，不同型号的浊度传感器记录到的浊度探测数据也存在差异，即海水背景值nb的探测结果不同，因此需要对浊度传感器进行系统偏差校正。系统偏差校正具体如下：
在同一条测线上，分别对n1个浊度传感器计算海水取背景值，选取最小的海水背景值作为该条测线的海水背景值标准，用nb_min表示；第n1层的浊度传感器的系统偏差校正值为：
delta(n1)＝nb(n1)-nb_min；
其中，nb(n1)表示第n1层浊度传感器在同一条测线上测得的海水背景值，delta(n1)表示第n1层浊度传感器的海水背景值与该条测线的海水背景值之间的差值。
系统偏差校正后的数据存入数据矩阵MAPR6(x6，y6)中，x6＝1，...，s，y6＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、系统偏差校正后的浊度和趋势线提取后的温度信息。
步骤4：匹配化学传感器探测资料的位置信息
化学传感器探测所得数据chem(u，v)，u＝1，...，nu，nu为Eh探测值总数，v＝1，...，6分别代表探测时间、深度、氧化还原电位Eh、酸碱度pH、溶解氧DO和硫化氢H2S，共6列探测数据。
(4.1)依次读取化学传感器探测数据chem(u，v)和加密后的超短基线定位系统信息TARGET(k，l)；将化学传感器探测数据的探测时间转换为具体的日期值；读取日期值的前100组数据，将最小时间差设为探测数据的采样周期c2。
(4.2)遍历加密后的超短基线定位系统时间序列TARGET(k，l)，计算每个TARGET(k，l)与所有化学传感器探测数据chem(u，v)的时间差值，若两者的时间差绝对值小于采样周期c2，且TARGET(k，1)时刻对应的k在c2/2+1至n-c2/2范围内，则将TARGET(k，l)时刻前半个采样周期取整后的时刻至TARGET(k，l)时刻后半个采样周期取整后的时刻时间范围内的位置信息取平均值，作为该时刻的位置信息Loc1，将位置信息与化学传感器探测数据一并存入chem1(u1，v1)，u1＝1，...，nu，v1＝1，...，8分别代表探测时间、位置信息Loc1中的经度、位置信息Loc1中的纬度、深度、氧化还原电位Eh、酸碱度pH、溶解氧DO和硫化氢H2S。
(4.3)若在同一条测线上有多个层位的化学传感器资料，则循环步骤(4.1)和(4.2)，直至将所有测线的位置信息匹配完成。
步骤5：化学传感器探测资料的处理
(5.1)分别对Eh、pH、DO和H2S四列参数进行滑动平均处理，去除噪声点。
(5.2)计算Eh值的梯度，用gEh(u1)表示。Eh值反应的是热液区氧化还原现象，因此在受热液影响的调查区，会出现Eh值的剧烈变化。可通过计算Eh值的梯度变化，来进一步观察调查区的异常情况，效果会更加明显。若gEh(u1)时刻对应的u1在至范围内，则设定time3(u1)时刻参与计算Eh值梯度的探测数据范围是至时间段对应的数据范围，那么time3(u1)时刻的前点的Eh探测值表示成：Eh(u1_1)，u1_1取值范围为至u1，后点的Eh探测值表示成：Eh(u1_2)，u1_2取值范围为u1至
其中，表示chem1(u1，1)-T/2取整后的时刻，表示chem1(u1，1)+T/2取整后的时刻；fix为取整函数。
据此，可定义每一个Eh探测值对应的梯度值计算公式为：
其中，mean为取平均值运算，Eh(u1)是对chem1(u1，5)的表达；time3(u1)是对chem1(u1，1)的表达。
将计算后的Eh探测值对应的梯度值存入矩阵chem2(u2，v2)，u2＝1，nu1，nuu1为chem2(u2，v2)总行数，v2＝1，...，9分别代表探测时间、位置信息Loc1中的经度、位置信息Loc1中的纬度、深度、滑动平均处理后的氧化还原电位Eh、滑动平均处理后的酸碱度pH、滑动平均处理后的溶解氧DO、滑动平均处理后的硫化氢H2S和Eh探测值对应的梯度值。
步骤6：匹配甲烷传感器探测资料的位置信息
甲烷传感器探测所得数据meth(e，f)，e＝1，...，ne，ne为甲烷传感器探测值总数，f＝1，...，3分别代表探测时间、深度和甲烷；
(6.1)依次读取甲烷传感器探测数据meth(e，f)和加密后的超短基线定位系统信息TARGET(k，l)；将甲烷传感器探测数据的探测时间转换为具体的日期值；读取日期值的前100组数据，将最小时间差设为探测数据的采样周期c3。遍历加密后的超短基线定位系统时间序列TARGET(k，l)，计算每个TARGET(k，l)与所有甲烷传感器探测数据meth(e，f)的时间差值，若两者的时间差绝对值小于采样周期c3，且TARGET(k，1)时刻对应的k在c3/2+1至n-c3/2范围内，则将TARGET(k，l)时刻前半个采样周期取整后的时刻至TARGET(k，l)时刻后半个采样周期取整后的时刻时间范围内的位置信息取平均值，作为该时刻的位置信息Loc2，将位置信息与甲烷传感器探测数据一并存入meth1(e1，f1)，e1＝1，...，ne，f1＝1，...，5分别代表探测时间、位置匹配后的经度、位置匹配后的纬度、深度和甲烷。
(6.2)若在同一条测线上有多个层位的甲烷传感器资料，则循环步骤(6.1)，直至将所有测线的位置信息匹配完成。
步骤7：甲烷传感器探测资料的处理
对甲烷参数进行噪声点处理，得到噪声处理后的甲烷传感器探测数据，存入矩阵meth2(e2，f2)，e2＝1，...，ne，f2＝1，...，5分别代表探测时间、位置匹配后的经度、位置匹配后的纬度、深度和噪声点处理后的甲烷。
步骤8：成图
通过MATLAB和Origin软件对处理好的最终浊度传感器探测数据、化学传感器探测数据和甲烷传感器探测数据绘制浊度-深度剖面图、温度-深度剖面图、Eh-深度剖面图、甲烷-深度剖面图，以及浊度、温度、Eh、甲烷异常综合分布图。根据异常出现的水深位置和平面位置，对比分析热液流体的异常分布情况。
本发明具有的有效效果：
本发明可对多种热液调查过程中采集的多类型外业数据进行融合并快速处理，提高了数据处理的效率，操作方便，费用相对较低，为寻找热液异常提供了更加有效的方法，可在热液调查现场和资料处理中得到广泛应用。
附图说明
图1是本发明的流程图；
图2是本发明实施例中第一条测线近底30m、150m及270m浊度传感器的处理好的最终探测数据中浊度与纬度的关系图(虚线框为异常区域)；
图3是本发明实施例中第二条测线近底150m浊度传感器的处理好的最终探测数据中温度与纬度的关系图；
图4是本发明实施例中第二条测线近底150m浊度传感器的处理好的最终探测数据中浊度与纬度的关系图(虚线框为异常区域)；
图5是本发明实施例中第二条测线近底0m化学传感器的处理好的最终探测数据中Eh值与纬度的关系图(虚线框为异常区域)。
具体实施方式
下面结合附图对本发明作进一步说明。
如图1所示，一种多类型探测传感器综合观测的海底热液探测方法，是指用浊度传感器、化学传感器及甲烷传感器进行综合观测，并对探测所得的资料(包括温度、浊度、氧化还原电位、甲烷探测信息)进行处理、融合并有效成图的方法，包括下列步骤：
步骤1：定位数据准备工作
(1.1)整理GPS定位数据
GPS定位数据，采样周期为1秒，有时间、经度、纬度三列数据。其中原始的GPS定位数据中，时间序列是字符串格式，表现形式为“20170930140500”，即“YYYYMMDDHHMMSS”，分别表示“年份”、“月份”、“日期”、“小时”、“分钟”和“秒”。时间序列是多参数探测数据进行融合的唯一标准，因此在数据准备阶段，通过MATLAB软件的datenum函数将字符串的时间序列转换成具体的日期值。根据实际测线的探测时间，选取该时间范围内的GPS定位数据存入矩阵ship(i1，j1)中，i1＝1，...，nx1，j1＝1，...，nyl，其中nx1为矩阵总行数，nyl为矩阵总列数，i1、j1、nx1、ny1均为自然数。在海上作业现场，往往只需要对刚结束调查的单独一条测线数据进行处理和分析，因此按测线逐条处理即可；而在实验室后期处理中，可将调查航次中其中一个航段的GPS数据进行汇总，统一处理，方便探测数据的位置匹配工作。
(1.2)整理超短基线定位系统(USBL)数据
在海底热液流体探测工作中，以温度、浊度、氧化还原电位、甲烷、硫化氢、pH参数为主要指标，分别在距离综合摄像拖体0m(即加载在综合摄像拖体本体上)、5m、30m、80m、L1、L2的高度加挂综合探测传感器(包括浊度传感器、化学传感器和甲烷传感器)，L1在200-300m中取值，L2在300-500m中取值。当综合摄像拖体近底作业时，超短基线定位系统开始记录综合摄像拖体的位置信息。超短基线定位系统的采样周期不固定，通常在9-11秒之间浮动，有时间、经度、纬度三列数据，与GPS数据的处理过程一样，将超短基线定位系统的时间序列通过MATLAB软件的datenum函数转换成具体的日期值，并将超短基线定位系统数据存入矩阵target(i2，j2)中，其中i2＝1，...，nx2，j2取值为1、2或3分别对应时间、经度、纬度，target(i2，j2)为第i2行第j2列的超短基线定位系统信息，nx2为target矩阵总行数，i2、j2、nx2均为自然数。
(1.3)利用GPS数据将超短基线定位系统数据加密处理
将GPS数据与超短基线定位系统数据进行加密处理，根据实际探测的测线有无超短基线定位系统数据的情况，分两种情况进行加密处理。
(1.3.1)若实际探测测线有超短基线定位系统数据，则进行以下操作：
a)若在超短基线定位系统记录时间以外的时间，一般都是综合摄像拖体垂直下放和回收的过程，此时认为综合摄像拖体位于科考船的正下方，将船的GPS信息等同综合摄像拖体的位置信息。若在超短基线定位系统记录到综合摄像拖体位置信息的时间范围内，则计算测线中前100组超短基线定位系统数据的采样周期，取最小的采样周期作为超短基线定位系统的采样周期t。
b)遍历time1(i1)＝ship(i1，1)，计算每个time1(i1)与所有i2对应的超短基线定位系统时间time2(i2)＝target(i2，1)的差值。若差值在超短基线定位系统的采样周期t范围内，即time1(i1)-time2(i2)＜t，则计算time1(i1)时刻对应的母船位置与time2(i2)时刻对应的综合摄像拖体位置之间的距离差，存入矩阵[Lon(i2)，Lat(i2)]，距离差包含经度差[Lon(i2)和纬度差Lat(i2)，其中Lon(i2)＝ship(i1，2)-target(i2，2)，Lat(i2)＝ship(i1，3)-target(i2，3)。
c)若在[dl，d2]之间，则将该时刻的Lon(i2)和Lat(i2)存入距离差矩阵LonLat(K，j3)，其中，d1的取值一般为0.005°，d2的取值一般为0.02°，j3取1或2，LonLat(K，1)＝Lon(i2)，LonLat(K，2)＝Lat(i2)，K初值为0，K用来计算time1(i1)时间点符合要求的time2(i2)时间点个数，当time1(i1)时间点对应的在[d1，d2]之间时K加1。
d)当K≥7时，对符合要求的最后7个距离差取平均值，经度差平均值和纬度差平均值分别用mean1和mean2表示，并通过平均值来校正GPS定位信息，形成加密后的超短基线定位系统位置信息，并存入矩阵TARGET(k，l)，其中k＝1，...，n，l取1、2或3分别对应时间、通过平均值来校正后的经度、通过平均值来校正后的纬度，TARGET(k，l)为第k行第l列的加密后的超短基线定位系统信息，n为TARGET矩阵总行数，l为TARGET矩阵总列数，k、l、n均为自然数；TARGET(k，1)＝ship(i1，1)，TARGET(k，2)＝ship(i1，2)-mean1，TARGET(k，3)＝ship(i1，3)-mean2。
若K＜7，则选取符合要求的最后一个距离差作为校正值，形成加密后的超短基线定位系统位置信息，并存入矩阵TARGET(k，l)，则此时TARGET(k，1)＝ship(i1，1)，TARGET(k，2)＝ship(i1，2)-LonLat(K，1)，TARGET(k，3)＝ship(i1，3)-LonLat(K，2)。
(1.3.2)若实际探测的测线无超短基线定位系统数据，则利用船、综合摄像拖体以及船下方海底三个位置形成的直角三角形推算出综合摄像拖体位置。在综合摄像拖体作业过程中，为保证整个作业的安全，需要将拖缆保持位于船体的竖直对称面上，因此我们认为综合摄像拖体在作业过程中无横向漂移，将作业现场记录的缆长与水深代入计算，利用直角三角形原理估算综合摄像拖体的位置，认为其等效于超短基线定位系统数据，因此将计算结果存入矩阵TARGET(k，l)。
步骤2：匹配浊度传感器探测资料的位置信息
(2.1)读取浊度传感器探测所得的信息MAPR(x，y)(包含探测时间、深度、浊度和温度四列数据)和加密后的超短基线定位系统信息TARGET(k，l)。x＝1，...，s，s为MAPR(x，y)的总行数，y取1、2、3或4分别表示时间、深度、浊度和温度信息的浊度传感器探测数据。默认MAPR矩阵数据的第一列均为探测时间，将其转换为具体的日期值；读取MAPR(x，y)矩阵数据的前100组数据，将数据记录间隔的最小值设为浊度传感器探测数据的采样周期c1。
(2.2)遍历加密后的超短基线定位系统时间序列TARGET(k，1)，计算每个TARGET(k，1)与所有浊度传感器探测数据MAPR(x，1)的时间差值，若时间差的绝对值小于探测数据的采样周期c1，且TARGET(k，1)时刻对应的k在c1/2+1至n-c1/2范围内，则将TARGET(k，1)时刻前半个采样周期取整后的时刻至TARGET(k，1)时刻后半个采样周期取整后的时刻时间范围内的位置信息取平均，作为该时刻的位置信息Loc，将位置信息与浊度传感器探测数据一并存入矩阵{MAPR1(x1，y1)}，其中，x1＝1，...，s，y1＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、浊度和温度信息。
(2.3)若在同一条测线上有n1(n1≥2)个层位的浊度传感器分别测得探测数据，则取第q(2≤q≤n1)个层位的传感器循环步骤(2.1)和(2.2)，且将循环后的MAPR1(x1，y1)导出并保存至第q个数据文档中。
步骤3：浊度传感器探测资料的处理
(3.1)噪声点处理
选取深度大于1000米的温度和浊度探测值，分别依次进行滤波圆滑、梯度判断和趋势线提取。
滤波圆滑中，滤波窗口大小m设置为12(因为在实际近底拖曳探测过程中，假定小范围(30m)内羽状流的温度和浊度变化较小，船速一般控制在1节左右，因此滤波窗口大小m设置为12比较合适)，滤波处理过程具体为：通过MATLAB软件的nanstd函数计算每一个滤波窗口的标准偏差std，通过MATLAB软件的nanmean函数计算每一个滤波窗口的平均值var；若某一点探测值与对应滤波窗口的平均值var之差是标准偏差std的三倍以上，则认定为噪声点，用平均值var替换该点数值，将计算结果存入矩阵MAPR2(x2，y2)中，x2＝1，...，s，y2＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、滤波圆滑处理后的浊度和滤波圆滑处理后的温度信息。
梯度判断具体如下：根据矩阵MAPR2(x2，y2)数据计算纬度方向(假定测线呈南北走向)和深度方向的变化梯度，分别记为G1和G2，即：
G1(x2)＝(NTU(x2)-NTU(x2-1))/(Lat2(x2)-Lat2(x2-1))；
G2(x2)＝(NTU(x2)-NTU(x2-1))/(Depth(x2)-Depth(x2-1))；
其中，进行温度探测值梯度判断时，NTU(x2)＝MAPR2(x2，6)，进行浊度探测值梯度判断时，NTU(x2)＝MAPR2(x2，5)；
Lat2(x2)＝MAPR2(x2，3)；Depth(x2)＝MAPR2(x2，4)。
根据变化梯度的整体分布情况设定阈值tl和t2，对变化梯度进行判断。若G1(x2)＜t1且G2(x2)＜t2，则保留NTU(x2)原始数据；若G1(x2)＞t1或G2(x2)＞t2，则视为噪声点去除，得到梯度判断后的探测数据矩阵MAPR3(x3，y3)，x3＝1，...，s，y3＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、梯度判断后的浊度和梯度判断后的温度信息。
趋势线提取具体如下：对MAPR1(x1，y1)数据每5个点取一个平均值，连接各个平均值得到最大数据密度包络线，对最大数据密度包络线进行线性插值，获得与MAPR1(x1，y1)数据个数相同的探测值趋势线；将探测值趋势线与梯度处理后的数据进行对比，选取在探测值趋势线上下浮动幅度为M的数据点作为该测线的浊度值或温度值，得到趋势线处理后的数据矩阵MAPR4(x4，y4)，x4＝1，...，s，y4＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、趋势线提取后的浊度和趋势线提取后的温度信息；当探测值为浊度时，M≤0.02NTU；当探测值为温度时，M≤0.05℃。
(3.2)背景海水处理
不同海区海水的成分和性质存在差别，即使在不受热液羽状流影响的深度，所有的浊度传感器仍然记录到一个较小的浊度值，这个值即为当地海水的浊度背景值nb(一般取深度1000-1200m处海水的浊度平均值作为该区域的浊度背景值)。将浊度传感器探测得到的热液羽状流浊度数据扣除海水本身的浊度背景值，得到背景海水处理后的数据矩阵MAPR5(x5，y5)，x5＝1，...，s，y5＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、背景海水处理后的浊度和趋势线提取后的温度信息。MAPR5(x5，y5)真实反映热液异常情况。
(3.3)若在同一条测线上有nl(nl≥2)个层位的浊度传感器分别测得探测数据，则取第q(2≤q≤n1)个层位的传感器循环步骤(3.1)和(3.2)，且将循环后的MAPR5(x5，y5)导出并保存至第q个数据文档中；由于浊度传感器本身性能的问题，即使在假定海水浊度值为定值的海水环境中(深度1000-1200m)作业，不同型号的浊度传感器记录到的浊度探测数据也存在差异，即海水背景值nb的探测结果不同，因此需要对浊度传感器进行系统偏差校正。系统偏差校正具体如下：
在同一条测线上，分别对n1个浊度传感器计算海水取背景值，选取最小的海水背景值作为该条测线的海水背景值标准，用nb_min表示；第nl层的浊度传感器的系统偏差校正值为：
delta(n1)＝nb(n1)-nb_min；
其中，nb(n1)表示第n1层浊度传感器在同一条测线上测得的海水背景值，delta(n1)表示第nl层浊度传感器的海水背景值与该条测线的海水背景值之间的差值。
系统偏差校正后的数据存入数据矩阵MAPR6(x6，y6)中，x6＝1，...，s，y6＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、系统偏差校正后的浊度和趋势线提取后的温度信息。
图2是第一条测线近底30m、150m及270m浊度传感器的处理好的最终探测数据中浊度与纬度的关系图(虚线框为异常区域)；图3是第二条测线近底150m浊度传感器的处理好的最终探测数据中温度与纬度的关系图；图4是第二条测线近底150m浊度传感器的处理好的最终探测数据中浊度与纬度的关系图(虚线框为异常区域)。
步骤4：匹配化学传感器探测资料的位置信息
化学传感器探测所得数据chem(u，v)，u＝1，...，nu，nu为Eh探测值总数，v＝1，...，6分别代表探测时间、深度、氧化还原电位Eh、酸碱度pH、溶解氧DO和硫化氢H2S，共6列探测数据。
(4.1)依次读取化学传感器探测数据chem(u，v)和加密后的超短基线定位系统信息TARGET(k，l)；将化学传感器探测数据的探测时间转换为具体的日期值；读取日期值的前100组数据，将最小时间差设为探测数据的采样周期c2。
(4.2)遍历加密后的超短基线定位系统时间序列TARGET(k，l)，计算每个TARGET(k，l)与所有化学传感器探测数据chem(u，v)的时间差值，若两者的时间差绝对值小于采样周期c2，且TARGET(k，1)时刻对应的k在c2/2+1至n-c2/2范围内，则将TARGET(k，l)时刻前半个采样周期取整后的时刻至TARGET(k，l)时刻后半个采样周期取整后的时刻时间范围内的位置信息取平均值，作为该时刻的位置信息Loc1，将位置信息与化学传感器探测数据一并存入chem1(u1，v1)，u1＝1，...，nu，v1＝1，...，8分别代表探测时间、位置信息Loc1中的经度、位置信息Loc1中的纬度、深度、氧化还原电位Eh、酸碱度pH、溶解氧DO和硫化氢H2S。
(4.3)若在同一条测线上有多个层位的化学传感器资料，则循环步骤(4.1)和(4.2)，直至将所有测线的位置信息匹配完成。
步骤5：化学传感器探测资料的处理
(5.1)分别对Eh、pH、DO和H2S四列参数进行滑动平均处理，去除噪声点。
(5.2)计算Eh值的梯度，用gEh(u1)表示。Eh值反应的是热液区氧化还原现象，因此在受热液影响的调查区，会出现Eh值的剧烈变化。可通过计算Eh值的梯度变化，来进一步观察调查区的异常情况，效果会更加明显。若gEh(u1)时刻对应的u1在至范围内，则设定time3(u1)时刻参与计算Eh值梯度的探测数据范围是至时间段对应的数据范围，那么time3(u1)时刻的前点的Eh探测值表示成：Eh(u1_1)，u1_1取值范围为至u1-1，后点的Eh探测值表示成：Eh(u1_2)，u1_2取值范围为u1+1至其中，表示chem1(u1，1)-T/2取整后的时刻，表示chem1(u1，1)+T/2取整后的时刻；fix为取整函数。
据此，可定义每一个Eh探测值对应的梯度值计算公式为：
其中，mean为取平均值运算，Eh(u1)是对chem1(u1，5)的表达；time3(u1)是对chem1(u1，1)的表达。
将计算后的Eh探测值对应的梯度值存入矩阵chem2(u2，v2)，u2＝1，nu1，nu1为chem2(u2，v2)总行数，v2＝1，...，9分别代表探测时间、位置信息Loc1中的经度、位置信息Loc1中的纬度、深度、滑动平均处理后的氧化还原电位Eh、滑动平均处理后的酸碱度pH、滑动平均处理后的溶解氧DO、滑动平均处理后的硫化氢H2S和Eh探测值对应的梯度值。
图5是第二条测线近底0m化学传感器的处理好的最终探测数据中Eh值与纬度的关系图(虚线框为异常区域)。
步骤6：匹配甲烷传感器探测资料的位置信息
甲烷传感器探测所得数据meth(e，f)，e＝1，...，ne，ne为甲烷传感器探测值总数，f＝1，...，5分别代表探测时间、经度、纬度、深度和甲烷；
(6.1)依次读取甲烷传感器探测数据meth(e，f)和加密后的超短基线定位系统信息TARGET(k，l)；将甲烷传感器探测数据的探测时间转换为具体的日期值；读取日期值的前100组数据，将最小时间差设为探测数据的采样周期c3。遍历加密后的超短基线定位系统时间序列TARGET(k，l)，计算每个TARGET(k，l)与所有甲烷传感器探测数据meth(e，f)的时间差值，若两者的时间差绝对值小于采样周期c3，且TARGET(k，1)时刻对应的k在c3/2+1至n-c3/2范围内，则将TARGET(k，l)时刻前半个采样周期取整后的时刻至TARGET(k， l)时刻后半个采样周期取整后的时刻时间范围内的位置信息取平均值，作为该时刻的位置信息Loc2，将位置信息与甲烷传感器探测数据一并存入meth1(e1，f1)，e1＝1，...，ne，f1＝1，...，5分别代表探测时间、位置匹配后的经度、位置匹配后的纬度、深度和甲烷。
(6.2)若在同一条测线上有多个层位的甲烷传感器资料，则循环步骤(6.1)，直至将所有测线的位置信息匹配完成。
步骤7：甲烷传感器探测资料的处理
对甲烷参数进行噪声点处理，得到噪声处理后的甲烷传感器探测数据，存入矩阵meth2(e2，f2)，e2＝1，...，ne，f2＝1，...，5分别代表探测时间、位置匹配后的经度、位置匹配后的纬度、深度和噪声点处理后的甲烷。
步骤8：成图
通过MATLAB和Origin软件对处理好的最终浊度传感器探测数据、化学传感器探测数据和甲烷传感器探测数据绘制浊度-深度剖面图、温度-深度剖面图、Eh-深度剖面图、甲烷-深度剖面图，以及浊度、温度、Eh、甲烷异常综合分布图。根据异常出现的水深位置和平面位置，对比分析热液流体的异常分布情况。
Claims (2)
1.一种多类型探测传感器综合观测的海底热液探测方法，其特征在于：该方法包括下列步骤：
步骤1：定位数据准备工作
(1.1)整理GPS定位数据
GPS定位数据，采样周期为1秒，有时间、经度、纬度三列数据；通过MATLAB软件的datenum函数将字符串的时间序列转换成具体的日期值；根据实际测线的探测时间，选取该时间范围内的GPS定位数据存入矩阵ship(i1，j1)中，i1＝1，...，nx1，j1＝1，...，ny1，其中nxl为矩阵总行数，nyl为矩阵总列数，i1、j1、nx1、ny1均为自然数；
(1.2)整理超短基线定位系统(USBL)数据
分别在距离综合摄像拖体0m、5m、30m、80m、L1、L2的高度加挂综合探测传感器，综合探测传感器包括浊度传感器、化学传感器和甲烷传感器，L1在200-300m中取值，L2在300-500m中取值；综合摄像拖体下放和回收的整个过程中，综合探测传感器持续探测数据；当综合摄像拖体作业时，超短基线定位系统开始记录综合摄像拖体的位置信息；超短基线定位系统的采样周期在9-11秒之间浮动，有时间、经度、纬度三列数据，与GPS数据的处理过程一样，将超短基线定位系统的时间序列通过MATLAB软件的datenum函数转换成具体的日期值，并将超短基线定位系统数据存入矩阵target(i2，j2)中，其中i2＝1，...，nx2，j2取值为1、2或3分别对应时间、经度、纬度，target(i2，j2)为第i2行第j2列的超短基线定位系统信息，nx2为target矩阵总行数，i2、j2、nx2均为自然数；
(1.3)利用GPS数据将超短基线定位系统数据加密处理
将GPS数据与超短基线定位系统数据进行加密处理，根据实际探测的测线有无超短基线定位系统数据的情况，分两种情况进行加密处理；
(1.3.1)若实际探测测线有超短基线定位系统数据，则进行以下操作：
a)若在超短基线定位系统记录时间以外的时间，认为综合摄像拖体位于科考船的正下方，将船的GPS信息等同综合摄像拖体的位置信息；若在超短基线定位系统记录到综合摄像拖体位置信息的时间范围内，则计算测线中前100组超短基线定位系统数据的采样周期，取最小的采样周期作为超短基线定位系统的采样周期t；
b)遍历time1(i1)＝ship(i1，1)，计算每个time1(i1)与所有i2对应的超短基线定位系统时间time2(i2)＝target(i2，1)的差值；若差值在超短基线定位系统的采样周期t范围内，即time1(i1)-time2(i2)＜t，则计算time1(i1)时刻对应的母船位置与time2(i2)时刻对应的综合摄像拖体位置之间的距离差，存入矩阵[Lon(i2)，Lat(i2)]，距离差包含经度差Lon(i2)和纬度差Lat(i2)，其中Lon(i2)＝ship(i1，2)-target(i2，2)，Lat(i2)＝ship(i1，3)-target(i2，3)；
c)若在[d1，d2]之间，则将该时刻的Lon(i2)和Lat(i2)存入距离差矩阵LonLat(K，j3)，其中，d1的取值一般为0.005°，d2的取值一般为0.02°，j3取1或2，LonLat(K，1)＝Lon(i2)，LonLat(K，2)＝Lat(i2)，K初值为0，K用来计算time1(i1)时间点符合要求的time2(i2)时间点个数，当timel(i1)时间点对应的在[d1，d2]之间时K加1；
d)当K≥7时，对符合要求的最后7个距离差取平均值，经度差平均值和纬度差平均值分别用mean1和mean2表示，并通过平均值来校正GPS定位信息，形成加密后的超短基线定位系统位置信息，并存入矩阵TARGET(k，l)，其中k＝1，...，n，l取1、2或3分别对应时间、通过平均值校正后的经度、通过平均值校正后的纬度，TARGET(k，l)为第k行第l列的加密后的超短基线定位系统信息，n为TARGET矩阵总行数，l为TARGET矩阵总列数，k、l、n均为自然数；TARGET(k，1)＝ship(i1，1)，TARGET(k，2)＝ship(i1，2)-mean1，TARGET(k，3)＝ship(i1，3)-mean2；
若K＜7，则选取符合要求的最后一个距离差作为校正值，形成加密后的超短基线定位系统位置信息，并存入矩阵TARGET(k，l)，则此时TARGET(k，1)＝ship(i1，1)，TARGET(k，2)＝ship(i1，2)-LonLat(K，1)，TARGET(k，3)＝ship(i1，3)-LonLat(K，2)；
(1.3.2)若实际探测的测线无超短基线定位系统数据，则利用船、综合摄像拖体以及船下方海底三个位置形成的直角三角形，将作业现场记录的缆长与水深代入计算，推算出综合摄像拖体位置，认为等效于超短基线定位系统数据，存入矩阵TARGET(k，l)；
步骤2：匹配浊度传感器探测资料的位置信息
(2.1)读取浊度传感器探测所得的信息MAPR(x，y)和加密后的超短基线定位系统信息TARGET(k，l)；MAPR(x，y)包含探测时间、深度、浊度和温度四列数据；x＝1，...，s，s为MAPR(x，y)的总行数，y取1、2、3或4分别表示时间、深度、浊度和温度信息的浊度传感器探测数据；默认MAPR矩阵数据的第一列均为探测时间，将其转换为具体的日期值；读取MAPR(x，y)矩阵数据的前100组数据，将数据记录间隔的最小值设为浊度传感器探测数据的采样周期c1；
(2.2)遍历加密后的超短基线定位系统时间序列TARGET(k，1)，计算每个TARGET(k，1)与所有浊度传感器探测数据MAPR(x，1)的时间差值，若时间差的绝对值小于探测数据的采样周期c1，且TARGET(k，1)时刻对应的k在c1/2+1至n-c1/2范围内，则将TARGET(k，1)时刻前半个采样周期取整后的时刻至TARGET(k，1)时刻后半个采样周期取整后的时刻时间范围内的位置信息取平均，作为该时刻的位置信息Loc，将位置信息与浊度传感器探测数据一并存入矩阵MAPR1(x1，y1)，其中，x1＝1，...，s，y1＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、浊度和温度信息；
(2.3)若在同一条测线上有n1个层位的浊度传感器分别测得探测数据，n1≥2，则依次取第q个层位的传感器循环步骤(2.1)和(2.2)，且将循环后的MAPR1(x1，y1)导出并保存至第q个数据文档中，2≤q≤nl；
步骤3：浊度传感器探测资料的处理
(3.1)噪声点处理
选取深度大于1000米的温度和浊度探测值，分别依次进行滤波圆滑、梯度判断和趋势线提取；
滤波圆滑中，滤波窗口大小m设置为12，滤波处理过程具体为：通过MATLAB软件的nanstd函数计算每一个滤波窗口的标准偏差std，通过MATLAB软件的nanmean函数计算每一个滤波窗口的平均值var；若某一点探测值与对应滤波窗口的平均值var之差是标准偏差std的三倍以上，则认定为噪声点，用平均值var替换该点数值，将计算结果存入矩阵MAPR2(x2，y2)中，x2＝1，...，s，y2＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、滤波圆滑处理后的浊度和滤波圆滑处理后的温度信息；
梯度判断具体如下：根据矩阵MAPR2(x2，y2)数据计算纬度方向和深度方向的变化梯度，分别记为G1和G2，即：
G1(x2)＝(NTU(x2)-NTU(x2-1))/(Lat2(x2)-Lat2(x2-1))；
G2(x2)＝(NTU(x2)-NTU(x2-1))/(Depth(x2)-Depth(x2-1))；
其中，进行温度探测值梯度判断时，NTU(x2)＝MAPR2(x2，6)；进行浊度探测值梯度判断时，NTU(x2)＝MAPR2(x2，5)；
Lat2(x2)＝MAPR2(x2，3)；Depth(x2)＝MAPR2(x2，4)；
根据变化梯度的整体分布情况设定阈值t1和t2，对变化梯度进行判断；若G1(x2)＜t1且G2(x2)＜t2，则保留NTU(x2)原始数据；若G1(x2)＞t1或G2(x2)＞t2，则视为噪声点去除，得到梯度判断后的探测数据矩阵MAPR3(x3，y3)，x3＝1，...，s，y3＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、梯度判断后的浊度和梯度判断后的温度信息；
趋势线提取具体如下：对MAPR1(x1，y1)数据每5个点取一个平均值，连接各个平均值得到最大数据密度包络线，对最大数据密度包络线进行线性插值，获得与MAPR1(x1，y1)数据个数相同的探测值趋势线；将探测值趋势线与梯度处理后的数据进行对比，选取在探测值趋势线上下浮动幅度为M的数据点作为该测线的浊度值或温度值，得到趋势线处理后的数据矩阵MAPR4(x4，y4)，x4＝1，...，s，y4＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、趋势线提取后的浊度和趋势线提取后的温度信息；当探测值为浊度时，M≤0.02NTU；当探测值为温度时，M≤0.05℃；
(3.2)背景海水处理
将浊度传感器探测得到的热液羽状流浊度数据扣除海水本身的浊度背景值，得到背景海水处理后的数据矩阵MAPR5(x5，y5)，x5＝1，...，s，y5＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、背景海水处理后的浊度和趋势线提取后的温度信息；MAPR5(x5，y5)真实反映热液异常情况；
(3.3)若在同一条测线上有n1个层位的浊度传感器分别测得探测数据，n1≥2，则依次取第q个层位的传感器循环步骤(3.1)和(3.2)，且将循环后的MAPR5(x5，y5)导出并保存至第q个数据文档中，2≤q≤n1；对浊度传感器进行系统偏差校正；系统偏差校正具体如下：
在同一条测线上，分别对n1个浊度传感器计算海水取背景值，选取最小的海水背景值作为该条测线的海水背景值标准，用nb_min表示；第n1层的浊度传感器的系统偏差校正值为：
delta(n1)＝nb(n1)-nb_min；
其中，nb(n1)表示第n1层浊度传感器在同一条测线上测得的海水背景值，delta(n1)表示第n1层浊度传感器的海水背景值与该条测线的海水背景值之间的差值；
系统偏差校正后的数据存入数据矩阵MAPR6(x6，y6)中，x6＝1，...，s，y6＝1，...，6分别表示探测时间、位置信息Loc中的经度、位置信息Loc中的纬度、深度、系统偏差校正后的浊度和趋势线提取后的温度信息；
步骤4：匹配化学传感器探测资料的位置信息
化学传感器探测所得数据chem(u，v)，u＝1，...，nu，nu为Eh探测值总数，v＝1，...，6分别代表探测时间、深度、氧化还原电位Eh、酸碱度pH、溶解氧DO和硫化氢H2S，共6列探测数据；
(4.1)依次读取化学传感器探测数据chem(u，v)和加密后的超短基线定位系统信息TARGET(k，l)；将化学传感器探测数据的探测时间转换为具体的日期值；读取日期值的前100组数据，将最小时间差设为探测数据的采样周期c2；
(4.2)遍历加密后的超短基线定位系统时间序列TARGET(k，l)，计算每个TARGET(k，l)与所有化学传感器探测数据chem(u，v)的时间差值，若两者的时间差绝对值小于采样周期c2，且TARGET(k，1)时刻对应的k在c2/2+1至n-c2/2范围内，则将TARGET(k，l)时刻前半个采样周期取整后的时刻至TARGET(k，l)时刻后半个采样周期取整后的时刻时间范围内的位置信息取平均值，作为该时刻的位置信息Loc1，将位置信息与化学传感器探测数据一并存入chem1(u1，v1)，u1＝1，...，nu，v1＝1，...，8分别代表探测时间、位置信息Loc1中的经度、位置信息Loc1中的纬度、深度、氧化还原电位Eh、酸碱度pH、溶解氧DO和硫化氢H2S；
(4.3)若在同一条测线上有多个层位的化学传感器资料，则循环步骤(4.1)和(4.2)，直至将所有测线的位置信息匹配完成；
步骤5：化学传感器探测资料的处理
(5.1)分别对Eh、pH、DO和H2S四列参数进行滑动平均处理，去除噪声点；
(5.2)计算Eh值的梯度，用gEh(u1)表示；若gEh(u1)时刻对应的u1在至范围内，则设定time3(u1)时刻参与计算Eh值梯度的探测数据范围是至时间段对应的数据范围，那么time3(u1)时刻的前点的Eh探测值表示成：Eh(u1_1)，u1_1取值范围为至u1，后点的Eh探测值表示成：Eh(u1_2)，u1_2取值范围为u1至
其中，表示chem1(u1，1)-T/2取整后的时刻，表示chem1(u1，1)+T/2取整后的时刻；fix为取整函数；
定义每一个Eh探测值对应的梯度值计算公式为：
其中，mean为取平均值运算，Eh(u1)是对chem1(u1，5)的表达；time3(u1)是对chem1(u1，1)的表达；
将计算后的Eh探测值对应的梯度值存入矩阵chem2(u2，v2)，u2＝1，nu1，nu1为chem2(u2，v2)总行数，v2＝1，...，9分别代表探测时间、位置信息Loc1中的经度、位置信息Loc1中的纬度、深度、滑动平均处理后的氧化还原电位Eh、滑动平均处理后的酸碱度pH、滑动平均处理后的溶解氧DO、滑动平均处理后的硫化氢H2S和Eh探测值对应的梯度值；
步骤6：匹配甲烷传感器探测资料的位置信息
甲烷传感器探测所得数据meth(e，f)，e＝1，...，ne，ne为甲烷传感器探测值总数，f＝1，...，3分别代表探测时间、深度和甲烷；
(6.1)依次读取甲烷传感器探测数据meth(e，f)和加密后的超短基线定位系统信息TARGET(k，l)；将甲烷传感器探测数据的探测时间转换为具体的日期值；读取日期值的前100组数据，将最小时间差设为探测数据的采样周期c3；遍历加密后的超短基线定位系统时间序列TARGET(k，l)，计算每个TARGET(k，l)与所有甲烷传感器探测数据meth(e，f)的时间差值，若两者的时间差绝对值小于采样周期c3，且TARGET(k，1)时刻对应的k在c3/2+1至n-c3/2范围内，则将TARGET(k，l)时刻前半个采样周期取整后的时刻至TARGET(k，l)时刻后半个采样周期取整后的时刻时间范围内的位置信息取平均值，作为该时刻的位置信息Loc2，将位置信息与甲烷传感器探测数据一并存入meth1(e1，f1)，e1＝1，...，ne，f1＝1，...，5分别代表探测时间、位置匹配后的经度、位置匹配后的纬度、深度和甲烷；
(6.2)若在同一条测线上有多个层位的甲烷传感器资料，则循环步骤(6.1)，直至将所有测线的位置信息匹配完成；
步骤7：甲烷传感器探测资料的处理
对甲烷参数进行噪声点处理，得到噪声处理后的甲烷传感器探测数据，存入矩阵meth2(e2，f2)，e2＝1，...，ne，f2＝1，...，5分别代表探测时间、位置匹配后的经度、位置匹配后的纬度、深度和噪声点处理后的甲烷。
2.根据权利要求1所述的一种多类型探测传感器综合观测的海底热液探测方法，其特征在于：还包括步骤8：成图，具体过程为：
通过MATLAB和Origin软件对处理好的最终浊度传感器探测数据、化学传感器探测数据和甲烷传感器探测数据绘制浊度-深度剖面图、温度-深度剖面图、Eh-深度剖面图、甲烷-深度剖面图，以及浊度、温度、Eh、甲烷异常综合分布图；根据异常出现的水深位置和平面位置，对比分析热液流体的异常分布情况。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810735835.XA CN109101996B (zh) | 2018-07-06 | 2018-07-06 | 一种多类型探测传感器综合观测的海底热液探测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810735835.XA CN109101996B (zh) | 2018-07-06 | 2018-07-06 | 一种多类型探测传感器综合观测的海底热液探测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109101996A true CN109101996A (zh) | 2018-12-28 |
| CN109101996B CN109101996B (zh) | 2021-08-03 |
Family
ID=64845663
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810735835.XA Active CN109101996B (zh) | 2018-07-06 | 2018-07-06 | 一种多类型探测传感器综合观测的海底热液探测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109101996B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115727822A (zh) * | 2022-11-29 | 2023-03-03 | 速度时空信息科技股份有限公司 | 一种海洋监测数据提取用样本分析系统 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1804827A (zh) * | 2006-01-14 | 2006-07-19 | 中国海洋大学 | 海底热液活动探测数据处理和信息管理方法 |
| CN101571599A (zh) * | 2009-06-08 | 2009-11-04 | 浙江大学 | 用于探测深海海底热液硫化物的磁探测系统 |
| CN201434909Y (zh) * | 2009-06-08 | 2010-03-31 | 浙江大学 | 用于探测深海海底热液硫化物的磁探测系统 |
| CN103389077A (zh) * | 2013-07-24 | 2013-11-13 | 国家海洋局第二海洋研究所 | 一种基于mbes的海底沙波地貌运动探测方法 |
| CN103605168A (zh) * | 2013-10-12 | 2014-02-26 | 国家海洋局第二海洋研究所 | 一种海底多金属硫化物综合信息快速找矿方法 |
| CN105510968A (zh) * | 2015-12-31 | 2016-04-20 | 中国海洋大学 | 一种基于地震海洋学的海水物性测量方法 |
| CN106202926A (zh) * | 2016-07-11 | 2016-12-07 | 河南大学 | 基于多节点协同探测的空间系统偏差配准优化方法 |
| CN107219529A (zh) * | 2017-06-07 | 2017-09-29 | 国家深海基地管理中心 | 一种高精度海底地形地貌图的获取方法及系统 |
-
2018
- 2018-07-06 CN CN201810735835.XA patent/CN109101996B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1804827A (zh) * | 2006-01-14 | 2006-07-19 | 中国海洋大学 | 海底热液活动探测数据处理和信息管理方法 |
| CN101571599A (zh) * | 2009-06-08 | 2009-11-04 | 浙江大学 | 用于探测深海海底热液硫化物的磁探测系统 |
| CN201434909Y (zh) * | 2009-06-08 | 2010-03-31 | 浙江大学 | 用于探测深海海底热液硫化物的磁探测系统 |
| CN103389077A (zh) * | 2013-07-24 | 2013-11-13 | 国家海洋局第二海洋研究所 | 一种基于mbes的海底沙波地貌运动探测方法 |
| CN103605168A (zh) * | 2013-10-12 | 2014-02-26 | 国家海洋局第二海洋研究所 | 一种海底多金属硫化物综合信息快速找矿方法 |
| CN105510968A (zh) * | 2015-12-31 | 2016-04-20 | 中国海洋大学 | 一种基于地震海洋学的海水物性测量方法 |
| CN106202926A (zh) * | 2016-07-11 | 2016-12-07 | 河南大学 | 基于多节点协同探测的空间系统偏差配准优化方法 |
| CN107219529A (zh) * | 2017-06-07 | 2017-09-29 | 国家深海基地管理中心 | 一种高精度海底地形地貌图的获取方法及系统 |
Non-Patent Citations (5)
| Title |
|---|
| 潘依雯等: "基于多参数化学传感器的海底热液探测方法研究", 《海洋学报》 * |
| 秦华伟等: "化学传感器链的集成及其在龟山岛海底热液区的实际应用", 《应用技术》 * |
| 秦华伟等: "热液喷口探测化学传感器的研制及应用", 《热带海洋学报》 * |
| 秦华伟等: "舟山海域垂直方向上水动力特征数值模拟分析", 《杭州电子科技大学学报》 * |
| 苏程等: "实时多波束海底地形数据三维瀑布图绘制方法", 《计算机应用与软件》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115727822A (zh) * | 2022-11-29 | 2023-03-03 | 速度时空信息科技股份有限公司 | 一种海洋监测数据提取用样本分析系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109101996B (zh) | 2021-08-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105865424B (zh) | 一种基于非线性模型的多光谱遥感水深反演方法及装置 | |
| Lipizer et al. | Qualified temperature, salinity and dissolved oxygen climatologies in a changing Adriatic Sea | |
| Dinehart et al. | Repeated surveys by acoustic Doppler current profiler for flow and sediment dynamics in a tidal river | |
| CN104180873B (zh) | 一种单波束测深仪水深粗差检测修正方法及系统 | |
| CN109425906B (zh) | 一种磁异常探测矢量磁目标识别方法 | |
| CN106052795B (zh) | 一种获取潮位的方法及装置 | |
| Zhou et al. | Intrusion pattern of the Kuroshio Subsurface Water onto the East China Sea continental shelf traced by dissolved inorganic iodine species during the spring and autumn of 2014 | |
| Wu et al. | Research on the relative positions-constrained pattern matching method for underwater gravity-aided inertial navigation | |
| CN112819249A (zh) | 基于走航adcp观测海流数据的潮流调和分析计算方法 | |
| CN103513286A (zh) | 一种地质模型约束下的滩坝结构单元判别方法 | |
| Fuhrmann et al. | High-resolution determination of the pH of seawater with a flow-through system | |
| CN108959741A (zh) | 一种基于海洋物理生态耦合模型的参数优化方法 | |
| CN111366936B (zh) | 多波束测深数据处理方法及装置 | |
| Kim et al. | Changes in seawater N: P ratios in the northwestern Pacific Ocean in response to increasing atmospheric N deposition: Results from the East (Japan) Sea | |
| CN109101996A (zh) | 一种多类型探测传感器综合观测的海底热液探测方法 | |
| Keiner et al. | Empirical orthogonal function analysis of sea surface temperature patterns in Delaware Bay | |
| Leong et al. | Monitoring harmful algal blooms in Singapore: Developing a HABs observing system | |
| Pascual | Combining new and conventional sensors to study the Balearic Current | |
| CN107407740A (zh) | 用于远程测量和量化来自海洋铁富集的二氧化碳封存的过程和方法 | |
| Ding et al. | The optimization of Ag/Ag2S electrode using carrier electroplating of nano silver particles and its preliminary application to offshore Kueishan Tao, Taiwan | |
| Jones et al. | A computer-based system for the acquisition and display of continuous vertical profiles of temperature, salinity, fluorescence and nutrients | |
| CN110617803A (zh) | 一种用于水下地形测量中的水深检校方法 | |
| CN107817332A (zh) | 一种表征潜在烃源岩后期保存效率的方法 | |
| Al Balushi | Detecting Submarine Groundwater Plumes in the Salalah Area: A Radon Tracer and Numerical Integral Plume Model Approach | |
| CN113836870B (zh) | 陆棚沉积亚相确定方法及装置 |
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 |