CN111931555B - 一种利用视频图像识别船舶ais是否开启的方法 - Google Patents

一种利用视频图像识别船舶ais是否开启的方法 Download PDF

Info

Publication number
CN111931555B
CN111931555B CN202010539454.1A CN202010539454A CN111931555B CN 111931555 B CN111931555 B CN 111931555B CN 202010539454 A CN202010539454 A CN 202010539454A CN 111931555 B CN111931555 B CN 111931555B
Authority
CN
China
Prior art keywords
ship
ais
target
image
coordinates
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.)
Active
Application number
CN202010539454.1A
Other languages
English (en)
Other versions
CN111931555A (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.)
China Academy of Transportation Sciences
Original Assignee
China Academy of Transportation Sciences
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 China Academy of Transportation Sciences filed Critical China Academy of Transportation Sciences
Priority to CN202010539454.1A priority Critical patent/CN111931555B/zh
Publication of CN111931555A publication Critical patent/CN111931555A/zh
Application granted granted Critical
Publication of CN111931555B publication Critical patent/CN111931555B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/40—Scenes; Scene-specific elements in video content
    • G06V20/41—Higher-level, semantic clustering, classification or understanding of video scenes, e.g. detection, labelling or Markovian modelling of sport events or news items
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • 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/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
    • G06V10/267—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion by performing operations on regions, e.g. growing, shrinking or watersheds
    • 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/40—Extraction of image or video features
    • G06V10/44—Local feature extraction by analysis of parts of the pattern, e.g. by detecting edges, contours, loops, corners, strokes or intersections; Connectivity analysis, e.g. of connected components
    • 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/40—Extraction of image or video features
    • G06V10/50—Extraction of image or video features by performing operations within image blocks; by using histograms, e.g. histogram of oriented gradients [HoG]; by summing image-intensity values; Projection analysis
    • G06V10/507—Summing image-intensity values; Histogram projection analysis
    • 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/40—Extraction of image or video features
    • G06V10/56—Extraction of image or video features relating to colour
    • 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/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/75—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
    • G06V10/751—Comparing pixel values or logical combinations thereof, or feature values having positional relevance, e.g. template matching
    • 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/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/467—Encoded features or binary features, e.g. local binary patterns [LBP]

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • Computational Linguistics (AREA)
  • Databases & Information Systems (AREA)
  • Computing Systems (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • General Engineering & Computer Science (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

本发明涉及一种利用视频图像识别船舶AIS是否开启的方法，包括以下步骤：S1、船舶视频图像识别步骤；S2、船舶AIS轨迹分析步骤；S3、船舶视频图像与船舶AIS轨迹融合步骤，判定未开启AIS的船舶。本发明的优点是：解决了海事执法中对于AIS违法关闭情形的判定依赖于人工比对和登船核实的问题。尤其是在船舶流量密集的航道，全部登船检查或者准确人工比对难度极大，严重影响执法效率，会对船舶航行安全造成巨大的隐患。

Description

一种利用视频图像识别船舶AIS是否开启的方法
技术领域
本发明涉及一种利用视频图像识别船舶AIS是否开启的方法，属于海事监管技术领域。
背景技术
AIS(Automatic identification System)又称船舶自动识别系统，船舶在航行时通过AIS的特定通信频段自动向周围发送并接受周围船舶的名称、位置、航速、航向等动静态信息，是船舶了解周边交通情况避免发生碰撞的关键设备，同时也是海事交管人员掌握辖区船舶交通情况，组织指挥船舶航行的主要参考依据。然而在实际情况中，船舶驾驶人员为了规避海事监管或者疏于设备维护，经常造成船舶在航行中AIS未开启的情况发生，给自身及围边船舶航行安全造成了极大的风险。目前，海事监管人员对船舶AIS开启状态只能在船舶靠泊时通过登船检查的方式或者通过AIS岸台、CCTV监控视频的人工比对方式进行核查。
一方面现有核查手段只能发现AIS设备因为损坏长期无法开启的情况，无法及时发现船舶在航行中主动关闭AIS设备的情况，监管效率较低；另一方面，在船舶流量密集的航道，全部登船检查或者准确人工比对难度极大，严重影响执法效率，会对船舶航行安全造成巨大的隐患。尽管海事监管中已经应用了AIS岸台、CCTV监控视频等多种信息采集手段，但是缺少自动将AIS信息与CCTV视频图像进行关联比较并识别出船舶AIS是否开启的方法。
发明内容
为克服现有技术的缺陷，本方案的目的是利用AIS岸台采集的船舶动、静态信息与CCTV视频服务器采集的船舶图像信息，通过视频图像识别结合船舶AIS轨迹数据，发明一种通过视频自动识别船舶AIS是否开启的方法，识别出船舶AIS的开启状态，避免其在未开启AIS的情况下进入交通流量较大的航段航行，能够极大的提高海事监管能力，保障航行安全避免造成重大损失，具有十分重要的现实意义。本发明提供一种利用视频图像识别船舶AIS是否开启的方法，本发明的技术方案是：
一种利用视频图像识别船舶AIS是否开启的方法，包括以下步骤：
S1、船舶视频图像识别步骤；
S2、船舶AIS轨迹分析步骤；
S3、船舶视频图像与船舶AIS轨迹融合步骤，判定未开启AIS的船舶。
所述的步骤S1为：利用CCTV视频服务器，对船舶进行前景提取和多特征关联，生成已识别待匹配的视频船舶目标集合，每个视频船舶目标进入视场时间、离开视场时间、位置坐标集、船舶长度、船舶速度和航行方向六个特征向量。
所述的步骤S2为：利用船舶AIS报文位置信息，通过对位置信息的坐标转换、航迹插值得到待匹配的AIS船舶目标集合以及每个AIS船舶目标进入视场时间、离开视场时间、位置坐标集、船舶长度、船舶速度和航行方向六个特征向量。
所述的步骤S3为：融合步骤S1和步骤S2获得的待匹配的视频船舶目标集合与待匹配的AIS船舶目标集合，在相同的坐标体系之下，对步骤S1以及步骤S2船舶集合的六个特征向量进行匹配，最后得到匹配船舶集合与未匹配船舶集合，其中，视频船舶目标集合中属于未匹配集合中的就是未开启AIS的船舶。
所述的步骤S1具体为：利用ViBe运动检测算法检测视频图像背景，包括以下步骤，
步骤S1-1、背景建模：对每个位置的像素点N上有代表性的颜色值vi进行合成，即背景模型Mo(x,y)＝{v1,v2,…,vN}，通常设定N＝5；
步骤S1-2、摄像机模型建立：采用带有一阶径向畸变的针孔摄像机模型：
(1)对于世界坐标系Ow-XwYwZw及其对应的摄像机坐标系Oc-XcYcZc与图像坐标系Ou-xuyu，世界坐标系下空间中任意一点P(Xw，Yw，Zw)，Pc(xc，yc)是其相应的摄像机坐标点，Pu(xu，yu)是其相应的理想图像坐标点，即针孔成像坐标，Pd(xd，yd)为畸变点图像坐标，只考虑摄像机的一阶径向畸变，无论畸变程度多少，从图像中心点Oc到图像畸变点Pd的向量OuPd的方向保持不变，而且与世界空间坐标系中心点Poz到P的向量PozP平行；由该条件可得：
其中，摄像机的等效焦距f对图像的x轴和y轴产生同样的作用，即对畸变点xd和yd的影响相同，所以焦距的大小不影响向量OuPd的方向；世界坐标系原点的定义位置与向量OuPd的方向无关；
步骤S1-3、标定摄像机参数：全部的摄像机参数，包括摄像机有效焦距f、平移矩阵T、旋转矩阵以及一阶畸变系数k，该摄像机参数标定的方法先用线性方法进行线性部分的求解，计算出摄像机旋转矩阵R参数的初始值，然后对包括摄像机有效焦距f、平移矩阵T、以及一阶畸变系数k进行非线性优化；
步骤1-3-1、线性部分的求解：
根据世界坐标和摄像机坐标之间的转化关系，设Zw＝0，则有：
其中，tx、ty为平移参数，对于将(2)两侧同时除以ty并整理可以得到：
上式中，令a1＝r1/ty，a2＝r2/ty，a3＝r3/ty，a4＝r4/ty，a5＝r5/Ty，则上式有五个未知数，对于每一个参考点通过上式得到一个方程，当有超过五个参考点时，采用最小二乘法得到上面五个未知数的解；定义一个矩阵如果该矩阵满秩矩阵，则有
否则的话，有
其中，而ai，aj是剩余不为零的两个参数，有：r1＝a1ty r2＝a2ty r4＝a4ty r5＝a5ty tx＝a3ty，根据所得外参数r1，r2，r4，r5，tx的值和已知世界坐标和图像坐标的对应点，估算出该点在摄像机坐标系下的坐标(Xc，Yc)；由于旋转矩阵R具有单位正交性质，则其余外参数计算方法如下：
则有：
r7＝r2r6-r3r5 r8＝r3r4-r1r6 r9＝r1r5-r2r4；
步骤1-3-2、非线性部分的求解：
基于上述计算的参数，根据一阶径向畸变摄像机模型得到如下计算公式：
先设定畸变系数为零，整理上式得到：
其中，F＝r4Xw+r5Yw+ty，G＝r7Xw+r8Yw，选取数个标定点，则上式形成一个超定线性方程组，利用线性最小二乘法求出f和tz，然后将求出的f、tz以及k1＝0作为初始值进一步进行非线性优化计算出f、tz以及k1的真实值；上述标定方法中主点坐标被认为是图像中心点坐标，图像的纵横比sx＝1；利用同一摄像机不同焦距拍摄同一场景得到主点坐标；计算如下：
其中，(u1,v1)和(u2,v2)分别为焦距1和焦距2下的同一个特征点的坐标，选择多的特征点及可求出主点(cx,cy)的坐标；由于摄像机焦距在图像x方向和y方向同时缩放，所以垂直拍摄一个规则物体分别计算图像水平和垂直方向上的像素直径比就可以计算图像的横纵比sx；
利用张正友模板标定法得到主点坐标和横纵比，通过以上步骤，求出摄像机的全部参数，然后导入5个以上图像位置像素坐标点及其在北京54坐标系下对应时刻的船舶AIS坐标点，就可以获得摄像机投影矩阵的全部参数，即摄像机有效焦距f、平移矩阵T、旋转矩阵R以及一阶畸变系数k，根据该投影矩阵，将摄像机图像坐标下的船舶目标的轨迹位置特征值进行转换。
步骤1-4、前景物体分割：
将当前像素值vt(x)与相应背景模型中的每个像素值进行比较，如果最少有1个像素值与当前像素值的距离小于给定阈值γ，其中，γ＝10，那么当前像素点被判断为背景，反之判断为像素值波动；
步骤1-5、背景模型更新：
如果一个像素点被判断为背景，则该像素点随机地替代自己背景模型中的一个像素值；
步骤1-6、前景去噪：该前景去噪包括腐蚀运算和膨胀运算的步骤：
所述的腐蚀运算的步骤包括：
E(f(x,y))＝{min(f(xi,yi)|(xi,yi)∈f(x,y)∩B} (10)
选择算子覆盖下的二值图像的最小像素值作为中心像素的大小，腐蚀运算在该处的作用是去除小面积噪声，B为结构算子，即参与计算的临近像素集；
所述的膨胀运算的步骤包括：
D(f(x,y))＝{max(f(xi,yi)|(xi,yi)∈f(x,y)∩B} (11)
选择算子覆盖下的二值图像的最大像素值作为中心像素的大小，膨胀运算在该处的作用是填补前景区域内小的空缺，B为结构算子，即参与计算的临近像素集。
所述的步骤S2具体为分离目标的合并和多轨迹关联：
步骤2-1、分离目标的合并：
在无遮挡情况下，每一个视频帧检测出的船舶目标相邻，通过聚类方法实现轨迹的连通，具体如下：连通域是由一个最小外接矩形来标记的，矩形框的表示方法为：
[x y width height] (12)
x,y表示的是矩形框的左上角的坐标，横轴x，纵轴y，width表示矩形框的宽度，即横轴长，height表示矩形框的高度，即纵轴长；
当两个矩形框满足公式(13)、(14)的条件时，默认这其中的两个前景块是属于同一个目标的将它们合并成一个大的矩形框；
|xi-xj|＜＝max(widthi,widthj)+width_threshold (13)
|yi-yj|＜＝max(heighti,heightj)+height_threshold (14)
对于合并后的框[X Y widthi,j heighti,j]，如果满足公式(15)、(16)对大小限制、公式(17)、(18)对景深限制、公式(19)对比例限制、公式(11)对有效前景比例中的某任一种限制，则取消合并：
步骤2-1-1、满足大小限制：
合并后的框在高度上不会超过图像高度的五分之二，在宽度上也不会超过图像宽度的五分之一，同时，面积也不会超过两万像素，默认标准图像帧大小归一化为288*352像素，具体是：
heighti,j＞(frame_height)/2.5 (15)
widthi,j＞(frame_width)/5 (16)
其中，[frame_width frame_height]微视频图像的宽度和高度；
步骤2-1-2、满足景深限制：
合并后的位置如果在图像的偏上部分，即处于远景范围内，此时如果面积过大，在真实场景中可能是水面经过带来的暂时性大面积前景或者是前景检测失效，显然不是船舶目标，具体是：
widthi,j*heighti,j＞20000 (17)
Y＜200&widthi,j*heighti,j＞15000 (18)
步骤2-1-3、满足比例限制：
由于船的形状有固定特点，船舶目标的长宽比在预估范围内，经验值是高宽比不会大于6或者小于1，具体是：
heighti,j/widthi,j＞6 or heighti,j/widthi,j＜1 (19)
步骤2-2、多轨迹关联：
对于合并目标形成的轨迹，如果中间出现遮挡或者丢失检测，那么重新出现的目标被认为是新的目标轨迹，造成轨迹完整性缺失，结合目标的内、外部的纹理特征、颜色特征、外部特征进行多轨迹关联：
步骤2-2-1、计算纹理特征LBP：
(1)找到前景图的各个连通域后，按照最小外接矩形的大小在原图上分割出样本，记为l1,对样本l1进行灰度化和大小变换，统一变换到32*16的大小，记为l2；
(2)用大小为3*3的模板遍历样本l2，对模板中心的像素，比较周围8个邻域像素和中心像素的大小，大于中心像素的取1，否则取0，求取LBP特征值，记为LBP1；
(3)对该特征值的二进制数进行循环移位，找到最小的那个值，为该中心像素的旋转不变LBP值，记为LBP2；
(4)模板遍历完样本l2之后，对所有样本的LBP2进行直方图统计，因为是旋转不变8邻域的LBP特征，所以是36维，在36维上统计直方图，得到H1；
(5)不同连通域m、n之间的灰度直方图距离为：
LBP(m,n)＝|Hm-LBP-Hn-LBP| (20)
步骤2-2-2、计算外部特征：
thre＝(α·height+βv)·Tdisappear (21)
thre表示对每一个连通区域设置不同的“近邻”阈值，α、β是高度和速度的权重，height是目标的高度，v是目标速度，目标最近10帧的平均路程计算所得：
n表示当前连通区域内的目标点个数，ft＝(xt,yt)序列点的中心位置。Tdisappear表示连通区域丢失更新的帧数；
则不同连通域m、n之间的外部特征距离为：
Hist(m,n)＝|threm-traj-thren-traj| (23)
步骤2-2-3、融合更新连通域：
循环比较连通区域集合Traj内各个连通区域之间的特征，进行融合，更新得到新的船舶目标集合Traj：
(1)对于每一个新生成的连通区域traji，记录进入视场时间t′i，搜索其近邻值threi内的k个前景域[blob1,blob2...blobk]，计算与每个前景域的特征匹配度：
Score＝pLBP*LBP(traji,blobj)+pHist*Hist(traji,blobj) (24)
其中，pLBP、pHist分别为加权因子，这里两者均取0.5；
(2)找到匹配度最高，即Score最小的前景域，分配给当前连通域，并更新该连通域的参数height、v；
(3)如果未找到匹配的前景域，则当前连通域traji已达到最大，记录其离开视场时间t″i；
(4)更新，重复融合更新连通域的步骤(1)至步骤(4)，直到找不到可以连通的区域。
所述的步骤S3具体为：
将北京54坐标系的平面直角笛卡尔坐标系，作为进行位置关联的绝对坐标系，CCTV视频图像的识别基于视频图像的相对坐标，AIS船舶位置报告通过GPS所测的经纬度来确定目标位置，发送的目标位置信息是基于WGS-84地理坐标系的经度和纬度；为了使两者进行空间位置的关联，将两种坐标系统一到基于北京54坐标系的平面直角笛卡尔坐标系下的坐标；具体为：
步骤3-1、轨迹转换：
(1)由W84下的大地坐标系(L84，B84，H84)转换为W84下的平面直角坐标(X84，Y84，Z84)，其转换方法为：
其中，W84坐标系下的第一偏心率：卯酉圈曲率半径椭球体长半轴：a1＝6378137m；椭球体短半轴：b1＝6356752.3142m；
(2)采用七参数法将W84下的平面直角坐标(X84，Y84，Z84)转换到北京54体系下的直角坐标(X54，Y54，Z54)，七参数方法的计算公式如下：
七参数法除了考虑三个坐标轴方向的平移(ΔX，ΔY，ΔZ)，也考虑了旋转(α，β，γ)和缩放比例k，ΔX，ΔY，ΔZ，α，β，γ这些参数的数值通过已知的WGS-84地理坐标系和北京54坐标系下的控制点获得；
(3)由北京54体系下的平面直角坐标(X54，Y54，Z54)进行椭球转换为北京54体系下的大地坐标(L54，B54，H54)：
其中，北京54体系下的第一偏心率：卯酉圈曲率半径
(4)由北京54体系下的大地坐标(L54，B54，H54)通过高斯-克吕格投影转换为北京54体系下的平面直角坐标(X54，Y54)，高斯投影的性质是投影后角度不变、长度比与点位有关，与方向无关、离中央子午线越远变形越大，高斯投影必须满足中央子午线投影后是直线且长度不变原则，同时必须满足正形投影条件；基于6度带分带的高斯投影的正算公式如下：
其中，x，y为高斯平面直角坐标系的横纵坐标值；L是从中央经线起算的经度值，B是纬度值，需要注意的是二者是以弧度形式表示的；X是赤道到纬度弧度为B处的子午线的长度；η为地球的第二偏心率，N为纬度弧度为B处的卯酉圈曲率半径；
步骤3-2、非线性插值：
在船舶航行中，AIS的报文发送频率随着船舶运行状态动态变化，对AIS航迹进行插值处理，由于船舶在实际航行中需要随航道走向随时调整航速航向，其船舶航迹近似于曲线，并且可以获取船舶在各个位置的航向，因此采用三次样条插值对船舶航迹进行插值拟合，以满足数据融合的条件；
设AIS报文在一个时间段[a，b]内的航迹进行插值，对于AIS航迹在平面直角坐标系上的投影的横坐标x0≤x1...≤xi≤...≤xn，有yi＝∫(xi)，则可以构造一个三次样条插值函数，使它满足3个条件：
(1)yi＝S(xi)，a≤i≤b；
(2)在每一个区间[xi，xi+1]内，S(x)是三次多项式；
(3)在[a，b]内S(x)具有二阶连续导数；
计算出被插值的函数在节点处的导数值，S′(xi)＝mi，则根据条件yi＝S(xi)，yi+1＝S(xi+1)，S′(xi)＝mi，S′(xi+1)＝mi+1，可以使用Hermmite插值公式，计算出每一个小区间[xi，xi+1]内的三次样条函数S(x)的计算公式：
其中hi＝xi+1-xi，对上式的左右两边求二阶导数可得：
因为S(x)在每一个插值点上同样具有二阶连续导数，则可以认为因此有：
根据上式，引入αi，βi，设
由此，对每一个采样点建立一个方程，则有如下方程：
(1-αi)mi-1+2mi+αimi＝βi (36)
上式中，对应n个采样点，得到包含n+1未知变量m0,m1,...,mn的n+1个线性方程组，因为未知变量个数大于方程个数，因此增加一个附加自然边界条件，即被插值函数在两个端点处二阶导数为零，设可得：
其中，α0＝1，α1＝1，至此，方程联立后求出唯一解m0,m1,...,mn，将解出的xi，yim0,m1,...,mn代入上式就可以得到三次样条插值函数；
步骤3-3、初次关联：
对于视频船舶目标traji，以其进入视场的时间ti′为中心,在初级时间关联阈值内进行AIS的数据初级匹配，满足初级关联判断的条件如下示：
kaj-kri＜PATT j＝1,2...,N i＝1,2,...,M (38)
Raj-Rri＜PADT j＝1,2...,N i＝1,2,...,M (39)
其中，kri表示视频图像目标i出现在视场内的时刻，kaj表示AIS目标j出现在视场内的时刻，Raj表示AIS目标j的位置，Rri表示视频图像目标i的位置，PATT表示初级时间关联阈值，PADT表示距离处级关联阈值，设置为穿越视场在航道最大投影的长度；如果没有满足上式的AIS目标船舶，则视频图像内船舶目标均为意思未开启AIS船舶；当上式满足时，认为AIS目标j和视频图像对象i在同一个时间范畴内满足初级关联条件，作为初级匹配的候选对象；
步骤3-4、二级关联：
对待匹配的目标集合进行二级匹配，通过各个特征值的匹配程度加权值进行排序，对排序优先的船舶纳入匹配集合，得到未匹配的视频船舶目标集合与未匹配的AIS船舶目标集合，其中未匹配的视频船舶目标集合即为未开启AIS的船舶目标的；AIS与视频图像中船舶数据关联的结果可以分为两个等级，形成关联结果评价集合W＝{1,0}，其中，1代表关联，0代表不关联；在直积F×W上定义从F到W的模糊关系矩阵就是因素评判矩阵F：
其中，f11ij(k)、f12ij(k)、f13ij(k)、f14ij(k)分别表示位置集合、长度、方向、速度四个单因素集合g1ij(k)、g2ij(k)、g3ij(k)、g4ij(k)在k时刻的关联隶属度，f11ij(k)、f12ij(k)、f13ij(k)、f14ij(k)分别表示单因素集合g1ij(k)、g2ij(k)、g3ij(k)、g4ij(k)在k时刻的不关联隶属度；
设A＝[a1,a2,a3,a4]为位置、长度、方向、速度四个单因素集合关联隶属度的权重集，视频图像船舶目标与AIS目标关联的综合因素关联隶属度可以由单因素权重模糊集A与F的合成运算结果获得；
合成运算的模型包括加权平均型、主观因素决定型以及混合型等。这里采用简单且易于工程实现的加权平均型，有
Cij(k)＝A×F (41)
求出与第i个视频图像目标满足初级关联条件第j个AIS目标的关联隶属度的值，其中隶属度最大的目标，即为匹配的目标：
Max(Cij(k))＝Max(Ci1(k),Ci2(k),...,CiN(k)),(j＝1,2,...,M) (42)
最大值对应的AIS目标即为关联度最大的船只；
当完成视频图像船舶集合与AIS船舶集合的匹配确定重复船舶之后，得到视频图像集合中未匹配船舶集合Traj∩(CU(Traj∩R))，即为未开启AIS的疑似船舶。
本发明的优点是：解决了海事执法中对于AIS违法关闭情形的判定依赖于人工比对和登船核实的问题。尤其是在船舶流量密集的航道，全部登船检查或者准确人工比对难度极大，严重影响执法效率，会对船舶航行安全造成巨大的隐患。
同时由于采集视频的CCTV的视频服务器是海事监管的现有设备，而AIS船台所有航行船舶都强制安装，因此该技术不需要额外的硬件设备成本，通过该方法直接提高了海事对于AIS违规关闭等情况的监管执法效率，节省了人力物力，间接的也提高了辖区水上交通安全的整体水平，经济性和推广前景良好。
具体实施方式
下面结合具体实施例来进一步描述本发明，本发明的优点和特点将会随着描述而更为清楚。但这些实施例仅是范例性的，并不对本发明的范围构成任何限制。本领域技术人员应该理解的是，在不偏离本发明的精神和范围下可以对本发明技术方案的细节和形式进行修改或替换，但这些修改和替换均落入本发明的保护范围内。
本发明涉及一种一种利用视频图像识别船舶AIS是否开启的方法，包括以下步骤：
S1、船舶视频图像识别步骤；
S2、船舶AIS轨迹分析步骤；
S3、船舶视频图像与船舶AIS轨迹融合步骤，判定未开启AIS的船舶。
所述的步骤S1为：利用CCTV视频服务器，对船舶进行前景提取和多特征关联，生成已识别待匹配的视频船舶目标集合，每个视频船舶目标进入视场时间、离开视场时间、位置坐标集、船舶长度、船舶速度和航行方向六个特征向量。
所述的步骤S2为：利用船舶AIS报文位置信息，通过对位置信息的坐标转换、航迹插值得到待匹配的AIS船舶目标集合以及每个AIS船舶目标进入视场时间、离开视场时间、位置坐标集、船舶长度、船舶速度和航行方向六个特征向量。
所述的步骤S3为：融合步骤S1和步骤S2获得的待匹配的视频船舶目标集合与待匹配的AIS船舶目标集合，在相同的坐标体系之下，对步骤S1以及步骤S2船舶集合的六个特征向量进行匹配，最后得到匹配船舶集合与未匹配船舶集合，其中，视频船舶目标集合中属于未匹配集合中的就是未开启AIS的船舶。
所述的步骤S1具体为：利用ViBe运动检测算法检测视频图像背景，包括以下步骤，
步骤S1-1、背景建模：对每个位置的像素点N上有代表性的颜色值vi进行合成，即背景模型Mo(x,y)＝{v1,v2,…,vN}，通常设定N＝5；
步骤S1-2、摄像机模型建立：采用带有一阶径向畸变的针孔摄像机模型：
(1)对于世界坐标系Ow-XwYwZw及其对应的摄像机坐标系Oc-XcYcZc与图像坐标系Ou-xuyu，世界坐标系下空间中任意一点P(Xw，Yw，Zw)，Pc(xc，yc)是其相应的摄像机坐标点，Pu(xu，yu)是其相应的理想图像坐标点，即针孔成像坐标，Pd(xd，yd)为畸变点图像坐标，只考虑摄像机的一阶径向畸变，无论畸变程度多少，从图像中心点Oc到图像畸变点Pd的向量OuPd的方向保持不变，而且与世界空间坐标系中心点Poz到P的向量PozP平行；由该条件可得：
其中，摄像机的等效焦距f对图像的x轴和y轴产生同样的作用，即对畸变点xd和yd的影响相同，所以焦距的大小不影响向量OuPd的方向；世界坐标系原点的定义位置与向量OuPd的方向无关；
步骤S1-3、标定摄像机参数：全部的摄像机参数，包括摄像机有效焦距f、平移矩阵T、旋转矩阵以及一阶畸变系数k，该摄像机参数标定的方法先用线性方法进行线性部分的求解，计算出摄像机旋转矩阵R参数的初始值，然后对包括摄像机有效焦距f、平移矩阵T、以及一阶畸变系数k进行非线性优化；
步骤1-3-1、线性部分的求解：
根据世界坐标和摄像机坐标之间的转化关系，设Zw＝0，则有：
其中，tx、ty为平移参数，对于将(2)两侧同时除以ty并整理可以得到：
上式中，令a1＝r1/ty，a2＝r2/ty，a3＝r3/ty，a4＝r4/ty，a5＝r5/Ty，则上式有五个未知数，对于每一个参考点通过上式得到一个方程，当有超过五个参考点时，采用最小二乘法得到上面五个未知数的解；定义一个矩阵如果该矩阵满秩矩阵，则有
否则的话，有
其中，而ai，aj是剩余不为零的两个参数，有：r1＝a1ty r2＝a2ty r4＝a4ty r5＝a5ty tx＝a3ty，根据所得外参数r1，r2，r4，r5，tx的值和已知世界坐标和图像坐标的对应点，估算出该点在摄像机坐标系下的坐标(Xc，Yc)；由于旋转矩阵R具有单位正交性质，则其余外参数计算方法如下：/>
则有：
r7＝r2r6-r3r5 r8＝r3r4-r1r6 r9＝r1r5-r2r4；
步骤1-3-2、非线性部分的求解：
基于上述计算的参数，根据一阶径向畸变摄像机模型得到如下计算公式：
先设定畸变系数为零，整理上式得到：
其中，F＝r4Xw+r5Yw+ty，G＝r7Xw+r8Yw，选取数个标定点，则上式形成一个超定线性方程组，利用线性最小二乘法求出f和tz，然后将求出的f、tz以及k1＝0作为初始值进一步进行非线性优化计算出f、tz以及k1的真实值；上述标定方法中主点坐标被认为是图像中心点坐标，图像的纵横比sx＝1；利用同一摄像机不同焦距拍摄同一场景得到主点坐标；计算如下：
其中，(u1,v1)和(u2,v2)分别为焦距1和焦距2下的同一个特征点的坐标，选择多的特征点及可求出主点(cx,cy)的坐标；由于摄像机焦距在图像x方向和y方向同时缩放，所以垂直拍摄一个规则物体分别计算图像水平和垂直方向上的像素直径比就可以计算图像的横纵比sx；
利用张正友模板标定法得到主点坐标和横纵比，通过以上步骤，求出摄像机的全部参数，然后导入5个以上图像位置像素坐标点及其在北京54坐标系下对应时刻的船舶AIS坐标点，就可以获得摄像机投影矩阵的全部参数，即摄像机有效焦距f、平移矩阵T、旋转矩阵R以及一阶畸变系数k，根据该投影矩阵，将摄像机图像坐标下的船舶目标的轨迹位置特征值进行转换。
步骤1-4、前景物体分割：
将当前像素值vt(x)与相应背景模型中的每个像素值进行比较，如果最少有1个像素值与当前像素值的距离小于给定阈值γ，其中，γ＝10，那么当前像素点被判断为背景，反之判断为像素值波动；
步骤1-5、背景模型更新：
如果一个像素点被判断为背景，则该像素点随机地替代自己背景模型中的一个像素值；
步骤1-6、前景去噪：该前景去噪包括腐蚀运算和膨胀运算的步骤：
所述的腐蚀运算的步骤包括：
E(f(x,y))＝{min(f(xi,yi)|(xi,yi)∈f(x,y)∩B} (10)
选择算子覆盖下的二值图像的最小像素值作为中心像素的大小，腐蚀运算在该处的作用是去除小面积噪声，B为结构算子，即参与计算的临近像素集；
所述的膨胀运算的步骤包括：
D(f(x,y))＝{max(f(xi,yi)|(xi,yi)∈f(x,y)∩B} (11)
选择算子覆盖下的二值图像的最大像素值作为中心像素的大小，膨胀运算在该处的作用是填补前景区域内小的空缺，B为结构算子，即参与计算的临近像素集。
所述的步骤S2具体为分离目标的合并和多轨迹关联：
步骤2-1、分离目标的合并：
在无遮挡情况下，每一个视频帧检测出的船舶目标相邻，通过聚类方法实现轨迹的连通，具体如下：连通域是由一个最小外接矩形来标记的，矩形框的表示方法为：
[x y width height] (12)
x,y表示的是矩形框的左上角的坐标，横轴x，纵轴y，width表示矩形框的宽度，即横轴长，height表示矩形框的高度，即纵轴长；
当两个矩形框满足公式(13)、(14)的条件时，默认这其中的两个前景块是属于同一个目标的将它们合并成一个大的矩形框；
|xi-xj|＜＝max(widthi,widthj)+width_threshold (13)
|yi-yj|＜＝max(heighti,heightj)+height_threshold (14)
对于合并后的框[X Y widthi,j heighti,j]，如果满足公式(15)、(16)对大小限制、公式(17)、(18)对景深限制、公式(19)对比例限制、公式(11)对有效前景比例中的某任一种限制，则取消合并：
步骤2-1-1、满足大小限制：
合并后的框在高度上不会超过图像高度的五分之二，在宽度上也不会超过图像宽度的五分之一，同时，面积也不会超过两万像素，默认标准图像帧大小归一化为288*352像素，具体是：
heighti,j＞(frame_height)/2.5 (15)
widthi,j＞(frame_width)/5 (16)
其中，[frame_width frame_height]微视频图像的宽度和高度；
步骤2-1-2、满足景深限制：
合并后的位置如果在图像的偏上部分，即处于远景范围内，此时如果面积过大，在真实场景中可能是水面经过带来的暂时性大面积前景或者是前景检测失效，显然不是船舶目标，具体是：
widthi,j*heighti,j＞20000 (17)
Y＜200&widthi,j*heighti,j＞15000 (18)
步骤2-1-3、满足比例限制：
由于船的形状有固定特点，船舶目标的长宽比在预估范围内，经验值是高宽比不会大于6或者小于1，具体是：
heighti,j/widthi,j＞6 or heighti,j/widthi,j＜1 (19)
步骤2-2、多轨迹关联：
对于合并目标形成的轨迹，如果中间出现遮挡或者丢失检测，那么重新出现的目标被认为是新的目标轨迹，造成轨迹完整性缺失，结合目标的内、外部的纹理特征、颜色特征、外部特征进行多轨迹关联：
步骤2-2-1、计算纹理特征LBP：
(1)找到前景图的各个连通域后，按照最小外接矩形的大小在原图上分割出样本，记为l1,对样本l1进行灰度化和大小变换，统一变换到32*16的大小，记为l2；
(2)用大小为3*3的模板遍历样本l2，对模板中心的像素，比较周围8个邻域像素和中心像素的大小，大于中心像素的取1，否则取0，求取LBP特征值，记为LBP1；
(3)对该特征值的二进制数进行循环移位，找到最小的那个值，为该中心像素的旋转不变LBP值，记为LBP2；
(4)模板遍历完样本l2之后，对所有样本的LBP2进行直方图统计，因为是旋转不变8邻域的LBP特征，所以是36维，在36维上统计直方图，得到H1；
(5)不同连通域m、n之间的灰度直方图距离为：
LBP(m,n)＝|Hm-LBP-Hn-LBP| (20)
步骤2-2-2、计算外部特征：
thre＝(α·height+β·v)·Tdisappear (21)
thre表示对每一个连通区域设置不同的“近邻”阈值，α、β是高度和速度的权重，height是目标的高度，v是目标速度，目标最近10帧的平均路程计算所得：
n表示当前连通区域内的目标点个数，ft＝(xt,yt)序列点的中心位置。Tdisappear表示连通区域丢失更新的帧数；
则不同连通域m、n之间的外部特征距离为：
Hist(m,n)＝|threm-traj-thren-traj| (23)
步骤2-2-3、融合更新连通域：
循环比较连通区域集合Traj内各个连通区域之间的特征，进行融合，更新得到新的船舶目标集合Traj：
(1)对于每一个新生成的连通区域traji，记录进入视场时间t′i，搜索其近邻值threi内的k个前景域[blob1,blob2...blobk]，计算与每个前景域的特征匹配度：
Score＝pLBP*LBP(traji,blobj)+pHist*Hist(traji,blobj) (24)
其中，pLBP、pHist分别为加权因子，这里两者均取0.5；
(2)找到匹配度最高，即Score最小的前景域，分配给当前连通域，并更新该连通域的参数height、v；
(3)如果未找到匹配的前景域，则当前连通域traji已达到最大，记录其离开视场时间t″i；
(4)更新，重复融合更新连通域的步骤(1)至步骤(4)，直到找不到可以连通的区域。
所述的步骤S3具体为：
将北京54坐标系的平面直角笛卡尔坐标系，作为进行位置关联的绝对坐标系，CCTV视频图像的识别基于视频图像的相对坐标，AIS船舶位置报告通过GPS所测的经纬度来确定目标位置，发送的目标位置信息是基于WGS-84地理坐标系的经度和纬度；为了使两者进行空间位置的关联，将两种坐标系统一到基于北京54坐标系的平面直角笛卡尔坐标系下的坐标；具体为：
步骤3-1、轨迹转换：
(1)由W84下的大地坐标系(L84，B84，H84)转换为W84下的平面直角坐标(X84，Y84，Z84)，其转换方法为：
其中，W84坐标系下的第一偏心率：卯酉圈曲率半径椭球体长半轴：a1＝6378137m；椭球体短半轴：b1＝6356752.3142m；
(2)采用七参数法将W84下的平面直角坐标(X84，Y84，Z84)转换到北京54体系下的直角坐标(X54，Y54，Z54)，七参数方法的计算公式如下：
七参数法除了考虑三个坐标轴方向的平移(ΔX，ΔY，ΔZ)，也考虑了旋转(α，β，γ)和缩放比例k，ΔX，ΔY，ΔZ，α，β，γ这些参数的数值通过已知的WGS-84地理坐标系和北京54坐标系下的控制点获得；
(3)由北京54体系下的平面直角坐标(X54，Y54，Z54)进行椭球转换为北京54体系下的大地坐标(L54，B54，H54)：
其中，北京54体系下的第一偏心率：卯酉圈曲率半径/>
(4)由北京54体系下的大地坐标(L54，B54，H54)通过高斯-克吕格投影转换为北京54体系下的平面直角坐标(X54，Y54)，高斯投影的性质是投影后角度不变、长度比与点位有关，与方向无关、离中央子午线越远变形越大，高斯投影必须满足中央子午线投影后是直线且长度不变原则，同时必须满足正形投影条件；基于6度带分带的高斯投影的正算公式如下：
其中，x，y为高斯平面直角坐标系的横纵坐标值；L是从中央经线起算的经度值，B是纬度值，需要注意的是二者是以弧度形式表示的；X是赤道到纬度弧度为B处的子午线的长度；η为地球的第二偏心率，N为纬度弧度为B处的卯酉圈曲率半径；
步骤3-2、非线性插值：
在船舶航行中，AIS的报文发送频率随着船舶运行状态动态变化，对AIS航迹进行插值处理，由于船舶在实际航行中需要随航道走向随时调整航速航向，其船舶航迹近似于曲线，并且可以获取船舶在各个位置的航向，因此采用三次样条插值对船舶航迹进行插值拟合，以满足数据融合的条件；
设AIS报文在一个时间段[a，b]内的航迹进行插值，对于AIS航迹在平面直角坐标系上的投影的横坐标x0≤x1...≤xi≤...≤xn，有yi＝∫(xi)，则可以构造一个三次样条插值函数，使它满足3个条件：
(1)yi＝S(xi)，a≤i≤b；
(2)在每一个区间[xi，xi+1]内，S(x)是三次多项式；
(3)在[a，b]内S(x)具有二阶连续导数；
计算出被插值的函数在节点处的导数值，S′(xi)＝mi，则根据条件yi＝S(xi)，yi+1＝S(xi+1)，S′(xi)＝mi，S′(xi+1)＝mi+1，可以使用Hermmite插值公式，计算出每一个小区间[xi，xi+1]内的三次样条函数S(x)的计算公式：
其中hi＝xi+1-xi，对上式的左右两边求二阶导数可得：
因为S(x)在每一个插值点上同样具有二阶连续导数，则可以认为因此有：
根据上式，引入αi，βi，设
由此，对每一个采样点建立一个方程，则有如下方程：
(1-αi)mi-1+2mi+αimi＝βi (36)
上式中，对应n个采样点，得到包含n+1未知变量m0,m1,...,mn的n+1个线性方程组，因为未知变量个数大于方程个数，因此增加一个附加自然边界条件，即被插值函数在两个端点处二阶导数为零，设可得：
其中，α0＝1，α1＝1，至此，方程联立后求出唯一解m0,m1,...,mn，将解出的xi，yim0,m1,...,mn代入上式就可以得到三次样条插值函数；
步骤3-3、初次关联：
对于视频船舶目标traji，以其进入视场的时间ti′为中心,在初级时间关联阈值内进行AIS的数据初级匹配，满足初级关联判断的条件如下示：
kaj-kri＜PATT j＝1,2...,N i＝1,2,...,M (38)
Raj-Rri＜PADT j＝1,2...,N i＝1,2,...,M (39)
其中，kri表示视频图像目标i出现在视场内的时刻，kaj表示AIS目标j出现在视场内的时刻，Raj表示AIS目标j的位置，Rri表示视频图像目标i的位置，PATT表示初级时间关联阈值，PADT表示距离处级关联阈值，设置为穿越视场在航道最大投影的长度；如果没有满足上式的AIS目标船舶，则视频图像内船舶目标均为意思未开启AIS船舶；当上式满足时，认为AIS目标j和视频图像对象i在同一个时间范畴内满足初级关联条件，作为初级匹配的候选对象；
步骤3-4、二级关联：
对待匹配的目标集合进行二级匹配，通过各个特征值的匹配程度加权值进行排序，对排序优先的船舶纳入匹配集合，得到未匹配的视频船舶目标集合与未匹配的AIS船舶目标集合，其中未匹配的视频船舶目标集合即为未开启AIS的船舶目标的；AIS与视频图像中船舶数据关联的结果可以分为两个等级，形成关联结果评价集合W＝{1,0}，其中，1代表关联，0代表不关联；在直积F×W上定义从F到W的模糊关系矩阵就是因素评判矩阵F：
其中，f11ij(k)、f12ij(k)、f13ij(k)、f14ij(k)分别表示位置集合、长度、方向、速度四个单因素集合g1ij(k)、g2ij(k)、g3ij(k)、g4ij(k)在k时刻的关联隶属度，f11ij(k)、f12ij(k)、f13ij(k)、f14ij(k)分别表示单因素集合g1ij(k)、g2ij(k)、g3ij(k)、g4ij(k)在k时刻的不关联隶属度；
设A＝[a1,a2,a3,a4]为位置、长度、方向、速度四个单因素集合关联隶属度的权重集，视频图像船舶目标与AIS目标关联的综合因素关联隶属度可以由单因素权重模糊集A与F的合成运算结果获得；
合成运算的模型包括加权平均型、主观因素决定型以及混合型等。这里采用简单且易于工程实现的加权平均型，有
Cij(k)＝A×F (41)
求出与第i个视频图像目标满足初级关联条件第j个AIS目标的关联隶属度的值，其中隶属度最大的目标，即为匹配的目标：
Max(Cij(k))＝Max(Ci1(k),Ci2(k),...,CiN(k)),(j＝1,2,...,M) (42)
最大值对应的AIS目标即为关联度最大的船只；
当完成视频图像船舶集合与AIS船舶集合的匹配确定重复船舶之后，得到视频图像集合中未匹配船舶集合Traj∩(CU(Traj∩R))，即为未开启AIS的疑似船舶。
本发明由于采集视频的CCTV的视频服务器是海事监管的现有设备，而AIS船台所有航行船舶都强制安装，因此该技术不需要额外的硬件设备成本，通过该方法直接提高了海事对于AIS违规关闭等情况的监管执法效率，节省了人力物力，间接的也提高了辖区水上交通安全的整体水平，经济性和推广前景良好。

Claims (4)

1.一种利用视频图像识别船舶AIS是否开启的方法，其特征在于，包括以下步骤：
S1、船舶视频图像识别步骤；
S2、船舶AIS轨迹分析步骤；
S3、船舶视频图像与船舶AIS轨迹融合步骤，判定未开启AIS的船舶；
所述的步骤S1具体为：利用ViBe运动检测算法检测视频图像背景，包括以下步骤，
步骤S1-1、背景建模：对每个位置的像素点N上有代表性的颜色值vi进行合成，即背景模型Mo(x，y)＝{v1，v2，...，vN}，设定N＝5；
步骤S1-2、摄像机模型建立：采用带有一阶径向畸变的针孔摄像机模型：
(1)对于世界坐标系Ow-XwYwZw及其对应的摄像机坐标系Oc-XcYcZc与图像坐标系Ou-xuyu，世界坐标系下空间中任意一点P(Xw，Yw，Zw)，Pc(xc，yc)是其相应的摄像机坐标点，Pu(xu，yu)是其相应的理想图像坐标点，即针孔成像坐标，Pd(xd，yd)为畸变点图像坐标，只考虑摄像机的一阶径向畸变，无论畸变程度多少，从图像中心点Oc到图像畸变点Pd的向量OuPd的方向保持不变，而且与世界空间坐标系中心点Poz到P的向量PozP平行；由该条件可得：
其中，摄像机的等效焦距f对图像的x轴和y轴产生同样的作用，即对畸变点xd和yd的影响相同，所以焦距的大小不影响向量OuPd的方向；世界坐标系原点的定义位置与向量OuPd的方向无关；
步骤S1-3、标定摄像机参数：全部的摄像机参数，包括摄像机有效焦距f、平移矩阵T、旋转矩阵以及一阶畸变系数k，该摄像机参数标定的方法先用线性方法进行线性部分的求解，计算出摄像机旋转矩阵R参数的初始值，然后对包括摄像机有效焦距f、平移矩阵T、以及一阶畸变系数k进行非线性优化；
步骤1-3-1、线性部分的求解：
根据世界坐标和摄像机坐标之间的转化关系，设Zw＝0，则有：
其中，tx、ty为平移参数，对于将(2)两侧同时除以ty并整理可以得到：
上式中，令a1＝r1/ty，a2＝r2/ty，a3＝r3/ty，a4＝r4/ty，a5＝r5/Ty，则上式有五个未知数，对于每一个参考点通过上式得到一个方程，当有超过五个参考点时，采用最小二乘法得到上面五个未知数的解；定义一个矩阵如果该矩阵满秩矩阵，则有
否则的话，有
其中，而ai，aj是剩余不为零的两个参数，有：r1＝a1ty r2＝a2ty r4＝a4ty r5＝a5ty tx＝a3ty，根据所得外参数r1，r2，r4，r5，tx的值和已知世界坐标和图像坐标的对应点，估算出该点在摄像机坐标系下的坐标(Xc，Yc)；由于旋转矩阵R具有单位正交性质，则其余外参数计算方法如下：
则有：
r7＝r2r6-r3r5 r8＝r3r4-r1r6 r9＝r1r5-r2r4；
步骤1-3-2、非线性部分的求解：
基于上述计算的参数，根据一阶径向畸变摄像机模型得到如下计算公式：
先设定畸变系数为零，整理上式得到：
其中，F＝r4Xw+r5Yw+ty，G＝r7Xw+r8Yw，选取数个标定点，则上式形成一个超定线性方程组，利用线性最小二乘法求出f和tz，然后将求出的f、tz以及k1＝0作为初始值进一步进行非线性优化计算出f、tz以及k1的真实值；上述标定方法中主点坐标被认为是图像中心点坐标，图像的纵横比sx＝1；利用同一摄像机不同焦距拍摄同一场景得到主点坐标；计算如下：
其中，(u1,v1)和(u2,v2)分别为焦距1和焦距2下的同一个特征点的坐标，选择多的特征点及可求出主点(cx,cy)的坐标；由于摄像机焦距在图像x方向和y方向同时缩放，所以垂直拍摄一个规则物体分别计算图像水平和垂直方向上的像素直径比就可以计算图像的横纵比sx；利用张正友模板标定法得到主点坐标和横纵比，通过以上步骤，求出摄像机的全部参数，然后导入5个以上图像位置像素坐标点及其在北京54坐标系下对应时刻的船舶AIS坐标点，就可以获得摄像机投影矩阵的全部参数，即摄像机有效焦距f、平移矩阵T、旋转矩阵R以及一阶畸变系数k，根据该投影矩阵，将摄像机图像坐标下的船舶目标的轨迹位置特征值进行转换；
步骤1-4、前景物体分割：
将当前像素值vt(x)与相应背景模型中的每个像素值进行比较，如果最少有1个像素值与当前像素值的距离小于给定阈值γ，其中，γ＝10，那么当前像素点被判断为背景，反之判断为像素值波动；
步骤1-5、背景模型更新：
如果一个像素点被判断为背景，则该像素点随机地替代自己背景模型中的一个像素值；
步骤1-6、前景去噪：该前景去噪包括腐蚀运算和膨胀运算的步骤：所述的腐蚀运算的步骤包括：
E(f(x,y))＝{min(f(xi,yi)|(xi,yi)∈f(x,y)∩B} (10)
选择算子覆盖下的二值图像的最小像素值作为中心像素的大小，腐蚀运算在该处的作用是去除小面积噪声，B为结构算子，即参与计算的临近像素集；
所述的膨胀运算的步骤包括：
D(f(x,y))＝{max(f(xi,yi)|(xi,yi)∈f(x,y)∩B} (11)
选择算子覆盖下的二值图像的最大像素值作为中心像素的大小，膨胀运算在该处的作用是填补前景区域内小的空缺，B为结构算子，即参与计算的临近像素集；
所述的步骤S2具体为分离目标的合并和多轨迹关联：
步骤2-1、分离目标的合并：
在无遮挡情况下，每一个视频帧检测出的船舶目标相邻，通过聚类方法实现轨迹的连通，具体如下：连通域是由一个最小外接矩形来标记的，矩形框的表示方法为：
[x y width height] (12)
x,y表示的是矩形框的左上角的坐标，横轴x，纵轴y，width表示矩形框的宽度，即横轴长，height表示矩形框的高度，即纵轴长；
当两个矩形框满足公式(13)、(14)的条件时，默认这其中的两个前景块是属于同一个目标的将它们合并成一个大的矩形框；
|xi-xj|＜＝max(widthi,widthj)+width_threshold (13)
|yi-yj|＜＝max(heighti,heightj)+height_threshold (14)
对于合并后的框[X Y widthi,j heighti,j]，如果满足公式(15)、(16)对大小限制、公式(17)、(18)对景深限制、公式(19)对比例限制、公式(11)对有效前景比例中的某任一种限制，则取消合并：
步骤2-1-1、满足大小限制：
合并后的框在高度上不会超过图像高度的五分之二，在宽度上也不会超过图像宽度的五分之一，同时，面积也不会超过两万像素，默认标准图像帧大小归一化为288*352像素，具体是：
heighti,j＞(frame_height)/2.5 (15)
widthi,j＞(frame_width)/5 (16)
其中，[frame_width frame_height]微视频图像的宽度和高度；
步骤2-1-2、满足景深限制：
合并后的位置如果在图像的偏上部分，即处于远景范围内，此时如果面积过大，在真实场景中可能是水面经过带来的暂时性大面积前景或者是前景检测失效，显然不是船舶目标，具体是：
widthi,j*heighti,j＞20000 (17)
Y＜200&widthi,j*heighti,j＞15000 (18)
步骤2-1-3、满足比例限制：
由于船的形状有固定特点，船舶目标的长宽比在预估范围内，经验值是高宽比不会大于6或者小于1，具体是：
heighti,j/widthi,j＞6 or heighti,j/widthi,j＜1 (19)
步骤2-2、多轨迹关联：
对于合并目标形成的轨迹，如果中间出现遮挡或者丢失检测，那么重新出现的目标被认为是新的目标轨迹，造成轨迹完整性缺失，结合目标的内、外部的纹理特征、颜色特征、外部特征进行多轨迹关联：
步骤2-2-1、计算纹理特征LBP：
(1)找到前景图的各个连通域后，按照最小外接矩形的大小在原图上分割出样本，记为l1,对样本l1进行灰度化和大小变换，统一变换到32*16的大小，记为l2；
(2)用大小为3*3的模板遍历样本l2，对模板中心的像素，比较周围8个邻域像素和中心像素的大小，大于中心像素的取1，否则取0，求取LBP特征值，记为LBP1；
(3)对该特征值的二进制数进行循环移位，找到最小的那个值，为该中心像素的旋转不变LBP值，记为LBP2；
(4)模板遍历完样本l2之后，对所有样本的LBP2进行直方图统计，因为是旋转不变8邻域的LBP特征，所以是36维，在36维上统计直方图，得到H1；
(5)不同连通域m、n之间的灰度直方图距离为：
LBP(m,n)＝|Hm-LBP-Hn-LBP| (20)
步骤2-2-2、计算外部特征：
thre＝(α·height+β·v)·Tdisappear (21)
thre表示对每一个连通区域设置不同的“近邻”阈值，α、β是高度和速度的权重，height是目标的高度，v是目标速度，目标最近10帧的平均路程计算所得：
n表示当前连通区域内的目标点个数，ft＝(xt,yt)序列点的中心位置，Tdisappear表示连通区域丢失更新的帧数；
则不同连通域m、n之间的外部特征距离为：
Hist(m,n)＝|threm-traj-thren-traj| (23)
步骤2-2-3、融合更新连通域：
循环比较连通区域集合Traj内各个连通区域之间的特征，进行融合，更新得到新的船舶目标集合Traj：
(1)对于每一个新生成的连通区域traji，记录进入视场时间t′i，搜索其近邻值threi内的k个前景域[blob1,blob2...blobk]，计算与每个前景域的特征匹配度：
Score＝pLBP*LBP(traji,blobj)+pHist*Hist(traji,blobj) (24)
其中，pLBP、pHist分别为加权因子，这里两者均取0.5；
(2)找到匹配度最高，即Score最小的前景域，分配给当前连通域，并更新该连通域的参数height、v；
(3)如果未找到匹配的前景域，则当前连通域traji已达到最大，记录其离开视场时间t″i；
(4)更新，重复融合更新连通域的步骤(1)至步骤(4)，直到找不到可以连通的区域；
所述的步骤S3具体为：
将北京54坐标系的平面直角笛卡尔坐标系，作为进行位置关联的绝对坐标系，CCTV视频图像的识别基于视频图像的相对坐标，AIS船舶位置报告通过GPS所测的经纬度来确定目标位置，发送的目标位置信息是基于WGS-84地理坐标系的经度和纬度；为了使两者进行空间位置的关联，将两种坐标系统一到基于北京54坐标系的平面直角笛卡尔坐标系下的坐标；具体为：
步骤3-1、轨迹转换：
(1)由W84下的大地坐标系(L84，B84，H84)转换为W84下的平面直角坐标(X84，Y84，Z84)，其转换方法为：
其中，W84坐标系下的第一偏心率：卯酉圈曲率半径椭球体长半轴：a1＝6378137m；椭球体短半轴：b1＝6356752.3142m；
(2)采用七参数法将W84下的平面直角坐标(X84，Y84，Z84)转换到北京54体系下的直角坐标(X54，Y54，Z54)，七参数方法的计算公式如下：
七参数法除了考虑三个坐标轴方向的平移(ΔX，ΔY，ΔZ)，也考虑了旋转(α，β，γ)和缩放比例k，ΔX，ΔY，ΔZ，α，β，γ这些参数的数值通过已知的WGS-84地理坐标系和北京54坐标系下的控制点获得；
(3)由北京54体系下的平面直角坐标(X54，Y54，Z54)进行椭球转换为北京54体系下的大地坐标(L54，B54，H54)：
其中，北京54体系下的第一偏心率：卯酉圈曲率半径
(4)由北京54体系下的大地坐标(L54，B54，H54)通过高斯-克吕格投影转换为北京54体系下的平面直角坐标(X54，Y54)，高斯投影的性质是投影后角度不变、长度比与点位有关，与方向无关、离中央子午线越远变形越大，高斯投影必须满足中央子午线投影后是直线且长度不变原则，同时必须满足正形投影条件；基于6度带分带的高斯投影的正算公式如下：
其中，x，y为高斯平面直角坐标系的横纵坐标值；L是从中央经线起算的经度值，B是纬度值，需要注意的是二者是以弧度形式表示的；X是赤道到纬度弧度为B处的子午线的长度；η为地球的第二偏心率，N为纬度弧度为B处的卯酉圈曲率半径；
步骤3-2、非线性插值：
在船舶航行中，AIS的报文发送频率随着船舶运行状态动态变化，对AIS航迹进行插值处理，由于船舶在实际航行中需要随航道走向随时调整航速航向，其船舶航迹近似于曲线，并且可以获取船舶在各个位置的航向，因此采用三次样条插值对船舶航迹进行插值拟合，以满足数据融合的条件；
设AIS报文在一个时间段[a，b]内的航迹进行插值，对于AIS航迹在平面直角坐标系上的投影的横坐标x0≤x1...≤xi≤...≤xn，有yi＝∫(xi)，则可以构造一个三次样条插值函数，使它满足3个条件：
(1)yi＝S(xi)，a≤i≤b；
(2)在每一个区间[xi，xi+1]内，S(x)是三次多项式；
(3)在[a，b]内S(x)具有二阶连续导数；
计算出被插值的函数在节点处的导数值，S′(xi)＝mi，则根据条件yi＝S(xi)，yi+1＝S(xi+1)，S′(xi)＝mi，S′(xi+1)＝mi+1，可以使用Hermmite插值公式，计算出每一个小区间[xi，xi+1]内的三次样条函数S(x)的计算公式：
其中hi＝xi+1-xi，对上式的左右两边求二阶导数可得：
因为S(x)在每一个插值点上同样具有二阶连续导数，则可以认为因此有：
根据上式，引入αi，βi，设
由此，对每一个采样点建立一个方程，则有如下方程：
(1-αi)mi-1+2mi+αimi＝βi (36)
上式中，对应n个采样点，得到包含n+1未知变量m0,m1,...,mn的n+1个线性方程组，因为未知变量个数大于方程个数，因此增加一个附加自然边界条件，即被插值函数在两个端点处二阶导数为零，设可得：
其中，α0＝1，α1＝1，
至此，方程联立后求出唯一解m0,m1,...,mn，将解出的xi，yi m0,m1,...,mn代入上式就可以得到三次样条插值函数；
步骤3-3、初次关联：
对于视频船舶目标traji，以其进入视场的时间t′i为中心,在初级时间关联阈值内进行AIS的数据初级匹配，满足初级关联判断的条件如下示：
kaj-kri＜PATT j＝1,2...,Ni＝1,2,...,M (38)
Raj-Rri＜PADT j＝1,2...,Ni＝1,2,...,M (39)
其中，kri表示视频图像目标i出现在视场内的时刻，kaj表示AIS目标j出现在视场内的时刻，Raj表示AIS目标j的位置，Rri表示视频图像目标i的位置，PATT表示初级时间关联阈值，PADT表示距离处级关联阈值，设置为穿越视场在航道最大投影的长度；如果没有满足上式的AIS目标船舶，则视频图像内船舶目标均为意思未开启AIS船舶；当上式满足时，认为AIS目标j和视频图像对象i在同一个时间范畴内满足初级关联条件，作为初级匹配的候选对象；
步骤3-4、二级关联：
对待匹配的目标集合进行二级匹配，通过各个特征值的匹配程度加权值进行排序，对排序优先的船舶纳入匹配集合，得到未匹配的视频船舶目标集合与未匹配的AIS船舶目标集合，其中未匹配的视频船舶目标集合即为未开启AIS的船舶目标的；AIS与视频图像中船舶数据关联的结果可以分为两个等级，形成关联结果评价集合W＝{1,0}，其中，1代表关联，0代表不关联；在直积F×W上定义从F到W的模糊关系矩阵就是因素评判矩阵F：
其中，f11ij(k)、f12ij(k)、f13ij(k)、f14ij(k)分别表示位置集合、长度、方向、速度四个单因素集合g1ij(k)、g2ij(k)、g3ij(k)、g4ij(k)在k时刻的关联隶属度，f11ij(k)、f12ij(k)、f13ij(k)、f14ij(k)分别表示单因素集合g1ij(k)、g2ij(k)、g3ij(k)、g4ij(k)在k时刻的不关联隶属度；
设A＝[a1,a2,a3,a4]为位置、长度、方向、速度四个单因素集合关联隶属度的权重集，视频图像船舶目标与AIS目标关联的综合因素关联隶属度可以由单因素权重模糊集A与F的合成运算结果获得；
合成运算的模型包括加权平均型、主观因素决定型以及混合型，这里采用简单且易于工程实现的加权平均型，有
Cij(k)＝A×F (41)
求出与第i个视频图像目标满足初级关联条件第j个AIS目标的关联隶属度的值，其中隶属度最大的目标，即为匹配的目标：
Max(Cij(k))＝Max(Ci1(k),Ci2(k),...,CiN(k)),(j＝1,2,...,M) (42)
最大值对应的AIS目标即为关联度最大的船只；
当完成视频图像船舶集合与AIS船舶集合的匹配确定重复船舶之后，得到视频图像集合中未匹配船舶集合Traj∩(CU(Traj∩R))，即为未开启AIS的疑似船舶。
2.根据权利要求1所述的一种利用视频图像识别船舶AIS是否开启的方法，其特征在于，所述的步骤S1为：利用CCTV视频服务器，对船舶进行前景提取和多特征关联，生成已识别待匹配的视频船舶目标集合，每个视频船舶目标进入视场时间、离开视场时间、位置坐标集、船舶长度、船舶速度和航行方向六个特征向量。
3.根据权利要求1所述的一种利用视频图像识别船舶AIS是否开启的方法，其特征在于，所述的步骤S2为：利用船舶AIS报文位置信息，通过对位置信息的坐标转换、航迹插值得到待匹配的AIS船舶目标集合以及每个AIS船舶目标进入视场时间、离开视场时间、位置坐标集、船舶长度、船舶速度和航行方向六个特征向量。
4.根据权利要求1所述的一种利用视频图像识别船舶AIS是否开启的方法，其特征在于，所述的步骤S3为：融合步骤S1和步骤S2获得的待匹配的视频船舶目标集合与待匹配的AIS船舶目标集合，在相同的坐标体系之下，对步骤S1以及步骤S2船舶集合的六个特征向量进行匹配，最后得到匹配船舶集合与未匹配船舶集合，其中，视频船舶目标集合中属于未匹配集合中的就是未开启AIS的船舶。
CN202010539454.1A 2020-06-14 2020-06-14 一种利用视频图像识别船舶ais是否开启的方法 Active CN111931555B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010539454.1A CN111931555B (zh) 2020-06-14 2020-06-14 一种利用视频图像识别船舶ais是否开启的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010539454.1A CN111931555B (zh) 2020-06-14 2020-06-14 一种利用视频图像识别船舶ais是否开启的方法

Publications (2)

Publication Number Publication Date
CN111931555A CN111931555A (zh) 2020-11-13
CN111931555B true CN111931555B (zh) 2023-08-08

Family

ID=73317598

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010539454.1A Active CN111931555B (zh) 2020-06-14 2020-06-14 一种利用视频图像识别船舶ais是否开启的方法

Country Status (1)

Country Link
CN (1) CN111931555B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113990108B (zh) * 2021-10-22 2023-01-20 苏交科集团股份有限公司 一种船舶优化识别和实时跟踪方法及防撞预警系统
CN114972502B (zh) * 2022-04-29 2024-03-15 湖北国际物流机场有限公司 基于拍摄成像的大型船舶海拔高度测算方法
CN115984842A (zh) * 2023-02-13 2023-04-18 广州数说故事信息科技有限公司 一种基于多模态的视频开放标签提取方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102915650A (zh) * 2012-09-21 2013-02-06 交通运输部科学研究院 基于交会摄影的桥梁水域船舶通航安全预警设备
WO2015112263A2 (en) * 2013-12-04 2015-07-30 Urthecast Corp. Systems and methods for processing distributing earth observation images
CN106816038A (zh) * 2017-03-17 2017-06-09 武汉理工大学 一种内河水域异常行为船舶自动识别系统及方法
CN107341828A (zh) * 2017-06-21 2017-11-10 深圳市置辰海信科技有限公司 一种cctv视频目标二维建模与gis投影方法
CN109460740A (zh) * 2018-11-15 2019-03-12 上海埃威航空电子有限公司 基于ais与视频数据融合的船舶身份识别方法
CN109856625A (zh) * 2019-03-06 2019-06-07 国网福建省电力有限公司莆田供电公司 一种基于多源数据融合的船舶位置识别方法
CN110363115A (zh) * 2019-06-28 2019-10-22 上海交通大学 基于ais轨迹数据的船舶作业异常半监督实时检测方法
CN111103977A (zh) * 2019-12-09 2020-05-05 武汉理工大学 一种船舶辅助驾驶数据的处理方法和系统

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102915650A (zh) * 2012-09-21 2013-02-06 交通运输部科学研究院 基于交会摄影的桥梁水域船舶通航安全预警设备
WO2015112263A2 (en) * 2013-12-04 2015-07-30 Urthecast Corp. Systems and methods for processing distributing earth observation images
CN106816038A (zh) * 2017-03-17 2017-06-09 武汉理工大学 一种内河水域异常行为船舶自动识别系统及方法
CN107341828A (zh) * 2017-06-21 2017-11-10 深圳市置辰海信科技有限公司 一种cctv视频目标二维建模与gis投影方法
CN109460740A (zh) * 2018-11-15 2019-03-12 上海埃威航空电子有限公司 基于ais与视频数据融合的船舶身份识别方法
CN109856625A (zh) * 2019-03-06 2019-06-07 国网福建省电力有限公司莆田供电公司 一种基于多源数据融合的船舶位置识别方法
CN110363115A (zh) * 2019-06-28 2019-10-22 上海交通大学 基于ais轨迹数据的船舶作业异常半监督实时检测方法
CN111103977A (zh) * 2019-12-09 2020-05-05 武汉理工大学 一种船舶辅助驾驶数据的处理方法和系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
田璐等.船舶AIS大数据资源管理及分析应用架构设计.《交通运输研究》.2020,(第05期),第31-40页. *

Also Published As

Publication number Publication date
CN111931555A (zh) 2020-11-13

Similar Documents

Publication Publication Date Title
CN111931555B (zh) 一种利用视频图像识别船舶ais是否开启的方法
CN112084869B (zh) 一种基于紧致四边形表示的建筑物目标检测方法
CN109255350B (zh) 一种基于视频监控的新能源车牌检测方法
CN109859226B (zh) 一种图形分割的棋盘格角点亚像素的检测方法
WO2023025236A1 (zh) 一种多通航要素数据融合方法
CN106845364B (zh) 一种快速自动目标检测方法
CN111882586B (zh) 一种面向剧场环境的多演员目标跟踪方法
CN104978567B (zh) 基于场景分类的车辆检测方法
CN110807355A (zh) 一种基于移动机器人的指针仪表检测与读数识别方法
CN111814686A (zh) 一种基于视觉的输电线路识别及异物入侵在线检测方法
CN108121991A (zh) 一种基于边缘候选区域提取的深度学习舰船目标检测方法
CN111681259B (zh) 基于无Anchor机制检测网络的车辆跟踪模型建立方法
CN106778633B (zh) 一种基于区域分割的行人识别方法
CN111738206A (zh) 基于CenterNet的用于无人机巡检的挖掘机检测方法
CN117036641A (zh) 一种基于双目视觉的公路场景三维重建与缺陷检测方法
CN103632376A (zh) 一种两级框架的车辆部分遮挡消除方法
CN110348307B (zh) 一种起重机金属结构攀爬机器人的路径边缘识别方法及系统
CN107169412B (zh) 基于混合模型决策的遥感图像靠港船只检测方法
Yong et al. Real-time traffic cone detection for autonomous vehicle
CN115240086A (zh) 基于无人机的河道船只检测方法、装置、设备及存储介质
CN105740814B (zh) 一种使用视频分析确定固废危废存放状态的方法
CN117649448A (zh) 一种隧道工作面渗漏水智能识别和分割方法
CN113092807A (zh) 基于多目标跟踪算法的城市高架道路车辆测速方法
CN109697418B (zh) 用于场景复原的针对遥感影像路网提取图像的后处理方法
CN115205564B (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