本发明提供了一种基于三角函数的矢栅地理数据可逆脱密方法，包括建立三角函数模型，根据原始矢栅地理数据计算三角函数模型参数并建立脱密模型，采用脱密模型对矢栅地理数据进行脱密处理，根据脱密模型参数建立恢复模型，采用恢复模型对脱密后矢栅地理数据进行恢复处理。本方法在维护矢量数据拓扑关系不变的前提下，充分利用三角函数优越的函数特性，通过调整振幅来控制脱密中误差；通过保持模型的周期性，保证脱密数据均匀变化且拓扑保持良好；通过保持其任意阶可导，使得脱密效果光滑渐进；该方法增加了模型抗攻击能力，保证脱密数据具有较高的安全性，同时具有脱密可逆、高效等特点，可为矢栅地理数据的安全保密、共享使用提供技术支持。
1.一种基于三角函数的矢栅地理数据可逆脱密方法，其特征在于：包括以下步骤：(1)建立三角函数模型
采用由正、余弦函数组成的在实数空间内任意阶连续可导的非线性函数模型，其表达式如下所示：f(x)＝b(a sin(ω1x)+(1-a)cos(ω2x)) (1)其中b为振幅参数，形态参数a用于控制函数的形状，ω1，ω2决定三角函数的周期；
根据式(1)的三角函数模型，以x坐标为自变量对y坐标进行三角函数变换，并对x坐标施加扰动 建立基于三角函数的脱密模型：同时，以y坐标为自变量对x坐标进行三角函数变换，并对y坐标施加扰动 建立基于三角函数的脱密模型：从二者中选择一个作为脱密模型；
(22)确定数据范围：输入具有相同空间参考的矢量地理数据和栅格地理数据，获取数据范围；
(23)解算模型参数：在保证三角函数图像均匀变化，没有突变的条件下，设定形态参数a，根据数据范围或者输入的周期T，按照式(4)计算周期参数ω1和ω2：式中，t用于控制三角函数图像形状；
根据输入的目标脱密精度RMS，由式(5)求取偏移量
将计算得到的参数带入式(2)或(3)中建立三角函数脱密模型；
(24)样本点选择：以矢量地理数据和栅格地理数据的图幅范围作为全局样本，取x方向的宽度为xlength，y方向的宽度为ylength，设定间隔为interval，在全局样本中均匀的选取num＝(xlength/interval)*(ylength/interval)个点作为计算全图中误差的样本点；
(25)中误差的迭代控制：采用步骤(23)建立的脱密模型对样本点ps(ssXi，ssYi)，i＝1，
2，...num做脱密处理，得到脱密后的样本点坐标Ts(tsXi，tsYi)，i＝1，2，....num，并由式(7)计算样本点的中误差RMSdecrypt：将RMSdecrypt与脱密指标d进行比较，如果|RMSdecrypt-d|＞ε则根据式(8)调整b来缩放全图中误差：迭代步骤(23)～(25)直到中误差达到目标脱密精度，执行步骤(26)；
(26)输出三角函数脱密模型用于脱密处理：输出模型的五参数a，b，ω1，ω2，(3)矢栅地理数据脱密处理(31)数据输入：将矢量地理数据的所有要素转化为点集Pi(xi，yi)，作为输入，并创建与原始矢量地理数据具有相同数据格式和坐标系统的新矢量地理数据来存储脱密后的数据；
对于栅格地理数据需获取数据的六参数包括左上角点坐标sP0(x0，y0)、像元的长度pX和宽度pY、栅格行数Row和列数Col以及波段数量Band，并创建与原数据具有相同数据格式、空间参考和六参数的新栅格地理数据；
(32)矢量地理数据脱密处理：将步骤(31)获取的原数据点集Pi(xi，yi)带入步骤(26)输出的脱密模型中计算得到脱密后坐标 然后根据脱密后的点集重新创建要素，输入到脱密数据中；
(33)栅格地理数据脱密处理：循环处理所有波段的每个像元，将新图像的行列号tN(i，j)，i＝1，2，...，Col；j＝1，2，...，Row带入式(9)中，计算其对应的坐标tPij(txij，tyij)：根据步骤(26)输出的脱密模型转换得到原图像像元坐标sPij(sxij，syij)，根据公式(10)计算原图像像元坐标对应的像元行列号sN(a，b)：式中， 表示向上取整；
如果0≤a≤Col&&0≤b≤Row，提取原图像中像元对应的像素值vij，将像素值vij赋给tN(i，j)所对应的像元；如果a＜0||a＞Col||b＜0||b＞Row，将0赋给tN(i，j)所对应的像元：(34)输出脱密后数据：将脱密后的矢量地理数据和栅格地理数据叠加输出；
将步骤(26)输出的脱密模型五参数代入式(11)和(12)中建立三角函数恢复模型：式(11)是x方向恢复，式(12)是y方向恢复，输出恢复模型用于脱密数据的恢复处理，根据脱密的方向选择采用x方向恢复还是y方向恢复；
(51)数据输入：将脱密矢量地理数据的所有要素转化为点集Pi(xi，yi)，作为输入，并创建与脱密矢量地理数据具有相同数据格式和空间参考的新矢量地理数据来存储恢复数据；
获取脱密栅格地理数据的六参数包括左上角点坐标sP0(x0，y0)、像元的长度pX和宽度pY、栅格行数Row和列数Col以及波段数量Band，并创建与原数据具有相同数据格式、空间参考和六参数的新栅格地理数据；
(52)矢量地理数据恢复处理：将步骤(51)获取的矢量地理数据点集Pi(xi，yi)带入步骤(4)输出的恢复模型进行计算，得到恢复后坐标Pi′(xi′，yi′)，根据恢复后的点集重新创建要素，并输入到矢量恢复数据中；
(53)栅格地理数据恢复处理：循环处理栅格地理数据的所有波段中的像元，将恢复数据的行列号tN(i，j)，i＝1，2，...，Col；j＝1，2，...，Row代入式(9)中，计算其对应的坐标tPij(txij，tyij)，根据步骤(4)输出的恢复模型将恢复数据坐标转换为脱密数据坐标sPij(sxij，syij)；根据公式(10)计算原数据像元坐标对应的像元行列号sN(a，b)，如果0≤a≤Col&&0≤b≤Row，提取脱密数据中该像素点对应的像素值vij，将像素值vij赋给tN(i，j)对应的像元；如果a＜0||a＞Col||b＜0||b＞Row，将0赋给tN(i，j)对应的像元；
(54)输出恢复数据：将恢复后的矢栅地理数据叠加输出。
一种基于三角函数的矢栅地理数据可逆脱密方法
技术领域
[0001] 本发明涉及地理信息安全领域，具体涉及一种矢栅地理数据可逆脱密方法。
背景技术
[0002] 我国相关法律法规和政策对基础地理数据共享提出了明确的保密要求，其中《测绘管理工作国家秘密范围的规定》、《公开地图内容表示补充规定》和《基础地理信息公开表示内容的规定》等文件对公开地图和地理信息的限制条件和保密内容做出了规定，特别强调基础地理信息及相关要素的空间位置精度的保密性。
[0003] 目前已有的几何精度脱密方法包括投影转换法、图幅变换法、基于移位圆的误差随机干扰法、线性变换法、非线性变换法等，这些方法尚存在一些弊端：一些模型和方法大多存在难以高精度恢复的问题；脱密后拓扑关系保持不够理想；支持矢栅一体脱密的模型和算法不多；部分模型的安全性不高，处理速度慢。
发明内容
[0004] 发明目的：本发明的目的在于针对现有技术的不足，提供一种基于三角函数的矢栅地理数据可逆脱密方法，产生了脱密效果好、脱密性能高且脱密安全性高的作用。
[0005] 技术方案：本发明提供了一种基于三角函数的矢栅地理数据可逆脱密方法，包括以下步骤：
[0006] (1)建立三角函数模型
[0007] 采用由正、余弦函数组成的在实数空间内任意阶连续可导的非线性函数模型，其表达式如下所示：
[0008] f(x)＝b(a sin(ω1x)+(1-a)cos(ω2x)) (1)
[0009] 其中b为振幅值，a用于控制函数的形状，ω1，ω2决定三角函数的周期；
[0010] 三角函数脱密模型并不局限于式(1)所示模型，可以通过对式(1)进行修改，实现不同正弦函数和余弦函数的组合，自定义脱密模型，使得模型多样化，提高模型安全性。
[0011] (2)建立三角函数脱密模型
[0012] (21)设计脱密模型
[0013] 根据式(1)的三角函数模型，以x坐标为自变量对y坐标进行三角函数变换，并对x坐标施加扰动 建立基于三角函数的脱密模型：
[0014]
[0015] 同时，以y坐标为自变量对x坐标进行三角函数变换，并对y坐标施加扰动 建立基于三角函数的脱密模型：
[0016]
[0017] 从二者中选择一个作为脱密模型；
[0018] (22)确定数据范围：输入具有相同空间参考的矢量地理数据和栅格地理数据，获取数据范围；
[0019] (23)解算模型参数：在保证三角函数图像均匀变化，没有突变的条件下，设定形态参数a，根据数据范围或者输入的周期T，按照式(4)计算周期参数ω1和ω2：
[0020]
[0021] 式中，t用于控制三角函数图像形状；
[0022] 根据输入的目标脱密精度RMS，由式(5)求取偏移量
[0023]
[0024] 式中，μ∈[0，1]；
[0025] 根据式(6)求取振幅参数b的初始值：
[0026]
[0027] 将计算得到的参数带入式(2)或(3)中建立三角函数脱密模型；
[0028] (24)样本点选择：以数据的图幅范围作为全局样本，取x方向的宽度为xlength，y方向的宽度为ylength，设定间隔为interval，在全局样本中均匀的选取num＝(xlength/interval)*(ylength/interval)个点作为计算全图中误差的样本点；
[0029] (25)中误差的迭代控制：采用步骤(23)建立的脱密模型对样本点ps(ssXi，ssYi)，i＝1，2，…num做脱密处理，得到脱密后的样本点坐标Ts(tsXi，tsYi)，i＝1，2，…num，并由式(7)计算样本点的中误差RMSdecrypt：
[0030]
[0031] 将RMSdecrypt与脱密指标d进行比较，如果||RMSdecrypt-d||＞ε则根据式(8)调整b来缩放全图中误差：
[0032]
[0033] 迭代步骤(23)～(25)直到中误差达到目标脱密精度，执行步骤(26)；
[0034] (26)输出三角函数脱密模型用于脱密处理：输出模型的五参数a，b，ω1，ω2，[0035] (3)矢栅地理数据脱密处理
[0036] (31)数据输入：将矢量地理数据的所有要素转化为点集，作为输入，并创建与原始矢量地理数据具有相同数据格式和坐标系统的新矢量地理数据来存储脱密后的数据；对于栅格地理数据需获取数据的六参数包括左上角点坐标sP0(x0，y0)、像元的长度pX和宽度pY、栅格行数Row和列数Col以及波段数量Band，并创建与原数据具有相同数据格式、空间参考和六参数的新栅格地理数据；
[0037] (32)矢量地理数据脱密处理：将步骤(31)获取的原数据点集Pi(xi，yi)带入步骤(26)输出的脱密模型中计算得到脱密后坐标 然后根据脱密后的点集重新创建要素，输入到脱密数据中；
[0038] (33)栅格地理数据脱密处理：循环处理所有波段的每个像元，将新图像的行列号tN(i，j)，i＝1，2，…，Col；j＝1，2，…，Row带入式(9)中，计算其对应的坐标tPij(txij，tyij)：
[0039]
[0040] 根据步骤(26)输出的脱密模型转换得到原图像像元坐标sPij(sxij，syij)，根据公式(10)计算原图像像元坐标对应的像元行列号sN(a，b)：
[0041]
[0042] 式中， 表示向上取整；
[0043] 如果0≤a≤Col&&0≤b≤Row，并提取原图像中该像素点对应的像素值vij，将像素值vij赋给tN(i，j)所对应的像元；如果a＜0||a＞Col||b＜0||b＞Row，将0赋给tN(i，j)所对应的像元：
[0044] (34)输出脱密后数据：将脱密后的矢量地理数据和栅格地理数据叠加输出；
[0045] (4)建立三角函数恢复模型
[0046] 将步骤(26)输出的脱密模型五参数代入式(11)和(12)中建立三角函数恢复模型：
[0047]
[0048]
[0049] 式(11)是x方向恢复，式(12)是y方向恢复，输出恢复模型用于脱密数据的恢复处理，根据脱密的方向选择采用x方向恢复还是y方向恢复；
[0050] (5)矢栅地理数据的恢复处理
[0051] (51)数据输入：将脱密矢量地理数据的所有要素转化为点集，作为输入，并创建与脱密矢量地理数据具有相同数据格式和空间参考的新矢量地理数据来存储恢复数据；获取脱密栅格地理数据的六参数包括左上角点坐标sP0(x0，y0)、像元的长度pX和宽度pY、栅格行数Row和列数Col以及波段数量Band，并创建与原数据具有相同数据格式、空间参考和六参数的新栅格地理数据；
[0052] (52)矢量地理数据恢复处理：将步骤(51)获取的矢量地理数据点集Pi(xi，yi)带入步骤(4)输出的恢复模型进行计算，得到恢复后坐标Pi(x′i，y′i)，根据恢复后的点集重新创建要素，并输入到矢量恢复数据中；
[0053] (53)栅格地理数据恢复处理：循环处理栅格地理数据的所有波段中的像元，将恢复数据的行列号tN(i，j)，i＝1，2，…，Col；j＝1，2，…，Row代入式(9)中，计算其对应的坐标tPij(txij，tyij)，根据步骤(4)输出的恢复模型将恢复数据坐标转换为脱密数据坐标sPij(sxij，syij)；根据公式(10)计算原数据像元坐标对应的像元行列号sN(a，b)，如果0≤a≤Col&&0≤b≤Row，提取脱密数据中该像素点对应的像素值vij，将像素值vij赋给tN(i，j)对应的像元；如果a＜0||a＞Col||b＜0||b＞Row，将0赋给tN(i，j)对应的像元；
[0054] (54)输出恢复数据：将恢复后的矢栅地理数据叠加输出。
[0055] 有益效果：本方法在维护矢量数据拓扑关系不变的前提下，充分利用三角函数优越的函数特性，通过调整振幅来控制脱密中误差；通过保持模型的周期性，保证脱密数据均匀变化且拓扑保持良好；通过保持其任意阶可导，使得脱密效果光滑渐进；该方法增加了模型抗攻击能力，保证脱密数据具有较高的安全性，同时具有脱密可逆、高效等特点，可为矢栅地理数据的安全保密、共享使用提供技术支持。
附图说明
[0056] 图1为基于三角函数的矢栅地理数据脱密流程图；
[0057] 图2为实施例选用的原始矢栅地理数据；
[0058] 图3为实施例中脱密前后矢量地理数据对比图；
[0059] 图4(a)(b)为实施例中脱密前后栅格地理数据对比图；
[0060] 图5(a)(b)为实施例中脱密前后矢栅地理数据叠加效果对比图。
具体实施方式
[0061] 下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。
[0062] 实施例：一种基于三角函数的矢栅地理数据可逆脱密方法，如图1所示，本实施例选择某地同一区域tif格式的栅格地理数据和shp格式的矢量地理数据作为原始矢栅地理数据，其坐标系统为WGS84投影坐标系，比例尺为1∶2000。栅格数据的分辨率为0.2m，数据尺寸为3667*3475，原始矢栅地理数据如图2所示，具体操作如下：
[0063] (一)建立三角函数模型
[0064] 采用一个正弦函数和一个余弦函数相加的方式建立三角函数模型，其表达式如下所示：
[0065] f(x)＝b(a sin(ω1x)+(1-a)cos(ω2x)) (1)
[0066] 其中b为振幅值，a用于控制函数的形状，ω1和ω2决定三角函数的周期。
[0067] (二)建立三角函数脱密模型
[0068] 步骤21：设计脱密模型
[0069] 实施例以x坐标为自变量对y坐标进行三角函数变换，并对x坐标施加扰动 建立基于三角函数的脱密模型如式(2)所示：
[0070]
[0071] 步骤22：确定数据范围。输入原始矢栅地理数据，获取坐标范围为x方向3937025.322～3939023.730m、y方向454740.343～456702.510m。
[0072] 步骤23：解算模型参数。设定形态参a＝0.5，设置周期T＝2000，计算得到周期参数ω1＝0.0628318，ω2＝0.0125664；设置目标脱密精度RMS＝50，求取偏移量 设置振幅参数初值为b＝10.1537。
[0073] 步骤24：样本点选择。以数据的图幅范围作为全局样本，取x方向的宽度为xlength＝1998.41m，y方向的宽度为ylength＝1962.167m，在全局样本中均匀的选取255个点作为计算全图中误差的样本点。
[0074] 步骤25：中误差的迭代控制。采用步骤23建立的脱密模型对样本点ps(ssXi，ssYi)，i＝1，2，…255做脱密处理，得到脱密后的样本点坐标Ts(tsXi，tsYi)，i＝1，2，…255，并由式(7)计算样本点的中误差RM5decrypt＝49.243。
[0075] 将RMSdecrypt与脱密指标d进行比较，由于||RMSdecrypt-d||＞0.01，增大b＝10.3082；迭代步骤23到步骤25直到中误差RMSdecrypt＝49.9901满足条件执行步骤26，此时b＝
21.8106。
[0076] 步骤26：输出三角函数脱密模型用于脱密处理。输出模型的五参数a＝0.5，b＝21.8106，ω1＝0.0628318，ω2＝0.0125664，
[0077] (三)矢栅地理数据脱密处理
[0078] 步骤31：数据输入。对于矢量地理数据将所有要素转化为点集，作为输入，并创建与原始矢量地理数据具有相同数据格式和坐标系统的新的矢量地理数据来存储脱密后的数据；对于栅格地理数据需获取数据的六参数，包括左上角点坐标sP0(11454762.91，3938930.714)、像元的长度pX＝0.5611474796975和宽度pY＝0.5611474796978、栅格行数Row＝3667和列数Col＝3475以及波段数量Band＝3，，并创建与原数据具有相同数据格式、空间参考和六参数的新栅格地理数据。
[0079] 步骤32：矢量地理数据脱密处理。以点p(11456114.75，3938186.47)为例进行说明，将其带入步骤26输出的脱密模型进行计算，得到恢复后坐标p′(11456163.74，3938184.031)，循环处理点集中所有的点，并由点集重新创建要素，并输入到矢量脱密数据中。
[0080] 步骤33：栅格地理数据脱密处理。循环处理所有波段的每个像元，以新图像像素点tN(100，129)为例，将其带入式(9)中，计算得到对应的坐标tP(11454819.02，3938858.326)。根据步骤26输出的脱密模型转换原图像像素坐标sP(11454770.03，
3938872.887)。根据公式(10)计算原图像像元坐标对应的像元行列号sN(13，103)，提取原图像中该像素点对应的像素值v＝6.92507127565e-039，将像素值v赋给tN所对应的像元。
[0081] 步骤34：输出脱密后数据。将脱密后的矢量地理数据和栅格地理数据叠加输出。
[0082] (四)建立三角函数恢复模型
[0083] 将步骤26输出的脱密模型五参数代入式(11)中建立三角函数恢复模型：
[0084]
[0085] (五)矢栅地理数据的恢复处理
[0086] 步骤51：数据输入。将脱密矢量地理数据的所有要素转化为点集，作为输入，并创建与脱密矢量地理数据具有相同数据格式和空间参考的新矢量地理数据来存储恢复数据；获取脱密栅格地理数据的六参数，包括左上角点坐标sP0(11454762.91，3938930.714)、像元的长度pX＝0.5611474796975和宽度pY＝0.5611474796978、栅格行数Row＝3667和列数Col＝3475以及波段数量Band＝3，并创建与原数据具有相同数据格式、空间参考和六参数的新栅格地理数据。
[0087] 步骤52：矢量地理数据恢复处理。以点p(11456163.74，3938184.031)为例进行说明，将其带入步骤41输出的恢复模型进行计算，得到恢复后坐标p′(11456163.74，3938184.035)，循环处理点集中所有的点，并由点集重新创建要素，并输入到矢量恢复数据中。
[0088] 步骤53：栅格地理数据恢复处理。
[0089] 循环处理所有波段的每个像元，以新图像的像素tN(13，103)为例，将其带入式(9)中，计算得到对应的坐标tP(11454770.03，3938872.887)。根据步骤26输出的脱密模型转换原图像像素坐标sP(11454819.02，3938858.326)。根据公式(10)计算原图像像素坐标对应的像素行列号sN(100，129)，提取原图像中该像素点对应的像素值v＝6.92507127565e-039，将像素值v赋给tN所对应的像元。
[0090] 步骤54：输出恢复数据。将恢复后的矢栅地理数据叠加输出。
[0091] 本发明可对矢栅地理数据进行几何精度脱密，脱密后的矢量地理数据可进行无损恢复。
[0092] 由图3可以看出脱密前后矢量地理数据的拓扑关系保持较好，图形变形不明显，说明模型对矢量数据的脱密效果较好。
[0093] 由图4可以看出脱密前后栅格地理数据的像素均匀变化，不存在突变，说明模型对栅格数据的脱密效果较好。
[0094] 由图5可以看到栅格地理数据中白色线状要素表示的道路与一些矢量面状要素的相对位置保持较好，尤其是居民地处栅格道路很清晰，大块湖泊的轮廓叠加较好。由此可见，矢量数据和栅格数据的叠加效果较好，也证明三角函数模型同时适用于矢栅地理数据的脱密。
