本发明涉及地理信息系统领域，特别是一种基于高精度dem数据的快速海陆判断方法。
背景技术：
海洋运输由于运量大、成本低、能耗小、污染少等优势，在我国乃至全球的综合交通运输体系中都占有非常重要的地位。航运管理信息化和智能化系统的相关研究一直是国际航运界关注的热点问题，而实时获得航船舶交通信息是构建智能航运系统的关键。船舶自动识别系统(ais)是目前获取船舶交通信息的一种重要手段，它主要由岸基设施和船载设备共同组成，是一种广播式的船－船、船－岸信息自动交互工具，是重要的船舶助航导航设备。海事监管部门通过岸基ais基站接收船台发送的ais报文，实现对船舶航行动态的实时跟踪与监控。ais船台发送的报文通常包括：mmsi编号、船舶类型、船舶在航状态、外形尺寸等静态信息和船舶速度、航向角、经纬度坐标等动态信息。由于ais报文包含的信息容量非常大，所以近年来成为分析水上交通行为与风险的重要数据来源。
然而，ais缺乏完整的信息验证机制，在实际应用中经常会由于ais传感器接收gps信号过程中产生的数据错误、ais信号传输过程中产生的数据错误、ais设备自身原因产生的数据错误、ais网络通信阻塞导致ais信号发射延迟等问题，导致其动态信息中含有某些错误信息。例如，船舶在正常航行状态下的经纬度坐标有时会异常出现在陆地部分，就是一个比较常见的错误。这种错误的坐标信息会干扰海事监管，导致对水面交通态势的误判，同时也会对以ais数据为基础的船舶交通行为可视化与分析产生困扰。因此，需要对实时接收的ais数据中的经纬度坐标进行海陆判断。此外，ais报文的更新频率非常高，正常情况下，船舶每隔2秒-6分钟就会广播1次ais信息，大量船舶不停广播，导致海量ais数据产生。仅我国海事部门搜集近海船舶动态数据，单日数据量峰值就高达千万条，3个月信息总量超过亿条。这就迫切需要提供一种准确高效的计算方法，有效解决大数据情况下的快速海陆判断问题。
技术实现要素：
本发明的目的是要克服现有技术中利用矢量地图数据进行海陆判断效率低下的问题，提供一种基于高精度dem数据的快速海陆判断方法，可有效解决大数据情况下的快速海陆判断问题。
为达到上述目的，本发明是按照以下技术方案实施的：
一种基于高精度dem数据的快速海陆判断方法，包括以下步骤：
s1、对全球astergdem数据文件进行插值重构获得重构后的全球dem数据文件：
s11、将全球平均分成4096×2048块，每100×100块合并为一个新的dem文件，块内顶点个数为256×256；
s12、在每个块内对256×256个顶点的高程值进行双线性插值，在用每个块内的四个原始顶点进行双线性插值的该过程中，一旦有两个或两个以上原始顶点为无效像素，则认为待插值点为无效像素点，直接将其高程值记录为-9999，以做标识；
s13、将双线性插值后的每100×100块存储为一个新的dem数据文件，后缀名为demx，全球共计861个这样的文件，每个demx文件分为三部分：文件头、块索引、块数据；其中，文件头存储数据的版本号和文件内的分块数目；每个块索引存储该块数据的字节总数和对应的块数据在文件的偏移地址；每个块数据存储该块的元信息和块内所有顶点的高程值；
s2、对插值重构后的全球dem数据文件进行二值化、压缩、合并：
s21、二值化：根据步骤s1中每个demx文件的块数据中顶点的高程值来判断该块的地形特征，分为四种情况：海洋、陆地、海陆接边、未知，依次用整数标识记为0、1、2、3；对于海陆接边块，对海陆接边块内每个顶点继续进行海陆二值化，如果高程值小于等于0，则认为该顶点位于海洋，标识为0；否则，认为其位于陆地，标识为1；
s22、压缩：对于海陆接边块中的每个顶点，每个顶点的标识占用1个bit，每8个相邻的顶点占用1个字节，然后再对二值化后的海陆接边块数据进行zip压缩；
s23、合并：将经过二值化和压缩后的861个demx文件合并为一个全球dem数据文件；
s3、利用合并后的全球dem数据文件进行实时海陆判断：
s31、将步骤s2中经过二值化、压缩、合并的全球dem数据文件存储到计算机的硬盘，在计算机内存中构建一级索引和二级索引，其中一级索引用于记录从文件中读取的所有块的类型标识，每块占据3个字节，即24个bit位；二级索引用于记录所有海陆接边块的原始索引信息，即块的字节数和块在文件中的偏移地址；根据块的地形类型，一级索引情况如下：
如果块为海洋，一级索引的24个bit位中最后3位记为000，其它位均为0，即3个字节从高位到低位依次为0、0、0；
如果块为陆地，一级索引的24个bit位中最后3位记为010，其它位均为0，即3个字节从高位到低位依次为0、0、2；
如果块为未知，一级索引的24个bit位中最后3位记为100，其它位均为0，即3个字节从高位到低位依次为0、0、4；
如果块为海陆接边，一级索引的24个bit位中最后1位记为1，其它23个bit位用于记录该块在二级索引中的序号，即该块在所有海陆接边中的自然顺序；
s32、在内存中构建完两级索引之后，进行快速海陆判断：
s321、根据输入点的经纬度坐标，计算其所在的地形块，并从内存中获取该块的一级索引信息；
s322、解析一级索引的3个字节，如果从高位到低位依次为0、0、0，则说明该块为海洋，从而知该点位于海洋，本次判断结束；如果从高位到低位依次为0、0、2，则说明该块为陆地，从而知该点位于陆地，本次判断结束；如果从高位到低位依次为0、0、4，则说明该块为未知，从而知该点地形未知，本次判断结束；否则，该块为海陆接边，该点位于海洋或者陆地，继续进行步骤③；
s323、根据海陆接边块的全球编号，到缓存中查找是否有该块解压后的海陆二值化数据，如果有，直接转向步骤④；否则，根据一级索引中前面23个bit位的值，反算其在二级索引中的序号，并根据该序号提取其在内存中的二级索引信息，根据二级索引信息中记录的文件偏移地址，到文件对应位置读取该块数据并解压为海陆二值化数据，然后存储到缓存中；
s324、根据该点的经纬度坐标，计算其在该海陆接边块中的行列号，根据行列号计算该点位于二值化数据的哪个字节，以及位于该字节的哪个bit位，如果该bit位为0，表示该点位于海洋，否则，表示该点位于陆地，本次判断结束。
与现有技术相比，本发明可对海量地理坐标进行快速的海陆判断，可广泛应用于航船舶交通信息坐标异常判断、船舶交通行为可视化与分析等领域。
附图说明
图1为本发明按照全球空间索引方案重构的demx文件格式。
图2为本发明合并后的全球demx文件格式。
图3为本发明合并后的全球demx文件栅格化结果。
图4为本发明基于两级索引的快速海陆判断算法流程图。
具体实施方式
下面结合具体实施例对本发明作进一步描述，在此发明的示意性实施例以及说明用来解释本发明，但并不作为对本发明的限定。
1、对全球astergdem数据文件进行插值重构
(1)全球空间索引方案
astergdem数据覆盖范围为北纬83°到南纬83°之间的所有陆地区域，比以往任何地形图都要广得多，达到了地球陆地表面的99％。其空间分辨率为1″，垂直精度达20米，水平精度达30米，是迄今为止可为用户提供的最完整的、精度最高的全球数字高程模型数据。astergdem在全球范围内按照1°×1°进行分片，每个分片包含两个压缩文件，一个数字高程模型(dem)文件和一个质量评估(qa)文件。所有dem文件均以geotiff格式存储，像素个数为3601×3601，海平面高程值为0，无效像素值为-9999。
按照astergdem数据的分片方式，其全球dem数据文件共计64800个，每个文件大小为24.7m字节，数据总量达1.5t字节。数据总量大但同时每个文件又很小，非常不便于进行频繁的io操作，而且不利于快速的海陆判断，这就需要对这些dem数据文件按照一定的空间索引方案进行插值重构。
本发明设计了如下的空间索引方案：将全球平均分成4096×2048块，每100×100块合并为一个新的dem文件，块内顶点个数为256×256，每个顶点的高程值用一个float型变量存储。按照这样的方案，数据水平分辨率为(360°÷4096)÷(256-1)≈1.2″，与原有水平分辨率1″基本相当；每个文件大小为2.44g字节，覆盖地理范围约为8.79°，可以大大提高频繁io操作的效率。
(2)插值重构
首先，在每个块内对256×256个顶点的高程值进行双线性插值。需要注意的是，原始astergdem数据中可能会因为云遮挡等原因出现无效像素值-9999，用该值参与双线性插值，是没有任何意义的。我们的做法是在用四个原始顶点进行双线性插值过程中，一旦有两个或两个以上顶点为无效像素，则认为待插值点为无效像素点，直接将其高程值记录为-9999，以做标识。
接下来，是要将插值后的每100×100块存储为一个新的dem数据文件，后缀名为demx，全球共计861个这样的文件。每个demx文件分为三部分：文件头、块索引、块数据(如图1所示)。其中，文件头存储数据的版本号和文件内的分块数目；每个块索引存储该块数据的字节总数和对应的块数据在文件的偏移地址；每个块数据存储该块的元信息和块内所有顶点(256×256个)的高程值。
2、对重构后的全球dem数据文件进行二值化、压缩、合并
(1)海洋和陆地块的二值化
重构后的全球dem数据分块进行存储，接下来需要根据块数据中顶点的高程值来判断该块的地形特征，分为四种情况：海洋、陆地、海陆接边、未知。
如果所有顶点的高程值全部小于等于0，则该块为海洋，块标识记为0；
如果所有顶点的高程值全部大于0，则该块为陆地，块标识记为1；
否则，该块为海陆接边，块标识记为2；
如果块中没有任何顶点，即字节数为0，则该块为未知(这些块位于最后一个demx文件的最后若干列，即该文件不满100×100块)，块标识记为3。
如果用n表示某顶点在块内的序号，用k表示其在8个相邻顶点中的序号，那么该点的海陆二值化结果vbit为：
其所在字节vbyte的计算方法为：vbyte＝vbyte|vbit。
如果块为海洋、陆地或未知，情况简单，只需要在原有块索引区将块标识记为一个整数即可，不再需要字节数约为256k字节的块数据，再加上全球海洋面积占地球表面积的71.8％，这种二值化做法就可以大大减少全球dem数据的字节数。
如果块为海陆接边，情况就要复杂一些，仅仅记录一个块标识还不够，还需要对块内每个顶点继续进行海陆二值化，以便明确每个顶点到底位于海洋还是陆地。
(2)海陆接边块的压缩
对于海陆接边块中的每个顶点，如果高程值小于等于0，则认为该顶点位于海洋，标识为0；否则，认为其位于陆地，标识为1。为减少数据量，每个顶点的标识占用1个bit，每8个相邻的顶点占用1个字节。这样，海陆接边块的数据量就从256k字节减少到了8k字节。为进一步减少数据量，还需要对二值化后的海陆接边块数据(8k字节)进行zip压缩，可使用开源库zlib完成该操作，将compress2()函数中的参数level设置为9，表示使用最高的压缩率。
实验结果表明，经过二值化和zip压缩之后，每个demx文件的大小从2.44g锐减到1m字节以下，这就大大减少了数据量，提高了io操作的效率，便于接下来的快速海陆判断。
(3)合并为一个数据文件
经过上述步骤处理之后的每个demx文件，不超过1m字节，为进一步提高io操作效率，可将全球861个demx文件合并为一个全球dem数据文件，全球dem数据文件格式如图2所示，将全球dem数据文件栅格化后的结果如图3所示。图2中的全球dem数据文件格式与图1中的文件格式很相似，主要区别在于块数据区不再需要记录海洋、陆地、未知这样三类块的数据了，只需要记录海陆接边块的zip压缩数据即可。注意：在头文件中需要记录四类块(海洋、陆地、海陆接边、未知)的总数，以简化后续的统计工作。
3、利用合并后的全球dem数据文件进行实时海陆判断
(1)在计算机内存中构建一级索引和二级索引
astergdem数据经过重构、二值化、压缩、合并之后，就可以存储到计算机硬盘，准备用于实时的海陆判断了。但要想快速判断出一个点位于海洋还是陆地，不可能频繁地到硬盘中去访问该文件，还需要将文件中的相关信息读取到计算机内存，因为计算机内存访问速度毕竟比io访问速度要快得多。
本发明在计算机内存中构建了两级索引，分别如下：
一级索引：用于记录从文件中读取的所有块的类型标识，每块占据3个字节，即24个bit位；
二级索引：用于记录所有海陆接边块的原始索引信息，即块的字节数和块在文件中的偏移地址。
根据块的地形类型，一级索引情况如下：
如果块为海洋，一级索引的24个bit位中最后3位记为000，其它位均为0，即3个字节从高位到低位依次为0、0、0；
如果块为陆地，一级索引的24个bit位中最后3位记为010，其它位均为0，即3个字节从高位到低位依次为0、0、2；
如果块为未知，一级索引的24个bit位中最后3位记为100，其它位均为0，即3个字节从高位到低位依次为0、0、4；
如果块为海陆接边，一级索引的24个bit位中最后1位记为1，其它23个bit位用于记录该块在二级索引中的序号，即该块在所有海陆接边中的自然顺序。
如果用k表示某个海陆接边块在所有海陆接边块中的序号，那么该块的一级索引值index的3个字节依次为：
v＝(k＜＜1)|0x01
index[0]＝(byte)(v＞＞16)
index[1]＝(byte)((v＜＜8)＞＞16)
index[2]＝(byte)((v＜＜16)＞＞16)。
(2)基于两级索引的快速海陆判断
在内存中构建完两级索引之后，就可以进行快速海陆判断了。如果是海洋块或陆地块，可直接认定待判断的点位于海洋或陆地，算法结束；如果是海陆接边块，还需要解析一级索引的3个字节，计算其在二级索引中的序号。如果分别用index[0]、index[1]、index[2]表示解析后的一级索引的3个字节，那么该块在二级索引中的序号indexsecond为：
indexsecond＝[((int)index[0]＜＜16)|((int)index[1]＜＜8)|(index[2])]＞＞1
根据该二级索引序号，即可到缓存或文件中获取该块的海陆二值化数据，进一步确定该点位于海洋还是陆地。如图4所示，为整个实时判断算法的流程图，主要包括以下几个步骤：
①根据输入点的经纬度坐标，计算其所在的地形块，并从内存中获取该块的一级索引信息；
②解析一级索引的3个字节，如果从高位到低位依次为0、0、0，则说明该块为海洋，从而知该点位于海洋，本次判断结束；如果从高位到低位依次为0、0、2，则说明该块为陆地，从而知该点位于陆地，本次判断结束；如果从高位到低位依次为0、0、4，则说明该块为未知，从而知该点地形未知，本次判断结束；否则，该块为海陆接边，该点位于海洋或者陆地，继续进行步骤③；
③根据海陆接边块的全球编号，到缓存中查找是否有该块解压后的海陆二值化数据，如果有，直接转向步骤④；否则，根据一级索引中前面23个bit位的值，反算其在二级索引中的序号，并根据该序号提取其在内存中的二级索引信息，根据二级索引信息中记录的文件偏移地址，到文件对应位置读取该块数据并解压为海陆二值化数据，然后存储到缓存中；
④根据该点的经纬度坐标，计算其在该海陆接边块中的行列号，根据行列号计算该点位于二值化数据的哪个字节，以及位于该字节的哪个bit位，如果该bit位为0，表示该点位于海洋，否则，表示该点位于陆地，本次判断结束。
本发明的技术方案不限于上述具体实施例的限制，凡是根据本发明的技术方案做出的技术变形，均落入本发明的保护范围之内。