[发明专利]一种DOM栅格数据脱密与恢复方法

摘要

本发明公开了一种DOM栅格数据脱密与恢复方法，属于地理信息安全领域。该方法包括以下步骤：（1）密钥生成过程：确定数据范围，确定数据变换量，计算线性和非线性变换参数；（2）脱密过程：读取密钥文件，根据密钥处理每个像元，保存脱密后的栅格数据；（3）恢复过程：打开脱密后的栅格数据循环处理每个像元。本方法针对DOM栅格数据的安全保护问题，在保证数据拓扑关系不发生改变的前提下，根据密钥可对数据进行脱密，脱密后的数据依据密钥可进行恢复。本方法具有随机性、渐变性、可逆性等特点，提高了DOM栅格数据脱密的可靠性，完善了地理信息安全保护的理论与方法体系，可用于DOM栅格数据的公开发布等方面。

主权项

1.一种DOM栅格数据脱密与恢复方法，其特征在于，包括如下过程：(一)密钥生成过程步骤11，确定数据范围：获取原始栅格数据Raster的最小外接矩形R，R左下角坐标为(x_min，y_min)，右上角坐标为(x_max，y_max)，根据公式(1)得数据中心点坐标(x_mid，y_mid)、数据长度XL和数据宽度YL；xmid=(xmin+xmax)/2ymid=(ymin+ymax)/2XL=xmax-xminYL=ymax-ymin---(1)]]>步骤12，确定数据变换量：具体步骤如下：输入数据总体变换量total，total＞0，非线性变换量nonlinear，0＜nonlinear＜＝total，根据公式(2)得到线性变换量linear；linear=total2-nonlinear2---(2)]]>步骤13，计算线性变换量linear引起的中误差，确定影响变换效果的参数：焦距f、航高H、偏角倾角ω、旋角κ，具体步骤如下：a)焦距f∈(0，1)，b)根据公式(3)计算航高H，HXL*YL/f---(3)]]>c)根据公式(4)计算线性变化量linear的扰动范围linearExtent，linearExtent=linear---(4)]]>d)生成控制点集合，具体步骤如下：在最小外接矩形R范围内生成m*n个均匀网格，在每个网格中随机选取一个控制点，组成源控制点集合FromPoints＝{(Fx_t，Fy_t)|t＝1，2，...m*n}，m*n＞＝10，根据公式(5)计算每个目标控制点坐标(Tx_t，Ty_t)组成目标控制点集合ToPoints＝{(Tx_t，Ty_t)|t＝1，2，...m*n}，Txt=Fxt+dir1×linear+random1×linearExtentTyt=Fyt+dir2×linear+random2×linearExtent---(5)]]>其中：方向参数dir₁∈{1，-1}、dir₂∈{1，-1}，控制点扰动参数random₁和控制点扰动参数random₂在[-1.0，1.0]范围内随机选取，e)坐标归一化，根据公式(6)对预案控制点集合FromPoints和目标控制点集合ToPoints进行归一化处理得到新坐标的源控制点集合FromPoints’＝{(Fx_t’，Fy_t’)|t＝1，2，...m*n}和目标控制点集合ToPoints’＝{(Tx_t’，Ty_t’)|t＝1，2，...m*n}，Fxt′=(Fxt-xmid)*f/HFyt′=(Fyt-ymid)*f/HTxt′=(Txt-mmid)*f/HTyt′=(Tyt-ymid)*f/H---(6)]]>f)计算偏角倾角ω和旋角κ，根据公式(7)利用最小二乘法对新坐标的源控制点集合FromPoints’中源控制点和目标控制点集合ToPoints’中目标控制点进行拟合解算得到偏角倾角ω、旋角κ，g)计算线性变换中误差accuracy₁，具体步骤如下：在最小外接矩形R范围内均匀选取s₁*s₂个样本点，s₁*s₂＞m*n，组成误差计算源控制点集合BPoints＝{(Bx_t，By_t)|t＝1，2，...s₁*s₂}，根据公式(8)和偏角倾角ω、旋角κ计算目标控制点坐标得到误差计算目标控制点集合APoints＝{(Ax_t，Ay_t)|t＝1，2，...s₁*s₂}，(Bx_t’，By_t’)是对误差计算源控制点(Bx_t，By_t)进行归一化处理得到的坐标，根据公式(9)计算中误差accuracy₁，accuracy1=Σ((Axt-Bxt)2+(Ayt-Byt)2)/s1*s2---(9)]]>h)调节目标控制点集合，具体步骤如下：如果|linear/accuracy₁-1|＞0.01，则根据公式(10)调节每个原目标控制点坐标(Tx_t，Ty_t)，得到新的目标控制点坐标(NTx_t，NTy_t)，用新的目标控制点替代原目标控制点即Tx_t＝NTx_t、Ty_t＝NTy_t，得到目标控制点集合ToPoints＝{(Tx_t，Ty_t)|t＝1，2，...m*n}，NTxt=Fxt+(linear/accuracy1)(Txt-Fxt)NTyt=Fyt+(linear/accuracy1)(Tyt-Fyt)---(10)]]>i)循环步骤e)-h)直至|linear/accuracy₁-1|＜＝0.01，得到最终的偏角倾角ω、旋角κ；步骤14，计算非线性变换量nonlinear引起的中误差，确定参数j₀-j₉，具体步骤如下：a)生成控制点高程值Fz_t，利用公式(11)计算控制点高程值最小值hMin，在(hMin，H)和(-H，-hMin)范围内给每个控制点随机选取一个值作为高程值Fz_t，生成三维源控制点集合FromPoints＝{(Fx_t，Fy_t，Fz_t)|t＝1，2，...m*n}，hMin=H*nonlinear/(xmid-xmin)2+(ymid-ymin)2---(11)]]>b)根据公式(12)对生成的三维源控制点集合FromPoints进行最小二乘解算，得到参数j₀-j₉，Fz_t＝j₀+j₁Fx_t+j₂Fy_t+j₃Fx_t²+j₄Fy_t²+j₅Fx_tFy_t+j₆Fx_tFy_t²+j₇Fx_t²Fy_t+j₈Fy_t³+j₉Fy_t³(12)c)计算非线性变换中误差accuracy₂，具体步骤如下：把步骤13中g)步中生成的BPoints作为误差计算源控制点集合，根据公式(13)和参数j₀-j₉解算每个源控制点(Bx_t，By_t)的Bz_t值，得到三维源控制点集合BPoints＝{(Bx_t，By_t，Bz_t)|t＝1，2，...，s₁*s₂}，Bz_t＝j₀+j₁Bx_t+j₂By_t+j₃Bx_t²+j₄By_t²+j₅Bx_tBy_t+j₆Bx_tBy_t²+j₇Bx_t²By_t+j₈By_t³+j₉By_t³(13)根据公式(14)对三维源控制点集合BPoints进行计算得到目标控制点集合APoints＝{(Ax_t，Ay_t)|t＝1，2，...，s₁*s₂}，Axt=(-f(Bxt-xmid)/(Bzt-H))*H/f+xmidAyt=(-f(Byt-ymid)/(Bzt-H))*H/f+ymid---(14)]]>根据公式(15)计算中误差accuracy₂，accuracy2=Σ((Axt-Bxt)2+(Ayt-Bt)2)s1*s2---(15)]]>d)循环步骤a)-c)，直至|nonlinear/accuracy₂-1|＜＝0.01，得到最终参数j₀-j₉；步骤15，焦距f、航高H、偏角倾角ω、旋角κ、数据中心点坐标(x_mid，y_mid)、参数j₀-j₉组成密钥Key，用非对称加密算法RSA对密钥Key进行加密并存入密钥文件Key.txt；(二)脱密过程步骤21，读取密钥文件Key.txt，解密后提取密钥Key；步骤22，打开原始栅格数据Raster，根据原始栅格数据Raster创建新栅格数据CRaster，新栅格数据CRaster的空间参照、栅格起点O(x_o，y_o)、像元长PX与宽PY、栅格行数row与列数column等属性与原始栅格数据Raster相同，但需添加一个波段，即如果原始栅格数据Raster有band个波段，则新栅格数据CRaster有band+1个波段，栅格像素类型改为PT_DOUBLE；步骤23，循环处理每个像元，具体步骤如下：a)获取像元坐标与像元值，设每个像元中心点为像元坐标，根据公式(16)得每个像元的坐标为p_i,j(x_i，j，y_i，j，z_i，j)，若波段值bs＜band+1，则获取原栅格数据像元p_i,j的像元值赋值给v_i,j，若波段值bs＝band+1，则此像元的像元值v_i，j＝z_i,j，xi,j=xo+i*PX+PX/2yi,j=yo-j*PY-PY/2zi,j=j0+j1xi,j+j2yi,j+j3xi,j2+j4yi,j2+j5xi,jyi,j+j6xi,jyi,j2+j7xi,j2yi,j+j8yi,j3+j9yi,j3---(16)]]>其中：i＝0，1，2，...，column-1；j＝0，1，2，...，row-1；b)转换像元坐标，根据公式(17)和密钥Key，对每个像元坐标p_i,j(x_i，j，y_i，j，z_i，j)进行计算，得到像元坐标p_i,j’(x_i,j’，y_i,j’，z_i,j)，c)坐标归一化，根据密钥Key和公式(18)对每个像元坐标p_i,j’进行归一化处理得到脱密后像元坐标p_i,j”(x_i,j”，y_i,j”，z_i,j)，xi,j′′=xmid+xi,j′*H/fyi,j′′=ymid+yi,j′*H/f---(18)]]>d)判断像元p_i,j”是否在新栅格数据CRaster范围内，具体步骤如下：如果x_o≤x_i，j″≤x_o+PX*column，y_o≤y_i，j″≤y_o+PY*row，则像元p_i,j”在新栅格数据CRaster内，根据公式(19)得到像元p_i,j”在新栅格数据CRaster中的对应像元cp_a，b(x_i，j”，y_i,j”)，给像元cp_a，b赋值即像元值cv_a，b＝v_i，j；其中，表示向下取整，步骤24，填充缺失像元，具体步骤如下：如果波段值bs＝band+1，则判断像元cp_a，b周围的8个像元是否每个波段像元值全为0，如果是则为缺失像元，将像元cp_a，b的每个波段的像元值cv_a，b赋值给缺失像元；步骤25，循环步骤23到24，直至每个波段处理完毕；步骤26，如果原栅格数据存在与之相邻的栅格数据，则打开此相邻的栅格数据，判断相邻的栅格数据脱密后是否有像元落在新栅格数据CRaster中，如果有，循环步骤23到步骤26，直至相邻数据处理完毕，保存脱密后的栅格数据CRaster；(三)恢复过程步骤31，读取密钥文件Key.txt，解密后提取密钥Key；步骤32，打开脱密后的栅格数据CRaster，根据脱密后的栅格数据CRaster创建新的栅格数据IRaster，新栅格数据IRaster的空间参照、栅格起点O(x_o，y_o)、像元长PX与宽PY、栅格行数row与列数column等属性与栅格数据CRaster相同，但需减少一个波段，即如果原始栅格数据CRaster有band+1个波段，则新栅格数据IRaster有band个波段；步骤33，循环处理每个像元，具体步骤如下：a)获取像元坐标与像元高程值，设每个像元中心点为像元坐标，根据公式(16)中x_i，j，y_i，j的计算公式得每个像元的坐标为p_i,j(x_i，j，y_i，j)，获取栅格数据CRaster波段值为band+1的像元的像元值并赋值给z_i,j，得到像元的坐标为p_i,j(x_i，j，y_i，j，z_i，j)，若波段值bs＜band+1，则获取原栅格数据像元p_i,j的像元值赋值给v_i,j，b)坐标归一化，根据公式(20)和密钥Key，对每个像元坐标p_i,j(x_i，j，y_i，j，z_i,j)进行计算，得到像元坐标p_i,j’(x_i,j’，y_i,j’，z_i,j)，xi,j′=(xi,j-xmid)*f/Hyi,j′=(yi,j-ymid)*f/H---(20)]]>c)转换像元坐标，根据公式(21)和密钥Key对每个坐标p_i,j’(x_i,j’，y_i,j’，z_i,j)进行计算，得到恢复后的像元坐标p_i,j”(x_i，j”，y_i,j”，z_i，j)，d)判断像元p_i,j”是否在新栅格数据IRaster范围内，具体步骤如下：如果x_o≤x_i，j″≤x_o+PX*column，y_o≤y_i，j″≤y_o+PY*row，则像元p_i,j”在新栅格数据IRaster内，根据公式(19)得到像元p_i,j”在新栅格数据IRaster中的对应像元ip_a，b(x_i，j”，y_i,j”)，给像元ip_a，b赋值即像元值iv_a，b＝v_i,j；步骤34，填充缺失像元，具体步骤如下：如果波段值bs＝band，则判断像元ip_a，b周围的8个像元是否每个波段像元值全为0，如果是则为缺失像元，将像元ip_a，b的每个波段的像元值iv_a，b赋值给缺失像元；步骤35，循环步骤33到34，直至每个波段处理完毕；步骤36，如果脱密栅格数据存在与之相邻的栅格数据，则打开此相邻的栅格数据，判断相邻的栅格数据恢复后是否有像元落在新栅格数据IRaster中，如果有，循环步骤33到步骤36，直至相邻数据处理完毕，保存恢复后的数据IRaster。