CN103093414A - 一种dom栅格数据脱密与恢复方法 - Google Patents

Abstract

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

Description

一种DOM栅格数据脱密与恢复方法

技术领域

本发明属于地理信息安全领域，具体涉及一种针对DOM栅格数据脱密与恢复方法。

背景技术

栅格数据作为一种基础地理数据，直接关系到国家安全和利益，其安全保护研究十分重要，在应用上是要加以限制的。栅格数据脱密主要包括降低栅格数据的分辨率和转换空间精度。降低栅格数据的分辨率可以使用重采样改变像元大小的方法，转换空间精度主要是指利用一定的数学模型转换像元在空间中的位置。栅格空间坐标转换主要包括多项式变换、仿射变换、样条变换等。一次多项式和仿射变换是线性变换，易于恢复；样条变换和二次、三次多项式变换是非线性变换，效果较显著，但是难以恢复，所以需要研究变换效果好，在没有密钥的情况下又难以恢复的栅格数据脱密与恢复方法。

发明内容

本发明的目的是针对现有脱密方法存在的缺陷，提供一种非线性混合模型对DOM栅格数据进行脱密与恢复处理的方法，具有误差的随机性、算法的可逆性和难以破解等特点。

本发明的目的是通过以下技术方案实现的：一种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\×\mathrm{linear}+{\mathrm{random}}_1\×\mathrm{linearExtent}\\ {\mathrm{Ty}}_t={\mathrm{Fy}}_t+{\mathrm{dir}}_2\×\mathrm{linear}+{\mathrm{random}}_2\×\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}^{\′}=({\mathrm{Fx}}_t-x_{\mathrm{mid}})*f/H\\ {{\mathrm{Fy}}_t}^{\′}=({\mathrm{Fy}}_t-y_{\mathrm{mid}})*f/H\\ {{\mathrm{Tx}}_t}^{\′}=({\mathrm{Tx}}_t-m_{\mathrm{mid}})*f/H\\ {{\mathrm{Ty}}_t}^{\′}=({\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{\mathrm{\Σ}({({\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{\mathrm{\Σ}({({\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}}^{\′\′}=x_{\mathrm{mid}}+{x_{i,j}}^{\′}*H/f\\ {y_{i,j}}^{\′\′}=y_{\mathrm{mid}}+{y_{i,j}}^{\′}*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}}^{\′}=(x_{i,j}-x_{\mathrm{mid}})*f/H\\ {y_{i,j}}^{\′}=(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。

本发明提出了一种非线性混合模型对DOM栅格数据进行脱密与恢复处理。本方法针对DOM栅格数据的安全保护问题，在保证数据拓扑关系不发生改变的前提下，根据密钥可对数据进行脱密，脱密后的数据依据密钥可进行恢复。本方法具有随机性、渐变性、可逆性等特点，提高了DOM栅格数据脱密的可靠性，完善了地理信息安全保护的理论与方法体系；本发明脱密和恢复后的栅格数据范围和原栅格数据图幅范围相同，如果原栅格数据存在与之相邻的栅格数据且其每个像元都有像元值，则脱密后的栅格数据不会存在黑边和接边问题，大大提高了数据的可用性；DOM栅格数据的恢复算法可用于对脱密后的数据进行实距量测等，本发明可用于DOM栅格数据的公开发布等方面。

附图说明

图1是本发明方法脱密过程流程图。

图2是本发明方法恢复过程流程图。

图3是本发明实施例选用的原始栅格数据。

图4是本发明实施例脱密后的栅格数据。

图5是本发明实施例恢复后的栅格数据

具体实施方式

下面结合附图对本发明的实施例做详细说明。

本实施例选择tif格式栅格数据，对数据进行读取、脱密与恢复操作，进一步详细说明本发明。本实施例选择某一栅格数据(如图3)作为原始栅格数据，包括以下步骤：

(一)密钥生成过程

步骤11，确定数据范围：获取原始栅格数据Raster的最小外接矩形R，R左下角坐标为(299000，354000)，右上角坐标为(299500.00486，354500.00486)，根据公式(1)得数据中心点坐标(x_mid，y_mid)、数据长度XL和数据宽度YL；

$\left\{\begin{array}{c}{x}_{\mathrm{mid}}=299250.00243\\ y_{\mathrm{mid}}=354250.00243\\ \mathrm{XL}=500.004859999986\\ \mathrm{Yl}=500.004859999986\end{array}\right.$

步骤12，确定数据变换量：具体步骤如下：输入数据总体变换量total＝50，非线性变换量nonlinear＝10，根据公式(2)得到线性变换量linear＝48.9897948556636；

步骤13，计算线性变换量linear引起的中误差，确定影响变换效果的参数：焦距f、航高H、偏角

倾角ω、旋角κ，具体步骤如下：

a)焦距f＝0.15，

b)根据公式(3)计算航高H＝3333.36573333324，

c)根据公式(4)计算线性变化量linear的扰动范围linearExtent＝6.99927102316117，

d)生成控制点集合，具体步骤如下：在最小外接矩形R范围内生成4*3个均匀网格，在每个网格中随机选取一个控制点，组成源控制点集合FromPoints＝{(Fx_t，Fy_t)|t＝1，2，...4*3}，根据公式(5)计算每个目标控制点坐标(Tx_t，Ty_t)组成目标控制点集合ToPoints＝{(Tx_t，Ty_t)|t＝1，2，...4*3}，

e)坐标归一化，根据公式(6)对源控制点集合FromPoints和目标控制点集合ToPoints进行归一化处理得到新坐标的源控制点集合FromPoints’＝{(Fx_t’，Fy_t’)|t＝1，2，...4*3}和目标控制点集合ToPoints’＝{(Tx_t’，Ty_t’)|t＝1，2，...4*3}，

f)计算偏角倾角ω和旋角κ，根据公式(7)利用最小二乘法对新坐标的源控制点集合FromPoints’中源控制点和目标控制点集合ToPoints’中目标控制点进行拟合解算得到偏角倾角ω＝0.0143859160009854、旋角κ＝0.00171113885383085，

g)计算线性变换中误差accuracy1，具体步骤如下：在最小外接矩形R范围内均匀选取100*100个样本点，组成误差计算源控制点集合BPoints＝{(Bx_t，By_t)|t＝1，2，...100*100}，根据公式(8)和偏角

倾角ω和旋角κ计算目标控制点坐标得到误差计算目标控制点集合APoints＝{(Ax_t，Ay_t)|t＝1，2，...100*100}，(Bx_t’，By_t’)是对误差计算源控制点(Bx_t，By_t)进行归一化处理得到的坐标，根据公式(9)计算中误差accuracy₁＝68.0488754948298，

${\mathrm{accuracy}}_1=\sqrt{\mathrm{\Σ}({({\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，...4*3}，

i)循环步骤e)-h)，当accuracy₁＝49.0043372348528时|linear/accuracy₁-1|＜＝0.01，得到最终的偏角

倾角ω＝0.0103597661835137、旋角κ＝0.00123380279867298；

步骤14，计算非线性变换量nonlinear引起的中误差，确定参数j₀-j₉，具体步骤如下：

a)生成控制点高程值Fz_t，利用公式(11)计算控制点高程值最小值hMin＝7.0533234616974，在(hMin，H)和(-H，-hMin)范围内给每个控制点随机选取一个值作为高程值Fz_t，生成三维源控制点集合FromPoints＝{(Fx_t，Fy_t，Fz_t)|t＝1，2，...4*3}，

b)根据公式(12)对生成的三维源控制点集合FromPoints进行最小二乘解算，得到参数j₀＝334993952、j₁＝-8706.607421875、j₂＝6016.97875976563、j₃＝-0.0133545268326998、j₄＝-0.015282241627574、j₅＝0.0288193542510271、j₆＝-3.75795892182396E-08、j₇＝1.34009283669911E-08、j₈＝3.11356203042124E-08、j₉＝7.85700038363757E-09，

c)计算非线性变换中误差accuracy₂，具体步骤如下：把步骤13中g)步中生成的BPoints作为误差计算源控制点集合，根据公式(13)和参数j₀-j₉解算每个源控制点(Bx_t，By_t)的Bz_t值，得到三维源控制点集合BPoints＝{(Bx_t，By_t，Bz_t)|t＝1，2，...，100*100}，根据公式(14)对三维源控制点集合BPoints进行计算得到目标控制点集合APoints＝{(Ax_t，Ay_t)|t＝1，2，...，100*100}，根据公式(15)计算中误差accuracy₂＝171.403130188463，

d)循环步骤a)-c)，直至|nonlinear/accuracy₂-1|＜＝0.01，得到最终参数j₀＝-61101631、j₁＝-1740.00621032715、j₂＝1563.04531860352、j₃＝0.00178630302252714、j₄＝-0.00287107258918695、j₅＝0.00299284062202787、j₆＝-1.17702832480582E-08、j₇＝6.1607428603061E-09、j₈＝-8.11415323820341E-10、j₉＝4.03725997077942E-09；

步骤15，焦距f、航高H、偏角

倾角ω、旋角κ、数据中心点坐标(x_mid，y_mid)、参数j₀-j₉组成密钥Key，用非对称加密算法RSA对密钥Key进行加密并存入密钥文件Key.txt；

(二)脱密过程

步骤21，读取密钥文件Key.txt，解密后提取密钥Key；

步骤22，打开原始栅格数据Raster，根据原始栅格数据Raster创建新栅格数据CRaster，新栅格数据CRaster的空间参照、栅格起点O(299000，354000.00014)、像元长PX＝0.500004859999986与宽PY＝0.500004859999986、栅格行数row＝1000与列数column＝1000等属性与原始栅格数据Raster相同，但需添加一个波段，即如果原始栅格数据Raster有3个波段，则新栅格数据CRaster有4个波段，栅格像素类型改为PT_DOUBLE；

步骤23，循环处理每个像元，具体步骤如下：

a)获取像元坐标与像元值，设每个像元中心点为像元坐标，根据公式(16)得每个像元的坐标为p_i,j(x_i，j，y_i，j，z_i，j)，若波段值bs＜4，则获取原栅格数据像元p_i,j的像元值赋值给v_i,j，若波段值bs＝4则此像元的像元值v_i,j＝z_i,j，

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)，

d)判断像元p_i,j”是否在新栅格数据范围CRaster范围内，具体步骤如下：如果299000≤x_i，j″≤299500.00486，354000≤y_i，j″≤354500.00486，则像元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＝4，则判断像元cp_a，b周围的8个像元是否每个波段像元值全为0，如果是则为缺失像元，将像元cp_a，b的每个波段的像元值cv_a，b赋值给缺失像元；

步骤25，循环步骤23到24，直至每个波段处理完毕；

步骤26，如果原栅格数据存在与之相邻的栅格数据，则打开此相邻的栅格数据，判断相邻的栅格数据脱密后是否有像元落在新栅格数据CRaster中，如果有，循环步骤23到步骤26，直至相邻数据处理完毕，本实施例中原栅格数据存在与之相邻的8个栅格数据，所以需对这8个栅格数据进行循环处理，全部处理完成后保存脱密后的栅格数据CRaster；

(三)恢复过程

步骤31，读取密钥文件Key.txt，解密后提取密钥Key；

步骤32，打开脱密后的栅格数据CRaster，根据脱密后的栅格数据CRaster创建新的栅格数据IRaster，新栅格数据IRaster的空间参照、栅格起点O(299000，354000.00014)、像元长PX＝0.500004859999986与宽PY＝0.500004859999986、栅格行数row＝1000与列数column＝1000等属性与栅格数据CRaster相同，但需减少一个波段，即如果原始栅格数据CRaster有4个波段，则新栅格数据IRaster有3个波段；

步骤33，循环处理每个像元，具体步骤如下：

a)获取像元坐标与像元高程值，设每个像元中心点为像元坐标，根据公式(16)中x_i，j，y_i，j的计算公式得每个像元的坐标为p_i,j(x_i，j，y_i，j)，获取栅格数据CRaster波段值为4的像元的像元值并赋值给z_i,j，得到像元的坐标为p_i,j(x_i，j，y_i，j，z_i，j)，若波段值bs＜4，则获取原栅格数据像元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)，

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范围内，具体步骤如下：如果299000≤x_i，j″≤299500.00486，354000≤y_i，j″≤354500.00486，则像元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＝3，则判断像元ip_a，b周围的8个像元是否每个波段像元值全为0，如果是则为缺失像元，将像元ip_a，b的每个波段的像元值iv_a，b赋值给缺失像元；

步骤35，循环步骤33到34，直至每个波段处理完毕；

步骤36，如果脱密栅格数据存在与之相邻的栅格数据，则打开此相邻的栅格数据，判断相邻的栅格数据恢复后是否有像元落在新栅格数据IRaster中，如果有，循环步骤33到步骤36，直至相邻数据处理完毕，本实施例中原栅格数据存在与之相邻的8个栅格数据，所以需对这8个栅格数据进行循环处理，全部处理完成后保存恢复后的数据IRaster。

本发明脱密前或恢复后可对数据进行重采样以降低栅格数据分辨率，本发明实施例在保证数据拓扑关系不发生改变的前提下对栅格数据进行脱密和恢复，可根据需求设定参数以达到所需的脱密效果，脱密后的数据根据密钥可进行恢复。

Claims (1)

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\×\mathrm{linear}+{\mathrm{random}}_1\×\mathrm{linearExtent}\\ {\mathrm{Ty}}_t={\mathrm{Fy}}_t+{\mathrm{dir}}_2\×\mathrm{linear}+{\mathrm{random}}_2\×\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}^{\′}=({\mathrm{Fx}}_t-x_{\mathrm{mid}})*f/H\\ {{\mathrm{Fy}}_t}^{\′}=({\mathrm{Fy}}_t-y_{\mathrm{mid}})*f/H\\ {{\mathrm{Tx}}_t}^{\′}=({\mathrm{Tx}}_t-m_{\mathrm{mid}})*f/H\\ {{\mathrm{Ty}}_t}^{\′}=({\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{\mathrm{\Σ}({({\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{\mathrm{\Σ}({({\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}}^{\′\′}=x_{\mathrm{mid}}+{x_{i,j}}^{\′}*H/f\\ {y_{i,j}}^{\′\′}=y_{\mathrm{mid}}+{y_{i,j}}^{\′}*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}}^{\′}=(x_{i,j}-x_{\mathrm{mid}})*f/H\\ {y_{i,j}}^{\′}=(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。

Images