一种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；$\left\{\begin{array}{c}{x}_{\mathrm{mid}}=(x_{\min }+x_{\max })/2\\ y_{\mathrm{mid}}=(y_{\min }+y_{\max })/2\\ \mathrm{XL}=x_{\max }-x_{\min }\\ \mathrm{YL}=y_{\max }-y_{\min }\end{array}\right.---\left(1\right)$步骤12，确定数据变换量：具体步骤如下：输入数据总体变换量total，total＞0，非线性变换量nonlinear，0＜nonlinear＜＝total，根据公式(2)得到线性变换量linear；$\mathrm{linear}=\sqrt{{\mathrm{total}}^2-{\mathrm{nonlinear}}^2}---\left(2\right)$步骤13，计算线性变换量linear引起的中误差，确定影响变换效果的参数：焦距f、航高H、偏角倾角ω、旋角κ，具体步骤如下：a)焦距f∈(0，1)，b)根据公式(3)计算航高H，$H\sqrt{\mathrm{XL}*\mathrm{YL}}/f---\left(3\right)$c)根据公式(4)计算线性变化量linear的扰动范围linearExtent，$\mathrm{linearExtent}=\sqrt{\mathrm{linear}}---\left(4\right)$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}，$\left\{\begin{array}{c}{\mathrm{Tx}}_t={\mathrm{Fx}}_t+{\mathrm{dir}}_1\times \mathrm{linear}+{\mathrm{random}}_1\times \mathrm{linearExtent}\\ {\mathrm{Ty}}_t={\mathrm{Fy}}_t+{\mathrm{dir}}_2\times \mathrm{linear}+{\mathrm{random}}_2\times \mathrm{linearExtent}\end{array}\right.---\left(5\right)$其中：方向参数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}，$\left\{\begin{array}{c}{{\mathrm{Fx}}_t}^{\prime }=({\mathrm{Fx}}_t-x_{\mathrm{mid}})*f/H\\ {{\mathrm{Fy}}_t}^{\prime }=({\mathrm{Fy}}_t-y_{\mathrm{mid}})*f/H\\ {{\mathrm{Tx}}_t}^{\prime }=({\mathrm{Tx}}_t-m_{\mathrm{mid}})*f/H\\ {{\mathrm{Ty}}_t}^{\prime }=({\mathrm{Ty}}_t-y_{\mathrm{mid}})*f/H\end{array}\right.---\left(6\right)$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₁，${\mathrm{accuracy}}_1=\sqrt{\Sigma ({({\mathrm{Ax}}_t-{\mathrm{Bx}}_t)}^2+{({\mathrm{Ay}}_t-{\mathrm{By}}_t)}^2)/s_1*s_2}---\left(9\right)$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}，$\left\{\begin{array}{c}{\mathrm{NTx}}_t={\mathrm{Fx}}_t+(\mathrm{linear}/\mathrm{accurac}{y}_1)({\mathrm{Tx}}_t-{\mathrm{Fx}}_t)\\ {\mathrm{NTy}}_t={\mathrm{Fy}}_t+(\mathrm{linear}/\mathrm{accurac}{y}_1)({\mathrm{Ty}}_t-{\mathrm{Fy}}_t)\end{array}\right.---\left(10\right)$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}，$\mathrm{hMin}=H*\mathrm{nonlinear}/\sqrt{{(x_{\mathrm{mid}}-x_{\min })}^2+{(y_{\mathrm{mid}}-y_{\min })}^2}---\left(11\right)$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₂}，$\left\{\begin{array}{c}{\mathrm{Ax}}_t=(-f({\mathrm{Bx}}_t-x_{\mathrm{mid}})/({\mathrm{Bz}}_t-H))*H/f+x_{\mathrm{mid}}\\ {\mathrm{Ay}}_t=(-f({\mathrm{By}}_t-y_{\mathrm{mid}})/({\mathrm{Bz}}_t-H))*H/f+y_{\mathrm{mid}}\end{array}\right.---\left(14\right)$根据公式(15)计算中误差accuracy₂，${\mathrm{accuracy}}_2=\sqrt{\Sigma ({({\mathrm{Ax}}_t-{\mathrm{Bx}}_t)}^2+{({\mathrm{Ay}}_t-B_t)}^2)s_1*s_2}---\left(15\right)$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，$\left\{\begin{array}{c}{x}_{i,j}=x_o+i*\mathrm{PX}+\mathrm{PX}/2\\ y_{i,j}=y_o-j*\mathrm{PY}-\mathrm{PY}/2\\ z_{i,j}=j_0+j_1{x}_{i,j}+j_2{y}_{i,j}+j_3{x_{i,j}}^2+j_4{y_{i,j}}^2+j_5{x}_{i,j}{y}_{i,j}+j_6{x}_{i,j}{y_{i,j}}^2+j_7{x_{i,j}}^2{y}_{i,j}+j_8{y_{i,j}}^3+j_9{y_{i,j}}^3\end{array}\right.---\left(16\right)$其中：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)，$\left\{\begin{array}{c}{x_{i,j}}^{\prime \prime }=x_{\mathrm{mid}}+{x_{i,j}}^{\prime }*H/f\\ {y_{i,j}}^{\prime \prime }=y_{\mathrm{mid}}+{y_{i,j}}^{\prime }*H/f\end{array}\right.---\left(18\right)$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)，$\left\{\begin{array}{c}{x_{i,j}}^{\prime }=(x_{i,j}-x_{\mathrm{mid}})*f/H\\ {y_{i,j}}^{\prime }=(y_{i,j}-y_{\mathrm{mid}})*f/H\end{array}\right.---\left(20\right)$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。