本发明属于数字海洋与渔业资源的结合技术领域，涉及一种项目用海的三维海底模型的建立方法以及其海洋生物损失量的计算方法，特别涉及填海造地项目对海洋生态环境造成损失的合理计算方法。
背景技术：
不管填海造地工程平面设计如何“科学”，填海本身对海洋生态环境的影响是不可逆的，涉海工程所用海域前期需要经过严格的工程可行性研究、海域使用论证、海洋环境影响评价，施工期需要海洋环境跟踪监测，竣工后要开展环保、海域验收工作，最终验收合格后才可获取海域使用证书。根据我国《海域使用管理法》，我国对海域采取有偿使用制度，用海单位除了缴纳海域使用金，还必须对涉海工程造成的生态损失进行赔偿，如何合理计算海洋生态损失直接关系到用海单位和国家利益，现行的海域使用论证报告书编制过程中，采纳中国水产行业标准《建设项目对海洋生物资源影响评价技术规程》(sc/t9110-2007)进行海洋生物损失和赔偿计算，通常将项目占用海域视作水平投影面来确定占用渔业水域面积或体积，从数字海洋角度出发，这是不科学的，需要更为正确的方法来确定占用渔业水域面积或体积。
技术实现要素：
本发明针对于现行惯用的海洋生物损失量计算方法所存在的上述不足，提供一种通过引入改进的四叉树算法来模拟海底表面，能够科学计算项目占用海底表面积与海域体积，为研究项目用海对海洋生态损失赔偿提供了行之有效的解决方法。
惯用的平面投影面积和海底表面是存在差异的，这些差异对于估算表面覆被生态系统储量、底栖生物资源量和水体中渔业资源量等带来较大的不确定性。
本发明的目的是通过如下技术方案来实现的。
一种基于改进的四叉树算法建立海底地形的德洛内三角网的方法，其特征在于，包括：
获取若干采样点集并创建海底地形的德洛内三角网，其中，所述采样点集为对待模拟海域的海底地形所获得的水深点，所述德洛内三角网是基于改进的四叉树算法建立；
获取不同采样率所对应的海底模型，其中，所述海底模型是根据对所述德洛内三角网中的高层四叉树模型进行可变分辨率插值获得。
在上述技术方案中，获得高层四叉树模型，具体包括：
对不同层级模型内对应的水深点进行采样，在采样时各层级所对应的采样水深点个数为行(2n+1)×列(2n+1)，三角面数为2(2n+1)，其中，n为对应相应层级数，n≥0，计算时保证不溢出即可；
层级间三角面数呈4倍关系，所述不同层级模型的建立方法为：从(2nmax+1)×(2nmax+1)个采样的水深点提取(2nmax-2+1)×(2nmax-2+1)个点建立低分辨率海底地形模型，从(2nmax+1)×(2nmax+1)个水深点提取(2nmax-1+1)×(2nmax-1+1)个点建立中分辨率海底地形模型，从(2nmax+1)×(2nmax+1)个水深点全部提取建立高分辨率海底地形模型，其中nmax为刚好不溢出的最高层级数；
确定所述不同层级模型内对应的水深点所构成的各不同层级采样点集合，保证各层间的视觉连续性，其中，所述采样点集合至少包括各层级模型内对应的水深点的连接点以及后续计算点，所述后续计算点的判定条件为：没有被使用的水深点与采样半径范围内水深点平均值之差绝对值大于设定阈值或者计算没有参与采样的水深点到相应海底地形表面的高差，且高差大于设定阈值，说明上述没有参与计算的水深点对所对应的层级模型地形模拟贡献度较大、特征明显，因此，对该水深点作保留并参加后续计算，直到小于设定阈值；保留高层四叉树模型中的水深最浅点与最深点；
所述不溢出是指水深采样点的个数刚好满足高层四叉树模型建立需求，再往高模拟则会出现漏洞或裂缝。
在上述技术方案中，根据对所述德洛内三角网中的高层四叉树模型进行可变分辨率插值获得海底模型，具体包括：
依据给定采样率，对高层四叉树模型进行重新采样；
对采样后的高层四叉树模型三角顶点进行重新三角剖分，其中，所述三角剖分规则同时满足下述准则a和b：
a.使得经过待剖分的每个三角顶点的外接圆均不包括其余任何点准则；
b.在两相邻三角形形成的凸四边形中，这两三角形中的最小内角一定大于交换凸四边形对角线后所形成的两三角形的最小内角的最大最小角准则。
在上述技术方案中，所述的采样率为40％～100％。
本发明还提供一种计算待模拟海域所对应的海底表面积或海域体积的方法，具体包括：
将组成上述的方法获得的海底模型的每个三角形的面积进行累加，获得待模拟海域所对应的海底表面积；
基于上述的法获得的海底模型，计算每个不规则三角形构成的立柱体体积后对每一立柱体体积进行体积累加，获取所述海域体积。
本申请还提供一种海洋生物损失量计算方法，包括：依据上述的方法计算出海底表面积与海域体积后结合单位面积底栖生物生物量和单位体积鱼卵、仔稚鱼数量，计算项目用海占用海床和水体的生物量分布数据。
本发明的有益效果：
本发明通过改进的四叉树算法建立海底地形的delaunay三角网，使用一定精度较高层四叉树模型进行可变分辨率插值，得到不同采样率的海底模型，并计算表面积和对应海域体积，提出一种项目用海对海洋生态损失赔偿的合理计算方法。本发明方法获得的高层四叉树模型的重新采样率可降低到40％，相对于100％，海底表面积和海域体积计算误差均很小，说明本发明方法能精准地模拟海底表面，建立的高层四叉树模型精确度高，使用较少重采样点即可达到进行准确计算生态损失量的目的。本发明对多分辨率海底地形三维模型实现算法的采纳，是对传统意义上项目用海申报时将项目占用海域基床视作水平投影面来确定占用水域表面积和体积的改进，从数字海洋角度出发，精确评估了涉海工程造成的生态损失额。
附图说明
图1表示高中低层四叉树采样示意图。
图2表示空间delaunay三角形。
图3表示平均海平面与海底表面形成的空间水体。
图4表示单个三角形立柱体。
图5表示基准面关系。
图6表示基于四叉树采样的delaunay算法模拟海底地形生成。
图7表示局部海域不同采样率三角化后面积体积计算差异。
具体实施方法
下述非限制性实施例可以使本领域的普通技术人员更全面地理解本发明，但不以任何方式限制本发明。下述实施例中，如无特殊说明，所使用的实验方法均为常规方法。
实施例1
结合附图1～4，对本发明作进一步说明。
一种基于改进的四叉树算法建立海底地形的delaunay三角网的方法，包括：
获取若干采样点集并创建海底地形的德洛内三角网，其中，所述采样点集为对待模拟海域的海底地形所获得的水深点，所述德洛内三角网是基于改进的四叉树算法建立；
获取不同采样率所对应的海底模型，其中，所述海底模型是根据对所述德洛内三角网图中的高层四叉树模型进行可变分辨率插值获得。
本发明所述的高层四叉树模型的构建，具体包括：
对不同层级模型内对应的水深点进行采样，在采样时各层级所对应的采样水深点个数为行(2n+1)×列(2n+1)，三角面数为2(2n+1)，其中，n为对应相应层级数，n≥0，计算时保证不溢出即可；
层级间三角面数呈4倍关系，所述待测海域层级模型包括低分辨率海底地形模型、中分辨率海底地形模型、以及高分辨率海底地形模型，所述不同层级模型的建立方法为：从(2nmax+1)×(2nmax+1)个采样的水深点提取(2nmax-2+1)×(2nmax-2+1)个点建立低分辨率海底地形模型，从(2nmax+1)×(2nmax+1)个水深点提取(2nmax-1+1)×(2nmax-1+1)个点建立中分辨率海底地形模型，从(2nmax+1)×(2nmax+1)个水深点全部提取建立高分辨率海底地形模型，其中nmax为刚好不溢出的最高层级数；
上述技术方案中，所述不溢出是指水深采样点的个数刚好满足高层四叉树模型建立需求，再往高模拟则会出现漏洞或裂缝。
使用delaunay算法时目标海域内地形水深数据中的每个水深点都会被采样，数据采样算法实现则是从最低分辨率地形模型开始生成，而且较低分辨率地形模型多边形的顶点会被合并到较高的分辨率地形模型中以保证相邻细节层次之间的平滑过渡，并能将模型数据库中相应的组节点组织成四叉树的形式。图1为n＝1,2,3层分辨率水深点采样示意图，其中第3层对应较高分辨率中形成三角面数为128个，9×9个的水深点被采样；第2层对应中分辨率中形成三角面数为32个，5×5个的水深点被采样；第1层对应最低分辨率中形成三角面数为8个，3×3个的水深点被采样。
确定所述不同层级模型内对应的水深点所构成的各组采样点集合，保证各层间的视觉连续性，其中，所述采样点集合至少包括各层级模型内对应的水深点的连接点以及后续计算点，所述后续计算点的判定条件为：没有被使用的水深点与采样半径范围内水深点平均值之差绝对值大于设定阈值或者计算没有参与采样的水深点到相应海底地形表面的高差，且高差大于设定阈值，说明上述没有参与计算的水深点对所对应的层级模型地形模拟贡献度较大、特征明显，因此，对该水深点作保留并参加后续计算，直到小于设定阈值；保留高层四叉树模型中的水深最浅点与最深点。
其中，所述设定阈值可以根据海底地形起伏程度，给出相应的值，一般地形平坦海域，阈值设置较大，地形复杂海域，阈值设置较小。
根据所述采样规则，经(2nmax+1)×(2nmax+1)水深点全部提取，得高层四叉树模型，其中(2nmax+1)×(2nmax+1)刚好小于或等于采样水深点总个数。
进一步，根据对所述高层四叉树模型进行可变分辨率插值获得海底模型，具体包括：
依据给定采样率，对高层四叉树模型进行重新采样；
对采样后的高层四叉树模型三角顶点进行重新三角剖分，其中，所述三角剖分规则同时满足下述准则a和b：
a.使得经过待剖分的每个三角顶点的外接圆均不包括其余任何点准则；
b.在两相邻三角形形成的凸四边形中，这两三角形中的最小内角一定大于交换凸四边形对角线后所形成的两三角形的最小内角的最大最小角准则。
本发明的技术方案中，所述给定采样率为40-100％，基于本发明方法获得的高层四叉树模型的重新采样率可降低到40％，相对于100％，针对实施例2，海底表面积计算误差为0.03％，海域体积计算误差为0.24％，说明本发明方法能精准的模拟海底表面，建立的高层四叉树模型精确度高，使用较少重采样点即可达到进行准确计算生态损失量的目的。
海底地形三维模型生成核心就是确定内插函数，使用有限水深点内插来模拟海底起伏状态。内插函数中要以delaunay(d-tin)算法在地形拟合方面表现最为出色，delaunay算法是一种基于delaunay三角形的地形生成算法，delaunay三角形是一种可以十分有效的对平面上离散点进行三角形剖分的方法，可以很好的反映数字高程模型在平面分布上的距离相关性，该算法遵循两个三角剖分准则，如上述给出的条件a和b。使用delaunay算法确定一个三角面，设任意三角形三个顶点a、b、c，内插点为p，四点坐标分别为(xa，ya，za)、(xb，yb，zb)、(xc，yc，zc)和(xp，yp，zp)，如图2所示：水深zp的计算公式为：
zp＝a0+a1xp+a2yp
根据如上的三角剖分准则，由剖分选定的已知水深点值，设定线性内插函数，计算，最终可获得未知点p点水深值，以此类推进而模拟待测海域其他未采样点水深。
另一方面，本发明中，采用高程中误差指标来评判基于四叉树采样的delaunay算法模拟海底三维地形的精度，所述误差指标判定模拟的delaunay三角网图相对于水深真值的偏离程度，计算公式如式ⅰ：
式中σ为高层四叉树模型的中误差，n为抽样检查点数量即为抽样水深真值的数量，zk为检查点的水深真值，rk为delaunay三角网图的水深模拟值。
由模拟值与验证真值来确定偏离程度，精度评定标准可参考美国地质调查局(usgs)在制定的数字高程模型数据要求中误差不允许超过1/3的等高距，也可参考我国测绘地理信息部门关于不同地形地貌生产数字高程模型数据的有关精度指标规定。
进一步，本发明根据上述得到的海底模型计算待模拟海域所对应的海底表面积与海域体积，具体为：
将组成上述的方法获得的海底模型的每个三角形的面积进行累加，获得待模拟海域所对应的海底表面积；
基于上述的法获得的海底模型，计算每个不规则三角形构成的立柱体体积后对每一立柱体体积进行体积累加，获取所述海域体积。
对于计算所述海底表面积，具体陈述为如下：
生成的海底三维地形(海底模型)由一系列三角形组成，其表面积为某海域内所有三角形面积之和，如图1所示，δabc所三边分别为a、b、c、s表面积计算公式如下：
δabc的面积依据海伦公式进行计算，其中p为三角形周长的一半：
海底三角网表面积为区域内所有三角形面积之和，即，如式ⅱ：
式ⅱ中，n为三角面总数，i为当前三角面，si为当前三角面的表面积。
对于计算所述待测海域体积，具体陈述为如下：
待测海域为所在地理论深度基准面与海底表面所围成的三维空间水体(如图3)，水体的体积计算公式如式ⅲ：
v＝∫∫a[h-ei(x，y)]dσ(式ⅲ)
式ⅲ中：v为水体体积；h为平均海平面高程；ei(x，y)为海底表面高程，考虑到海底表面是个复杂的曲面(用a表示)，理论上将其细分为一些足够小的三角形，dσ为海底表面足够小的面积微元。h-ei(x，y)表示面积微元到平均海平面的高度，用h代替，则公式可写为：
由d-tin模拟海底表面构成三维空间水体，其实际体积计算思路是：先求出每个不规则三角形构成立柱体体积，再把所有三角形立柱体(图4)的体积累加，从而得出整个空间水体体积，每个三角形立柱体体积计算公式，如式ⅵ：
vi＝ai(ha+hb+hc)/3
式ⅵ中：v为水体体积；vi为单个三角形立柱体体积；ai为海底表面三角单元的水平投影面积；a、b、c为海底表面三角单元的三个顶点；ha、hb、hc为海底表面三角单元的三个顶点水深值。
实施例2
以某港区东防波堤工程海域使用对海洋生态损失赔偿的合理计算为例：
1、用海项目与西防波堤工程形成对该港区整个东作业区港池水域的合围，形成未来可供15万吨级单向通行的620m宽的口门，本工程不仅整体起到掩护功能，有利于港区内船舶泊稳作业，而且能够明显起到拦沙导流的作用。
2、用海项目主要建筑内容包括：东防波堤总长度为5827.4m，由四段组成，自港区南围堰工程南段和东段交汇点，沿着规划的三突堤外海侧向外延伸依次分为e-1段、e-2段、e-3段和e-4段，长度分别为1956.4m、510.5m、1815.5m和1545.0m。
3、按照《海籍调查规范》(国海管字[2008]273号)的相关规定执行：“堤坝等非透水构筑物用海，以非透水构筑物(含基床)及其防护设施的水下外缘线为界”。东防波堤的征海范围同时考虑到与已完工的港区南围堰工程相互衔接，最终确定水工结构护坡外缘线水平投影面积为s水平投影面积＝44.4083hm2。
4、本海区属正规半日潮，平均涨潮历时5h，平均落潮历时7～8h，落潮历时大于涨潮历时。工程水深点的深度值基于当地理论深度基准面(图5所示)，据当地水文站1952～1972年和2003年资料统计(以当地理论最低潮面为基准面起算)，工程附近海域主要潮位特征值如下：最高潮位：5.20m(1956年9月4日)，最低潮位：-0.30m(1968年11月10日)，平均高潮位：3.32m，平均低潮位：0.64m，平均潮位(平均海面高度)：2.00m(2003年)，最大潮差：4.46m(2003年6月)，平均潮差：2.74m(2003年)。
依据本发明所述的方法，使用工程附近9321个水深点来生成数字海底地形模型(图6)，东防波堤坡脚线范围内水深点数为1649个，最大水深hmax＝-9.96m，平均水深havg＝-6.41m，计算工程海底水平投影面与平均海面形成海域体积为：
v水平投影-平均海面＝s水平投影面积×(-havg+h海平面)＝44.4083hm2×8.41m＝3734738.03m3
计算工程占用海域基于四叉树采样的delaunay算法三角网表面积(采样率100％)，s表面积＝44.4653hm2，计算工程海底表面积与平均海面形成海域体积为：
5、2011年9月21～23日和2012年5月19日，分别于项目附近海域进行了春秋两次底栖生物调查，采用0.05m2hnm-1型采泥器采集，区域内均匀布设12个站位，每个站位取4个平行样品，取样面积为0.2m2。将采集的样品用0.5mm分样筛淘洗，挑选出所有生物，装入标本瓶内，并标注标签，用5％福尔马林固定液固定。在实验室内用双目立体解剖镜观察、分类、计数，用感量0.01g的电子天平称重(湿重)，具体操作均按《海洋调查规范》(gb12763-2007)执行。将每个站位底栖生物生物量进行换算求出该站位单位面积生物量，取12个站位的平均值，即为该区域底栖生物单位面积生物量值。其中2011年9月航次调查结果：本次调查该海域底栖动物平均生物量为68.571g/m2，其中软体动物最高为44.159g/m2，占平均生物量的64.40％，其次腔肠动物为21.066g/m2，占平均生物量的30.72％；环节动物为1.389g/m2，占平均生物量的2.03％；甲壳动物为1.314g/m2，占平均生物量的1.92％；其它动物为0.643g/m2，占平均生物量的0.94％；2012年5月航次调查结果：本次调查该海域底栖动物平均生物量为31.48g/m2，其中腔肠动物最高为19.33g/m2，占平均生物量的61.41％，其次环节动物为2.69g/m2，占平均生物量的8.53％；软体动物为1.23g/m2，占平均生物量的3.91％；甲壳动物为0.34g/m2，占平均生物量的1.09％；其它动物为7.89g/m2，占平均生物量的25.05％。根据底栖生物春秋两季调查数据，最终将两季调查均值bmean＝50.0255g/m2作为工程附近海区底栖生物生物量，据此估算本工程对底栖生物资源长期损失量。本工程对底栖生物的长期损失主要集中在东防波堤占用海底空间上，在估算损失量时，对底栖生物群落的破坏是不可逆转的，影响期限为20年。将底栖生物平均生物量分别与水平投影面积和表面积进行叠加计算获得东防波堤底栖生物损失量w：
w水平投影面积＝bmean×s水平投影面积×年限＝50.0255(g/m2)×44.4083×104(m2)×20(a)＝444.3095t
w表面积＝bmean×s表面积×年限＝50.0255(g/m2)×44.4653×104(m2)×20(a)＝444.8798t
6、20.11年6月9日～2011年6月14日和2011年9月2日～2011年9月10日，分别于项目附近海域进行了春秋两次鱼卵、仔稚鱼调查，均匀布设12站位，采用浅水i型浮游生物网垂直拖曳和水平拖曳两种方式采集鱼卵。网口流量计计量滤水体积。水平拖网速度大约为1.5nmile/h，采集的样品经5％甲醛海水溶液固定保存后在实验室进行分类鉴定和计数。6月份垂直拖网鱼卵平均密度为0.228粒/m3，仔稚鱼平均密度为5.219尾/m3，优势种为斑尾复鰕虎鱼synechogo-biusommaturus；9月份垂直拖网鱼卵平均密度为0.004粒/m3，仔稚鱼平均密度为0.054尾/m3。优势种为沙氏下鱵鱼hyporhm2mphussajori。根据该海域2011年现状调查资料，本海域春秋两季鱼卵平均密度pmean＝0.116粒/m3，仔稚鱼平均密度emean＝2.6365尾/m3，东防波堤占用海域将不再具有渔业资源生产的功能，按照20a计算鱼卵、仔稚鱼损失量w：
w水平投影-平均海面＝pmean×v水平投影-平均海面×年限＝0.116(粒/m3)×3734738.03(m3)×20(a)
w表面积-平均海面＝pmean×v表面积-平均海面×年限＝0.116(粒/m3)×3608613.62(m3)×20(a)
w水平投影-平均海面＝emean×v水平投影-平均海面×年限＝2.6365(尾/m3)×3734738.03(m3)×20(a)
w表面积-平均面＝emean×v表面积-平均海面×年限＝2.6365(尾/m3)×3608613.62(m3)×20(a)
7、根据《建设项目对海洋生物资源影响评价技术规程》(sc/t9110-2007)的要求，鱼卵和仔稚鱼折算为鱼苗的换算比例，鱼卵生长到商品鱼苗按1％成活率计算，仔稚鱼生长到商品鱼苗按5％成活率计算，底栖生物资源价格按1.0万元/t考虑，鱼苗价格按1元/尾计算。项目用海对海洋生态损失赔偿为：
底栖生物损失额：
m水平投影面积＝m水平投影面积×单价＝444.3095t×1.0(万元/t)＝444.3095万元
m表面积＝m表面积×单价＝444.8798t×1.0(万元/t)＝444.8798万元
鱼卵损失额：
m水平投影-平均海面＝w水平投影-平均海面×1％×单价＝866.46×104粒×1％×1(元/尾)＝8.6646万元
m表面积-平均海面＝w表面积-平均海面×1％×单价＝837.20×104粒×1％×1(元/尾)＝8.3720万元
仔稚鱼损失量：
m水平投影-平均海面＝w水平投影-平均海面×5％×单价＝19693.27×104(尾)×5％×1(元/尾)＝984.6635万元
m表面积-平均海面＝w表面积-平均海面×5％×单价＝19028.22×104(尾)×5％×1(元/尾)＝951.411万元
8、东防波堤用海两种基于两种方法(常规与本发明)统计出的经济补偿金额相差32.9748万元(表1)，常规基于水平投影计算将导致项目建设方多付32.9748万元用于补偿生态损失，数据结果表明：通过改进的四叉树算法建立海底地形的delaunay三角网来计算工程海底表面积与体积并将其应用于资源量估算的计算方法是可行的，在资源量统计比水平投影方案的结果更加准确，这也为资源量统计计算提供了新的科学方案。
表1损失额度比较