本发明涉及土方施工领域，特别是涉及一种土方量估算方法。
背景技术：
1、在国土整治项目中，准确的国土整治土方量估算应大致满足挖填方平衡的原则，即挖方量与填方量应大致相等。但在实际应用中，由于国土整治项目施工区具有面积小、地形复杂的特点，大气噪声、飞行干扰、硬件误差等因素使得获取的lidar dem通常存在不可忽略的相对系统误差，即相对系统误差产生的土方量变化相对于国土整治施工产生的土方变化量不可忽略，导致施工前后lidar dem直接做差的差值法估算施工土方量结果不合理，表现为估算所得的挖方量和填方量存在较大差异。这种差异背离了土方平衡原则(挖方量与填方量应大致相等)，不符合国土整治项目施工的基本规律。所以现有的差值法不能很好适应国土整治项目土方量估算。
2、为解决由于相对系统误差存在而导致差值法估算国土整治项目施工土方量挖填方不平衡的问题，本发明提出了一种融入专有lidar dem配准算法的国土整治项目施工土方量估算方法。专有lidar dem配准算法使用非工作区的相对系统误差通过多项式拟合出工作区的相对系统误差，进而在差值法估算施工土方量前消除相对系统误差对估算结果的影响。
技术实现思路
1、为解决上述技术问题，本发明提供了一种土方量估算方法，本发明利用非工作区的相对系统误差多项式拟合出工作区的相对系统误差，进而在挖填方计算前消除相对系统误差的影响，从而准确估算国土整治项目施工土方量，实现了土方量估算中的挖填方平衡。
2、本发明解决其技术问题所采用的技术方案是：一种土方量估算方法，包括步骤如下：
3、s1采集两期lidar dem；
4、s2获取非工作区的相对系统误差；
5、s3建立误差拟合模型；
6、s4使用重加权最小二乘法求解误差拟合模型
7、s5根据误差拟合模型计算工作区相对系统误差；
8、s6对施工后工作区消除相对系统误差；
9、s7统计每个像素对应挖填方；
10、s8统计工作区的施工土方量。
11、优选地，s1使用机载激光雷达在工作区及周围区域采集施工前后两期lidar数据，然后将lidar数据处理为高精度dem数据。
12、优选地，s2获取非工作区的相对系统误差：
13、e(x,y)＝dem2(x,y)-dem1(x,y)(x,y)∈non-working area
14、其中e(x,y)是非工作区的相对系统误差；
15、dem2(x,y)是非工作区施工后的lidardem数据；
16、dem1(x,y)是非工作区施工前的lidar dem数据；
17、s3利用非工作区的相对系统误差通过多项式拟合工作区的相对系统误差；
18、e(x,y)＝θ0x3+θ1y3+θ2x2y+θ3xy2+θ4x2+θ5y2+θ6xy+θ7x+θ8y+θ9
19、其中e(x,y)为曲面上每个像素的相对系统误差；
20、(x,y)为每个像素对应的坐标；
21、θ为代求解的参数；
22、将工作区域周围的非工作区域相对系统误差数据带入上式子，得到如下的矩阵：
23、e＝aθ
24、
25、θ＝[θ0 θ1 θ2 θ3 θ4 θ5 θ6 θ7 θ8 θ9]t
26、e＝[e0 e1 ... en-2 en-1]t
27、其中θ为待求解参数矩阵；
28、e为非工作区域n个点的相对系统误差；
29、a为非工作区域n个点的坐标矩阵。
30、优选地，s4迭代重加权最小二乘法主要包含两个步骤，第一个步骤为使用加权最小二乘求解模型，第二个步骤为利用计算误差更新权值矩阵w。
31、优选地，s4使用加权最小二乘法，首先需要设置权值矩阵，如下公式：
32、
33、wi表示非工作区域n个点的权值；
34、得到如下式子的加权最小二乘法的目标函数：
35、
36、目标函数求导：
37、
38、将上式子转化为下式：
39、
40、令导数为0，求解精确结果为下式：
41、atwaθ-atwe＝0
42、θ＝(atwa)-1atwe
43、在迭代过程中具体求解公式为如下公式：
44、θk+1＝(atwka)-1atwke。
45、优选地，s4，更新权值，每次更新权值为如下公式：
46、
47、其中，wik为第k次迭代的权值的第i个元素，其初值为全1的对角矩阵；
48、ei为非工作区第i点的相对系统误差，
49、ai,j为非工作区第i点的坐标系数，
50、θjk为第k次迭代结果的第j个元素，
51、ξ＝103为迭代权值上限，避免部分权值过大导致过拟合。
52、优选地，对于s4，重复迭代重加权最小二乘法的两个步骤，直至权值w满足条件或迭代次数达到预设上限。
53、优选地，s5根据误差拟合模型计算工作区相对系统误差：
54、e(x,y)＝θ0x3+θ1y3+θ2x2y+θ3xy2+θ4x2+θ5y2+θ6xy+θ7x+θ8y+θ9
55、其中e(x,y)是工作区曲面上每个像素的相对系统误差；
56、(x,y)是每个像素对应的坐标；
57、θ是已求解的参数。
58、优选地，s6使用迭代重加权最小二乘法求解得到非工作区域相对系统误差曲面；通过下式在挖填方计算前消除工作区的相对系统误差，从而使挖填方结果更平衡:
59、dem2'＝dem2-e
60、dem2是整治后的工作区的dem数据；
61、dem2'是消除相对系统误差后的整治后工作区的dem数据；
62、e是工作区两期dem数据的相对系统误差。
63、优选地，s7使用配准后两期lidar dem数据工作区做差乘对应像素面积，计算每个像素的土方变化量:
64、e(x,y)＝area*[dem2'(x,y)-dem1(x,y)](x,y)∈working area
65、e(x,y)是工作区每个像素对应的土方变化量；
66、area是每个像素对应的实地面积；
67、dem2'(x,y)是消除相对系统误差后的整治后工作区的lidar dem数据；
68、dem1(x,y)是施工前工作区的lidar dem数据；
69、s7统计工作区的施工土方量：
70、
71、
72、fill是工作区施工填方量；
73、cut是工作区施工填方量。
74、本发明的有益效果：
75、减小了施工前后两期lidar dem数据之间的相对系统误差，提高了差值法估计施工土方量的可靠性。因为工作区周围的非工作区施工前后只有相对系统误差而没有施工土方量的变化，因此可以将非工作区施工前后的lidar dem差视为该区域的相对系统误差，本发明利用非工作区的相对系统误差多项式拟合出工作区的相对系统误差，进而在挖填方计算前消除相对系统误差的影响，从而准确估算国土整治项目施工土方量，实现了土方量估算中的挖填方平衡。
1.一种土方量估算方法，其特征在于，包括步骤如下：
2.根据权利要求1所述的一种土方量估算方法，其特征在于：
3.根据权利要求1所述的一种土方量估算方法，其特征在于：
4.根据权利要求1所述的一种土方量估算方法，其特征在于：
5.根据权利要求4所述的一种土方量估算方法，其特征在于：
6.根据权利要求4所述的一种土方量估算方法，其特征在于：
7.根据权利要求4所述的一种土方量估算方法，其特征在于：
8.根据权利要求1所述的一种土方量估算方法，其特征在于：
9.根据权利要求1所述的一种土方量估算方法，其特征在于：
10.根据权利要求1所述的一种土方量估算方法，其特征在于：
