【关键词】视频图像；图像放大；FPGA；双线性插值算法
1.引言
随着计算机、电子和数据通信技术的高速发展，实时图像放大技术在夜间辅助驾驶、医学图像分析、工业生产监控等领域得到越来越广泛的应用。大量的图像处理工作需要依靠高实时性且画面清晰的局部动态图像信息来完成，这就对图像放大处理系统的实时性及其放大后的图像品质提出了较高的要求。
图像的放大通常是通过插值完成，在图像放大插值算法中，双线性插值算法放大后的图像比较平滑，具有很好的视觉效果，但由于算法本身计算数据的复杂性，该算法主要通过软件编程实现。然而传统的冯・诺依曼结构的计算机属于标准串行机，编程语言单指令单数据的串行处理过程，将耗用大量的CPU时间，只能应用于一些实时性要求不高的场合。随着可编程逻辑技术的不断进步和创新，现场可编程门阵列FPGA（FieldProgrammableGateArray）具有数以万计的逻辑单元，性能高达500MHz，可通过硬件描述语言，灵活实现各种不同的硬件结构。而其突出的并行处理能力，打破了顺序执行的模式，在每个时钟周期内完成更多的处理任务，其超高强度运算能力，为用硬件实现数字图像处理创造了条件。
本文设计将双线性插值算法应用在以FPGA作为核心处理器的视频图像实时处理系统，使得该系统可在开启放大显示功能后，可以根据外界输入的控制信号，实现对视频图像中不同区域变倍率放大显示的功能。
在这样的背景下，本文综合考虑了图像处理速度和图像质量，针对视频图像的特点，提出了一种基于FPGA的双线性插值算法的硬件实现方案，并通过相关控制模块的设计将其应用于视频图像实时处理系统，完成对视频图像中不同区域实时放大显示的功能。
2.双线性插值算法及其硬件实现框架
2.1双线性插值算法
图像放大算法中，图像的插值过程就是一个图像数据再生的过程。双线性插值是利用新分辨率图像中的像素映射到的原始分辨率图像水平和垂直方向坐标最近的四个点的图像像素，通过加权求和确定，这四个像素的权值与其到需插入像素点的距离成反比。其算法示意图如图1所示。
采用双线性插值算法求新生成坐标H点的数值，需要先对横坐标进行插值，再对纵坐标进行插值，算法计算公式如下：
（1）先对横坐标进行插值
（1）
（2）
（2）再对纵坐标进行插值
（3）
其中，、为点H到点A的权值。
2.2算法硬件实现框架设计
二维线性插值算法的硬件实现可采用量化插值精度的思想。算法中每个固定的插值点所对应的权值为一组固定值，可经过离线计算得到，并将其存储在LUT（查找表）中。对于经过计算不在固定插值位置的，则圆整到固定的插值点中。在实现二维插值算法时，就可根据插值点的位置调出存储在LUT中相应的权值，再通过简单的乘加运算实现。其框架图如图2所示。
基于该思想构造的二维线性插值模块框架在耗用大量的存储资源的同时也难以保证系统的实时性。因而本文考虑将二维线性插值转化为双一维线性插值，即将两个方向上的同时插值转换为分别在水平和垂直两个方向进行插值。按照这个思想，先对整场图像做行内插值，然后再对其做行间插值。其框架图如图3所示。
基于该思想构造的二维线性插值模块框架在大大减少了硬件资源的使用的同时可以较大程度地提高图像插值运算的速度。
考虑到双一维框架结构较为繁杂，本文针对应用系统的便携性，提出一种FPGA实时求取权值的算法硬件实现方法，从而舍弃量化插值和查找对应权值表的过程，构造出更加简洁快速的双一维线性插值模块框架。在权值求取方法设计合理的情况下，该框架省去了权值的片内存储，大大简化了插值控制器的设计的同时也提高了数字实时放大系统的实时性。优化后的双一维线性插值模块框架如图4所示。
3.视频图像实时放大处理系统整体描述
3.1系统结构
本文设计应用的系统具有根据系统外界控制信号实现视频图像中不同区域实时放大显示的功能。其中，系统外界控制信号包括图像放大功能开启信号、待放大区域位置及放大倍数设定信号。系统采用FPGA作为核心处理器，此外还包括视频编、解码芯片、SRAM图像存储器以及针对红外夜视系统的CCD摄像机和LCD显示器。
整个视频放大处理系统的硬件结构图如图5所示。
3.2系统工作原理
系统上电后，FPGA模拟I2C总线对视频编、解码芯片进行初始化配置；CCD摄像头摄入的PAL制模拟电视信号经视频解码芯片SAA7113H后转换为ITU-RBT656格式的数字信号；该信号进入FPGA内部数字视频解码模块，进行有效视频信号中亮度信号的捕捉和行同步信号、场同步信号和奇偶场标志信号的提取；在放大显示功能开启时，视频图像放大处理模块，在同步化控制信号配合下，根据外接输入信号对图像中有效亮度信号进行待放大区域的截取和放大处理；处理后的图像数据送入数字视频编码模块，和同步化控制信号一起合成ITU-RBT.656标准的PAL制式数字视频数据流，该数据流经过两片SRAM乒乓缓存后，送视频编码芯片编码转换为PAL制模拟电视信号输出。
4.视频图像放大处理的实现
4.1视频图像放大处理模块设计
为了使双线性插值放大模块合理的嵌入到系统中根据用户使用需求完成系统功能，视频图像放大处理模块除了设计实现了双线性插值放大模块外，还设计了以下用于控制、显示处理的模块：数据选择模块，用于控制是否对视频图像进行放大处理；放大区域位置设定模块，用于设定需要放大的图像区域在一帧图像中的位置；格式化模块，用于将放大后不满一屏的图像用黑色像素补齐；缓存读写控制模块，用于控制两片SRAM的读写。
4.2双线性插值放大处理模块设计
模块的复杂性使得该模块依然要遵循层次化的设计思想，将其的内部结构进一步划分为时钟管理模块、图像大小计算模块、输入视频缓存控制模块、行内插值模块、行缓存模块、行间插值模块以及双口Ram模块。
时钟管理模块：根据系统输入的放大倍数，对原始数据像素时钟锁相环倍频得到控制时钟。
输入视频缓存控制模块：巧妙地设计行、场同步信号分别作为像素计数器和行计数器复位信号，实现行同步信号HS和像素时钟信号的计数，完成一场图像中待放大区域的截取和存储控制。
双口RAM模块：采取“串并转换”的思想，将同一行原始数据分别储存于两块双口RAM中，利用读端口的地址差，在同一个读时钟，得到计算所需的前后两个像素值，从而以很小的资源牺牲获得了可观的算法实现速度。
行内插值模块：用于从上级缓存模块读入数据完成行内插值，并将计算的结果存入下一级缓存模块中。对于行内插值运算，主要在于其关键参数q和的求取设计：
（4）
（0
其中，j表示一行新的像素的坐标值，表示向下取整运算，n表示放大倍数
设计中考虑将公式中的耗费大量资源的除法运算转换成乘以其倒数的乘法运算，最终由乘加运算完成行内插值过程。行内插值模块实现后的顶层RTL级视图如图6所示。
行缓存功能模块：用于存储两行行内插值模块计算的结果，为后续的行间插值模块提供原始数据，设计的亮点在于行间缓存模块的思路的提出，整个行缓存功能模块顶层原理图如图7所示。
行间插值模块：用于从上级缓存模块读入数据完成行间插值，求取思路同行内插值一致。
5.结论
本文结合视频图像处理系统的特点和FPGA应用技术，对双线性插值算法的硬件实现方案进行了一系列的研究、设计和改进，最终构建出简洁、合理且满足系统实时性要求的算法硬件实现框架，为采用硬件实现理论性的算法提供了思路。
本文完成了视频图像处理系统及视频图像放大处理模块的整体框架设计，将硬件实现后的双线性插值算法模块合理地嵌入到系统中，在各个控制模块的配合下，完成了系统功能，实现了设定区域视频图像的实时放大处理。
本文在设计中充分发挥FPGA并行处理技术优势，实现图像采集与图像的放大处理过程同时进行、并行完成，并将该处理思想贯穿于整个系统设计始终。数字视频解码模块设计避免数据分组处理、减轻存储和算法负担、简化了后续图像处理模块的设计、提高了系统的处理效率；视频图像放大处理模块中双口RAM设计，实现一个读时钟得到计算所需的前后两个原始数据；插值模块关键参数的求取方案，行缓存功能模块中行间缓存模块的设计等都大大提高了系统实时性。同时，设计中视频格式的转换思想的提出，在后续系统升级时，简化了图像去噪、增强等预处理模块的设计。
本文提出的系统整体设计方案综合考虑了图像放大效果和硬件实现代价，在较好地解决放大图像的锯齿问题的同时，也很好地保证了视频图像处理的实时性。经过对该系统的全面仿真和验证，该设计方案能够很好的保证整个系统的实时性和放大后的图像质量，并行处理是保证系统实时性的有效方法。
参考文献
[1]纪红.红外技术基础与应用[M].北京：科学出版社，1979.
[2]王岳环，彭晓.基于FPGA的数字图像实时放大设计[J].计算机工程与应用，2005，8.
[3]赵雅兴.FPGA原理、设计与应用[M].天津：天津大学出版社，1999，4.
[4]俞斯乐.电视原理（第6版）[M].北京：国防工业出版社，2007.
[5]李开宇，张焕春，经亚枝.基于FPGA动态可重构的高速、高质量的图像放大[J].中国图像图形学报，2005，10（1）：69-74.
[6]RafaelC.GonzalezRichardE.Woods，数字图像处理（第二版）[M].北京：电子工业出版社，2003.
[7]PRATTWK，Digitalimageprocessing[M].NewYork：WileyInterscience，1991：207-312.
【关键词】多波束数字测深；正态插值算法；地图重建
NormalInterpolationReconstructionofMultibeamDigitalBathymetryMap
YUJia-cheng1XUXue-qiang2
（1.DepartmentofElectronicsandInformationEngineering，NorthChinaInstituteofScienceandTechnology，HebeiLangfang，065201，China；2.ChinaSouthIndustriesResearchAcademy，Beijing，100089，China）
【Abstract】Themethodthatrealtimemultibeamdigitalbathymetrymapisreconstructedaccordingtothelawthatseabedelevationfollowsnormaldistributionisproposed.Atfirst，inversedistanceweighted（IDW）istransformedtonormaldistributionweighted（NDW），thentherawdepthdataisinterpolatedusingNDWtoimitateactualterrain.Secondly，theimitatedterrainismeasuredalong45degreesnortheastdirectionandrealtimedataarecorrectedbyfivefactors；Thirdly，thecorrecteddataareinterpolatedbyNDWtoconstructrealtimeterrainmap.Theconstructionprocessisverifiedbyusingpracticalwaterdepthdataofcertainseafloor.Simulationresultsshowthatthereconstructedmaperrorisdecreasedby3.3mrelativetoprevious，thusthemapreconstructionprecisionisimprovedobviously.
【Keywords】Multibeamdigitalbathymetry；Normaldistributionweighted（NDW）；Mapreconstruction
多波束水深测量的关键在于建立高分辨率的海底数字地形模型（digitalterrainmodelDTM）[1]，这是目前多波束数据处理研究中的一个重点。为实现大数据量离散深度数据的快速海底数字地形建模和可视化，国内外研究人员通常采用对离散深度数据进行格网化处理[2]，其实质是对深度数据进行内插，使离散分布的数据点变为规则栅格分布。根据数据量的大小，一般分为中小规模数据的插值和大规模数据的插值大类[3]。本文以面换能器多波束系统为例建立仿真模型，对波束的立体宽度，从横纵两个方向对水深测量的影响进行改正，以提高水深测量精度，在此基础上，采用以二维正态分布构建权重的算法来构建规则格网（GRID）的DTM。
1距离反比权重法的改进
1.1算法原理
距离反比权重法（inversedistanceweighted，IDW）基于如下的假定：相互之间距离较近的对象相对于距离较远的对象之间有更强的相似性。在连续的自然地形表面上所采集的离散点高程属性符合这种假定。IDW采用周边已经测定的采样点（已知点）来确定未知点的数值，距离该点位近的已知点对其数值的影响大于距离较远的点。因此，IDW假定数据集中的每一个已知点都对该未知点存在局部的影响力，这种影响力的大小随着距离的增大而逐渐消失。IDW给距离较近的点分配一个比较大的权值，而给较远的点分配一个较小的权值以体现这种影响力[1]。即[4]
设空间待插点为点P（xP，yP，zP），P点邻域内有已知离散点Qi（xi，yi，zi），i=1，2，…，n，利用距离反比权重法对P点的属性值z进行插值。待插点的属性值f（x，y）是待插点邻域内已知离散点属性值的加权平均，
zp=w（d）z/w（d）（1）
式中，w（di）是权函数，其值与待插点与邻域内离散点之间的距离有关，是距离的k次方的倒数，w（di）=1/d2i。k值越低，插值结果越平滑，k值越高，插值结果细节越详细，通常取为2。zi是第i个离散点的数值，di是该点到网格点之间的距离[5]。
1.2算法改M
大量研究表明：自然地形的起伏满足随机分布规律，而DEM是表示地面高程空间分布的数字地面模型，反映自然地形的起伏变化，确实服从一定的统计分布规律，一般可以用正态分布代替[6]，即以未知点为中心，未知点与周围点之间随位置的远近呈现出一种正态分布的关系。未知点对自身的影响最大，具有最大的权值，周围点随着距离的增大权值逐渐减小，符合正态分布。对海底地形而言，点相对未知点所具有的权重为
w（d）=f（x，y）=exp=exp（2）
由于距离越远，影响力越小。随着距离的增大，己知点对于未知点的影响实际上已经可以忽略不计，为了加快运算速度，设定这些点的权重为0。因此，需要指定一个未知点的搜索邻域，只考虑位于该邻域内的已知点对未知点所造成的影响。邻域的形状与数据集的特性有关：如果权值确定不需要考虑方向性的影响，则邻域的形状可以是一个圆形[6]
1.3搜索半径的确定
在文献[7]中，使用一矩形作为搜索区域，对落在矩形窗内的有效数据，都进行插值。为了防止失真，设立有效数据个数门限，当矩形窗内的有效数据个数小于此下限值，则此网格节点为无效值。这个门限缺省值为1，即只有矩形窗内存在数据，该网格节点即为有效点，但是为了确保插值的准确性，推荐值为3。在本文中，为了符合正态权重公式w（d）=exp，采用一圆形区域，如果所用网格为不完全规则的矩形网格，为了兼顾运算速度和确保插值的准确，把半径取为最大网格的对角线长度，与文献[1]不同的是，网格点为已知数据点，格内的任一点为待插点，这样任意待插点都至少有4个已知点为其插值作出贡献，也为网格大小不一时保留一定的余度，以防在待插点周围搜索不到已知点，如图1所示，|AB|=max（经度网格间距），|AC|=max（纬度网格间距），取搜索半径R=|AD|。如果最大经度网格间距与最大纬度网格间距不从属同一个网格，仍然这样取搜索半径。
2多波束系统水深扫描
2.1多波束数字测深模型
如图2所示，设沿船体横向（垂直龙骨方向）第i个波束投射到海底到换能器所处方向的横向倾角为β，沿船体运动方向（平行龙骨方向）的纵向倾角为，换能器所在位置P。设船体从原点开始沿方向航行，航速V，采样时间T，由于存在横纵两个方向的波束宽度，第i个波束照射的海底离散点形成一个矩形点集Vi。
2.2测深搜索算法设计
设船体沿以角度东偏北45°，船速6m/s，匀速直线航行50s，采样时间1s。在进行地图扫描时，首先把换能器中心归结到最近网格点，以该点为中心，以2max|z|・tg（iθ）为边长建立正方形作为搜索区域，搜索点数为int+1，（grid为网格间距）；得到ri（xn，ym），从而得到Zi、Li，脚印位置（xi，yi），按照方向，以足够大矩形区域进行搜索，中间波束脚印位置进行最近网格点归结后作为向外边缘的下一个波束的搜索区域的起点。搜索完成后产生的16×51水深数据地图如图3所示，由该图可以发现，测船扫描过一段左低右高的慢坡曲面。水深在-980米～-1060米之间。
3水深测量校正
多波束系统的测深换能器发射的声波具有一个立体波束角，声波到达海底或碰到不同的介质界面时返射，再被换能器接收，接收换能器收到的信号可能是在这个波束角所覆盖的任何水深。正确深度应该是波束角所覆盖区域的中心点水深值，显然，它不同于声波覆盖区域的边缘点水深值，因此产生水深测量误差。此项误差在海底水深变化较大的地方尤为明显[4]，尤其在实测水深较深或海底坡度较大的地形时，由于存在波束宽度，往往产生地形图等深线向断面方向前后凸凹、弯曲的情形[8]，形成虚假地形，从而使绘制的海底地图与实际地形存在差异。在水深测量中，是以最短距离作为记录长度，由于存在波束宽度及地形的起伏，并且随着波束投射角的增大，测量的误差越来越大，尤其边缘波束更趋明显。在以往的实际工作中，一般认为约占10%的外侧波束（左右各占5%）定义为边缘波束，为了保证数据的精度满足要求，常常强制性地将这些波束舍弃[9]。为了减小测量误差，需要对测得的数据进行地形补偿。
在多波束水深测量中，由于存在横向波束宽度，随着如图2所示中α，β的影响。为了减小测量误差，需要对测得的数据进行地形补偿。以因子A1=cos（（i-1）θ-α）・sec（（i-）θ-α），A2=cos（α-iθ）・sec（α-（i-）θ），A3=sec（|（i-）θ-α|）来校正沿船体横向（垂直龙骨）的倾斜角；以B1=，B2=来校正沿船体纵向（平行龙骨）的倾斜角。校正前后的水深误差如图4（a）（b）所示，比较发现，校正后的水深误差极大减小，精度提高4倍左右。
4实时图重建
4.1时图重插值
对经过校正后的水深值，按照校正后的水深进行插值，确定插值区域的原则是，插值区域为实时扫描图的中间矩形区域，并且最大程度的覆盖扫描区。在扫描过程中，相邻波束脚印之间的距离往往较大，并且随着水深变化、地形起伏也不断变化，在建立实时扫描地图时需要在中间进行插值填补以连接起来，沿航行方向及其垂直方向建立规则格网，如图5中所示矩形区域，从换能器所在航线的中点开始，向左右、上下蔓延，为保持一致，每次蔓延的网格间距与初始插值时相同（实际插值时没此限制，可以根据需要设置）。使用ceil取整，以使得待插点都是整数个网格间距，同时又在插值区域内，其中，“×”为已知脚印点，格网的交叉点为待插值点。
使用本文1中的方法计算校正后水深的平均值和标准差，进行正态插值，建立新的水深实时测量图如图6（a）所示，相应的建图误差如图6（b）所示。
4.2重建误差分析
将地图重建后的水深值、水深测量值、测量校正值，分别与同位置原始水深进行比较，同位置原始水深值通过对初始10×10格网进行正态插值获得，它们之间的误差如表1所示，由表1发现，改正后的平均误差与最大误差都显著减小，地图重建后的平均误差及最大误差相比测量改正后的误差进一步减小，这主要是因为在选取插值区域时选择了扫描地区的中间矩形区域，边缘波束被舍掉，而边缘波束测量及改正后的误差都在整个扫描图中最大所致。
5结论
本文通过正态权重的距离反比权重插值算法，对某海域的已有实际水深数据进行插值，继而使用多波束进行扫描、校正，完成实时地图的重建。海底地形变化多样，细节复杂，使用5个因子的校正是近似的。同时，为了航行安全起见，可以保留适度的浅误差，但如果取得的实时图是为了匹配，则要提高精度，要根据实际情况，二者兼顾。
【参考文献】
[1]马建林，金菁，来向华.多波束测深海底数字地形模型的建立.海洋测绘，2005，25（5）：15-17.
[2]陆丹，李海森，魏玉阔，邓平.多波束测深系统中的海底地形可视化技术研究[J].仪器仪表学报，2012，33（2）：450-456.
[3]何磊，袁德奎，陶建华.一种适用波浪潮流计算的地形插值算法[J].港工技术，2004，（4）：4-6.
[4]李章林，张夏林.距离平方反比法矿产资源储量计算模块设计与实现[J].地质与勘探，2007，43（6）：92-97.
[5]靳国栋，刘衍聪，牛文杰.距离加权反比插值法和克里金插值法的比较[J].长春工业大学学报，2003，24（3）：53-57.
[6]刘晓萌，常占强.数字高程模型数据分布规律的研究[J].河北师范大学学报（自然科学版），2006，30（4）：473-477.
[7]高佩兰，关永贤.多波束海底地形图的制作方法探讨[J].海洋测绘，2012，32（4）：71-73.
关键词：地面沉降；信息化管理；数据预测；空间插值；天地图
目前，中国19个省份中超过50个城市发生了不同程度的地面沉降，重灾区主要为长江三角洲地区、华北平原和汾渭盆地这3个区域[1]。地面沉降作为一种累进性的地质灾害，影响范围广，防治难度大，已成为困扰城市建设、影响居民安全生活的一种重要地质灾害，并且严重影响着城市地区的经济增长和可持续发展[2]。为了有效地控制城市地面沉降，杭嘉湖地区沿海县市政府城市规划、地质等管理部门每年都投入大量人力和财力对地面沉降及地下水位进行不间断的动态监测，并积累了大量丰富的观测数据，为城市规划和建设发展提供了重要的科学依据。然而由于历史条件或技术原因的限制，这些大量的宝贵资料都是通过手工整理，工作效率低，查阅资料也极为不便，大批的数据没有得到科学的管理、分析和有效的利用。如何有效管理这些数据，并进行合理分析，进而掌握沿海县市的地面沉降变化规律，分析地面沉降机理，对防震减灾、城市规划建设等都具有十分重要的意义。
1需求分析
为科学地管理多年历史沉降监测资料，建立一套合理科学的基于GIS的沉降信息管理和共享服务平台，实现不同时期、不同观测点的多源多时相沉降观测数据的自动入库、数值分析、可视化显示、趋势分析与预测、各类报表、图形的自动生成与输出，及时快速地掌握地面沉降变化规律，提供处理措施，为城市规划、城市建设与管理及社会公众信息需求提供服务，提高城市规划和城市空间的管理开发与利用水平，笔者所在单位组织相关技术力量进行了沿海县市级城区地面沉降信息管理系统的研究与开发，并推广至杭嘉湖多个县市，取得了很好的用户评价。
2系统设计
2.1系统总体框架
本系统采用B/S三层架构，PLSQLDeveloper进行Oracle11g数据库管理，IIS网站，实现对地面沉降数据的管理、统计查询分析和预测分析等。系统的技术框架如图1所示。
2.2关键技术
2.2.1数据库设计
基于B/S框架的沉降信息的管理与可视化系统设计中，数据库设计是整个设计中较为重要的部分之一，其设计的好坏对系统性能有很大影响。本系统运用Oracle实现与沉降相关数据交换与存储。与沉降相关的数据主要有站点的基本数据及站点的多期时间序列沉降数据，在该系统的数据库设计中主要分沉降站点基本数据表和多期历史沉降数据表，相关E-R如图2所示。系统根据该E-R图进行了沉降数据相关表(沉降观测点表、历年高程记录表)、底图配置表、用户表、图片或文件与观测点对应表等详细的数据表设计。
2.2.2沉降数据预测方法研究
利用数学模型将多期观测的地面沉降量数据拟合成曲线，进而计算后期沉降的方法在国内外得到了广泛的应用。这类方法包括双曲线法、三点法、Asaoka法、泊松曲线法、灰色理论和人工神经网络法等。各种预测模型的适用性不同，如双曲线法适合长期预测，时间越长，准确度越高[3]；而Asaoka法的优点在于它可利用较短时间内的观测资料，就能得到较为可靠的预测结果[4]；灰色模型要求原始数据变化较为平缓[5]。城市范围内地面观测点较多，分布范围广，地质条件不尽相同；任何一种地面沉降模型都是在一定地质条件下推导出来的，在城市沉降预测分析中，不可能预先为观测点选择适用的模型，也不可能设计一种通用的地面沉降模型[6]。因此，沉降预测模型计算出结果后，还需对结果进行精度验证，可以选用平均绝对百分误差(MAPE)进行模型精度评价[1]。通过这种精度评价可以得出的结论是：越多的原始序列期数越能反映观测点的沉降趋势；在实际应用中，应优先信任短期预测结果；可在观测期数5、6期以上时优先采用灰度模型；在缺少数据样本的支持下，应当优先采用Asaoka法该法进行预测；三点法在当前样本下比较适合长期预测；双曲线法适应性较强，其预测精度不高；观测点在不同预测模型的预测精度也不相同，可以对模型进行改进。
2.2.3天地图与沉降数据融合
地面沉降信息管理系统包括多期沉降观测数据的整理与和建库、沉降观测数据的空间可视化与查询，以及沉降观测数据时空分析，为城市地面沉降的预测、防控、减灾提供支持。沉降观测数据通常为单一离散观测点，因此，对沉降观测信息的管理、分析和预测需要城市级基础地理信息数据支撑。本系统提出以采用智慧城市地理信息共享平台架构如县市天地图作为城市离散沉降信息管理的空间展示基础。
2.2.4基于交叉验证的空间插值方法
在城市沉降监测中，常用的单点测量手段获取的单站观测数据是离散的、单一的。由于观测站点个数受到限制，获取的观测资料只能反映局部沉降监测信息，无法从整体上对观测区进行监测。为了弥补这一不足，近年来，多个专家学者利用空间插值的方法进行了地面沉降由点到面的监测研究，并在城市地面沉降观测中得到了较多应用。常见的空间插值方法有：IDW法、Kriging和样条函数方法[7]。对于不同的观测点分布、不同的地形环境和观测数据量，其插值计算结果在一定范围内存在差异。IDW插值把估算值和实测点间的距离作为权重因子，权重值为距离反比[8-9]；Kriging插值是通过一组具有z值的分散点生成估计表面的高级地统计过程[10-11]；样条函数插值是通过控制估计方差，利用一些特征节点，用多项式拟合的方法来产生平滑的插值曲线[12-14]。在城市地面沉降插值计算时，3种空间插值模型表现出各自的特点[15]：IDW法简单易行，效率高，但是易受极值影响；Kriging法考虑区域空间特性，但是算法复杂度高并且耗时；样条函数法保凸性、逼真性和平滑性好，但是误差估计困难要求点数量足够且分布均匀。
进行空间插值后需要对插值结果进行精度验证。首先假定某些观测点高程值未知，通过周围站点实测值插值来估算，然后计算验证站点的实测值与估算值之间的误差值。被估算高程值的站点称为验证站点，用于插值的站点称为采样站点。常用的评估方法有平均绝对误差(MAE)和平均误差平方的平方根(RMSIE)。采用这两种评估标准能很好地反映不同插值方法的相对精度，以此选择最适宜的空间插值方法进行城市沉降计算。
2.2.5沉降数据导出工具设计
导出工具用于将沉降趋势预测结果转换为天地图支持的格式。它是基于VisualStudio2015集成开发环境，采用C#语言，利用二次开发平台ArcGISEngine10.1进行开发的。导出工具的地图浏览功能基于桐乡地图的MXD格式的地图文档；数据库使用沉降信息管理系统的数据库，通过服务器的数据库服务进行通信。导出工具根据用户选择的条件对观测点进行沉降趋势预测并将结果以专题图的形式展现在地图中，且可以生成相应的图例。最后可以将生成的展示结果以MXD的格式保存至本地硬盘。其工作流程如图3所示。
2.3功能模块
本系统具体的功能模块包括：
(1)观测点与高程信息管理。存放观测点的位置信息和高程信息，是本系统进行沉降统计与分析的基础数据。
(2)文档管理。存放与观测点相关的信息，如观测点的照片信息和点之记等。文档管理可以按照文档内容的性质对文档进行分门别类，提高检索的便捷度。支持上传、编辑、下载、标注、条件查询等文件管理操作。
(3)底图管理。用户可以增加、删除、修改底图信息，切换使用不同的底图。
(4)数据导入与导出。数据导入可以随时更新系统基础数据，使得系统能够存储更多观测点的历年高程信息，导入沉降数据时还可以验证其有效性、正确性；数据导出可以选择一段数据，得到观测点的位置信息、高程信息及沉降统计信息等。
(5)用户管理。系统采用的用户分为一般用户、数据管理员及超级管理员等3种角色。不同类型的用户分别有不同的权限：一般用户只有查看预览与检索查询沉降信息的功能；数据管理员不但有预览沉降信息的功能，而且还有增加、编辑、底图配置、统计分析、预测分析、InSAR数据分析等权限；超级管理员具有系统的最高权限，除了有数据管理员的权限之外，还有对用户管理的权限。
(6)报告导出。用户可以根据需要，选择任一时间段的控制点信息，生成包含有趋势预测信息、沉降等值线信息、DEM变化信息、淹没信息、沉降量统计信息的报告并可以打印。
(7)空间查询和模糊查询。空间查询指用户可以在地图上框选观测点，查询其高程信息和属性信息。模糊查询主要是指用户可以根据观测点中名字的关键字查询到观测点，并在底图上高亮显示。
(8)图层管理。用户可以根据自己的需要选择哪些观测点可以出现在底图上，哪些工具需要被隐藏。
(9)趋势预测。将双曲线法、三点法、Asaoka法及灰度模型等4种方法分别应用于已知观测年份的沉降数据预测未知年份的沉降数据，并将预测结果以折线图或沉降图的形式展示出来。
(10)沉降统计。根据数据库中已经存在的观测点数据，对某一时间段或所有观测点的沉降量进行统计，最后以折线图和柱状图的形式展现出来。
(11)沉降等值线。根据已知高程信息，得到沉降等值线，并在地图上展现出来。
(12)影像分析。影像分析主要是根据InSAR数据进行分析，分为沉降分析和淹没分析。将分析的结果进行渲染，并在地图上显示出来。
(13)决策分析。已知某地区内一些点的高程信息和沉降信息，提供一组预备建设的规划图进行检测(可以是点、线、面)，不同的建设物对其所占区域的高程信息和沉降信息等有不同的要求，要求对其中不适合进行建设的地方进行标注提示，以辅助最终决策的制定。
3系统实现与推广应用
本系统开发完成后在杭嘉湖沿海地区产生了积极的响应，如在桐乡、海宁和海盐、嘉善地区都得到了很好的应用。系统中沉降等值线功能、沉降点预测、沉降面预测分别如图4和图5所示。
1RGB到Lab颜色空间的转换
扫描仪RGB到Lab颜色空间的转换，经常采用的方法为三维插值法、BP神经网络法以及多项式回归算法[3]。三维插值法主要包括三次线性插值、三棱柱插值、金字塔插值和四面体插值[4]。三维插值法的主要原理是将源颜色空间分成多个三维几何体，再对三维几何体进行分析，若输入的待测值落在三维几何体内部，则采用三维几何体的栅格点来进行插值计算[5]。其优点是针对规则的几何体插值时，分割的块越小，精度越高，计算速度越快；缺点是对于不均匀的几何体，分割出的块不均匀，查找栅格点较困难，并且计算量较大。BP神经网络法通过对标准测量数据和输出数据的训练，得到待转换空间和目标转换空间之间的映射关系，实现空间转换[6]。该方法具有较强的灵活性及较高的精度，但其对训练方法与节点数选择，需要投入较多的人为计算，花费时间较长。多项式回归算法主要通过寻找源颜色空间和目标颜色空间数据之间的对应关系，计算数据拟合系数，建立一个多项式模型[7]。多项式回归法要求选取的样本点数大于多项式回归的项数[8]。该算法相对简单，可得到较好的转换精度，使用也较方便，但并不是说，多项式回归项数越多，其转换精度就越高，有时项数增多，转换精度提升的幅度也不大。由上文对表1的分析可知，扫描仪扫描后获得的建模数据R,G,B值无规律可循，故无法使用插值法将其分割成规则的立方体；由于IT8.7/2色靶总共288个，相对较少，若是采用BP神经网络法，则需要大量的数据进行训练；从转换精度上考虑，采用多项式回归法最为合适。因此，本研究采用多项式回归法，以实现扫描仪RGB到Lab颜色空间的色彩校正。
1.1IT8.7/2色靶单元多项式回归法
多项式回归算法的核心是求解系数，采用不同的多项式回归项数，其模型精度会有所变化，一般情况下，采用的多项式回归项数增加，其模型精度会有所提高[9]，但是当项数达到一定数量时，即使再增加项数，精度也几乎不会发生变化。根据建模数据，求得L,a,b3组多项式中的系数，然后，将系数代入公式（1）计算，求得对应的L,a,b值。本研究采用单元多项式回归法，基于多项式原理[10]，针对IT8.7/2色靶，并以5级分割RGB颜色空间为例，以实现扫描仪RGB到Lab颜色空间的转换。1）分割RGB颜色空间自定义RGB子空间查找表，将RGB空间按照5级均等分割，分割点分别为0,50,100,150,200,255，将RGB空间共分割成53=125个子空间，获得如下自定义RGB子空间查找表。本研究主要采用手动5级和自定义5级均等分割中心点，分别使用6项12点、6项30点、13项30点和13项50点，对扫描仪的建模数据和测试数据进行测试，并对采用不同方法得到的结果进行色差分析，找到最小色差对应的方法。综上分析可知，当采用不同方法分割中心点时，所得到的结果差异较大，经过比较，发现均等分割中心点得到的色差分布较好，因此，针对IT8.7/2色靶，采用5级自定义均等分割中心点，并采用13项50点进行多项式回归，此方法效果最好。此时，最大色差为6.25，最小色差为0.09，平均色差为1.22。
1.2自定义色靶单元多项式回归法
针对自定义色靶，建模数据为1331个色块，测试数据为512个色块，实验采用6级、8级和9级均等分割中心点，并分别采用16项130点、16项180点、19项130点和19项180点，对扫描仪的测试数据512个色块进行测试，对采用不同方法得到的结果进行色差分析，找到最小色差对应的方法。自定义色靶扫描仪RGB到Lab颜色空间的转换程序界面见图4。由扫描仪测试数据512个色块的测试结果可知，当采用8级均等分割中心点及采用19项180点进行单元多项式回归时，精度最高，效果最好。此时，最大色差为5.78，最小色差为0.16，平均色差为1.92。
2结语
关键词:DEM;土方量;精度分析;误差
中图分类号:TP11文献标识码:A文章编号:1009-8631(2010)04-0052-02
数字高程模型DEM是用于表示地面特征形态的多种信息空间分布的有序数值阵列。与传统的地形图件相比较,DEM作为地形的数字表达形式具有其特点:首先,DEM可以很方便地以各种比例尺和多种形式表达地形信息,如多种比例尺的地形图、纵横断面图及立体透视图、立体景观图等;其次DEM具有稳定的精度,不会因为时间的推移等种种因素而导致精度的损失;再次,因为是数字形式的,可以方便地增加和改变信息,便于存储管理、更新维护、传输共享、计算机自动处理,适合于各种定量分析与三维建模,DEM容易实现自动化和实时化。
当前利用全站仪、计算机和测图软件等先进软硬件设备自动生成DEM模型已广泛应用,利用DEM进行土方量计算是建设规划和许多工程项目中经常遇到的问题。
一、DEM土方量计算的精度分析
(一)DEM数据的精度分析
DEM数据的精度是多因素共同作用的结果,主要因素有:原始数据点精度、原始数据点密度、原始数据点分布、地表特征和建模方法等。对DEM数据的精度进行研究,应充分考虑上述主要因素的影响。DEM数据的精度涉及该数据的使用者和生产者,因此具有非常重要的意义,DEM数据精度的评估主要通过理论和实验研究的方式进行。
70年代起,DEM的研究方向从内插技术的发展转移到了对DEM精度的评估和控制方面,提出了数字离程模型精度分析的数学模型等。然而,对这些精度模型进行评估的数据比较单一的数据,对精度分析模型的检验缺乏大量空间数据的支持,因此其精度评估缺乏说服力。由于很难有一种模型能够完全适合复杂多变的地形表面,因此,在实际的应用中,针对各种复杂的地形表面,如何提高空间数据的精度也就显得非常迫切。
(二)土方量计算精度分析
1.原始数据密度对DEM计算体积精度的影响
数据密度是原始数据的一个重要属性,也是影响原始数据形成DEM模型精度的一个主要因素。因此,分析数据密度对体积计算精度的影响程度,如何在原始数据较少的情况下进一步提高。土方量计算的精度,其讨论结果可以推广到实际土方量的计算应用中。网格间距对不同地形特征的表面精度的影响趋势相似,但影响的程度却有差别,地形越平坦且粗糙度越小,网格间距对其精度影响的变化越缓慢。可见,实际采样过程中,在达到相同精度要求的条件下,较平坦且粗糙度较小的地表与较陡峭且粗糙度较大的地表相比,前者的采集间隔可以放宽,减少采样点的数据量。当原始数据密度非常大时,原始数据的精度是影响体积精度的主要因素。
2.特征点对DEM计算体积精度的影响
如果所选点邻域内的点高程均高于或低于所选点的高程,则选为特征点,可用下式表示:
Pi=(hi-h)/h
式中Pi是第i个邻域点的高程变化梯度,hi是第i个邻域点的高程,h是所选点的高程。可以根据Pi的符号进行特征点的判断,并依据Pi值的大小可以判断高程变化的梯度。
特征点是影响DEM进行体积计算的另一个主要因素,在原始数据较稀疏的情况下,有无特征点对DEM数据进行体积计算的精度影响很大。
我们通过反复计算得出:(1)增加特征点可以提高体积的计算精度;(2)原始数据的密度越稀疏,增加特征点对体积计算精度的提高越显著,随着原始数据密度的增加,提高的程度越小,到最后几乎没有影响;(3)添加特征点对体积计算精度的影响和原始数据网格的间距之间没有明显的变化规律可循,在数据密度较大时,特征点对精度的影响比较小,但当数据密度稀疏到一定的程度时(16cm的网格间距时),增加特征点对精度的影响变得显著。因此,在原始数据较稀疏的情况下可以使体积的计算达到很高的精度水平。
(三)不同内插方法的精度的分析
原始数据是影响误差的主要因素,但由内插方法产生的误差也是不容忽视的,因为在内插的过程中不仅有原始数据误差的传递。内插过程中本身也会产生一定的误差,不同的内插方法所产生的误差也是不一样的。对于相同的网格间距数据,在保持其它因素不变的情况,采用不同的内插方法构建DEM模型。在高程数据中,由于距离远的点权重虽然小,但高程值较大,也会造成很大的影响,通过计算我们得出的结论是:原始数据密度越大,精度越高。
二、DEM土方量计算误差探讨
理论上说,基于DEM的土方量计算法适用于任何情形,DEM的精度是影响土方量计算准确与否的主要因素。大量实验表明,由于实际地形的非平稳性,DEM的精度主要取决于原始采样点的密度和分布以及地形特征顾及与否。
(一)DEM土方量计算的误差来源
DEM误差的一种来源是对自然真实表面的采样过程。这种误差出现在从原始资料产生地面点的过程中,误差是由原始资料本身的不合适性和使用的仪器引起的。我们可以采取一些措施尽量减小它,使之达到误差允许的范围内。这些措施有:用仿射变换消除图纸变形引起的误差;采用高精度的数字化仪、扫描仪;将数字地图回放成纸张地图与出版地图进行比较等。
DEM误差的另一种来源是重新采样,即在保留了与原始地面较为逼近的情况下,将由原始数据派生的数据压缩成易于管理的过程中。这是一种性质不同的处理,因为提取的信息在很大程度上受采样区间和所用插值方法的影响。
DEM误差的表达和研究涉及到DEM与实际地面的平均偏差,也涉及误差的分布和误差的非随机空间分量。对于一个实际应用问题,要确定一个合适的DEM,则基本上依赖于研究对象所要求的精度,采样方法以及对地形变化的敏感度。
总之,DEM的精度取决于采样密度、测量误差(偶然误差,系统误差和粗差)、地形类别、高程点数目和位置等。影响精度的主要因素是数据获取,通过选择适当的内插法,可以获取基本相同的精度。
(二)DEM的模型误差
模型误差可定义为所建立的模型与客观现实之间的差异。可表示为:
F=MW
式中,F是模型的真误差;M是所利用的数学模型计算值;W是未知的客观现实真值,且M不等于W。由于客观现实是未知的,在实际的应用中,我们只能利用精度较高的数据来代替客观现实的真值进行误差分析。
DEM插值是DEM生产的核心问题,按插值点的分布范围,可将插值分为三类:分块插值、逐点插值和整体插值。而按二元函数附近数学面和参考点的关系,插值又可分为纯二维插值和曲面拟合插值两种。
每种内插方法都是关于如何估计内插值的一种假设,不同的内插方法可以产生对实际表面的不同的模拟结果,结果的可靠程度取决于建模的采样数据点的数量、分布和实际的地面状况。对于任何一种内插方法,数据点越多、分布越合理、实际地面越平滑,所产生的模型结果越可靠。
由于实际工作中,人力、物力以及财力的限制,不可能采集到高密度的采样点高程坐标数据,这就需要进行空间内插。DEM的空间内插,就是根据采样点数据(原始数据),按照数据分布规律的一种假设,构建样本点数据外推模型,生成一定密度的、均匀分布的高程数据集。如何选择、利用现有的内插模型生成拟合精度较高的大地曲面高程坐标数据集成为测量工作中的一个关键问题。如何评价多种内插模型的拟合精度,在实际工作中存在困难.因为大地曲面高程坐标数据没有真值,多种内插模型的拟合数据集就缺乏比较的基准,因而就无法判断内插模型精度的高低。
三、结论
基于DEM的土方计算法实现了土方量计算的计算机自动化,流行的GIS主流软件都可以生产不同地形类型的DEM产品,且精度都符合国家标准规定的精度要求,但不同的软件对不同的地形类型在精度上存在一定的差异。所以在计算土方量时,必须根据具体要求,选择合适的DEM产品。
经实例计算后发现:
(1)基于数字高程模型的土方计算方法理论严密,计算准确,若在数模生成时能给出土与石的分层、分类信息,则可真正实现土石方数量计算;
(2)DEM土方数量计算方法完全可以代替传统的计算方法,其实现不算复杂。
参考文献:
[1]黄杏元等.地理信息系统概论[M].北京:高等教育出版社,2008.
[2]张正禄等.工程测量学[M].武汉:武汉大学出版社,2005.
关键词Delaunay三角剖分质心原理三维插值矩阵变换
中图分类号：TP301.6文献标识码：ADOI：10.16400/ki.kjdks.2016.12.015
AbstractInthispaper，basedonthethree-dimensionalinterpolationandimagingofthree-dimensionalinterpolationproblems，modeling，solvingandrelatedanalysisandgivethetwomethodsandtheiradvantagesanddisadvantages.
KeywordsDelaunay；Triangularsubdivision；Centroidprinciple；threedimensionalinterpolation；Matrixtransformation
第一种方法是根据分段插值的原理，采用基于Delaunay三角剖分的三维插值算法进行求解；为了减小距离所求点较远的点插值时因增加插值次数而导致龙格现象，本文所给出的第二种方法是基于质心原理的三维插值算法，有效减小节点间的函数震荡，相较于普通的插值方法具有更好的收敛效果。
1问题重述
1.1问题背景
在实际的科研或工程研究中，常常需要在已有数据点的情况下获得这些数据点之间的中间点的数据，就需要使用不同的插值方法进行数据插值，而对数据的插值加密与成像问题，因为要对繁杂大量的数据进行插值处理以及涉及部分数字图像的相关处理，所以在生活以及军工业方面具有十分重要的现实意义与研究价值。
1.2待解决的问题
在以上提及的背景，有4851组相关数据有待解决：根据所给的4851组数据的特点，建立两种适当的数学模型，保证用该方法插值计算后电阻率数据的极值及其对应的位置不发生改变，并计算空间点所对应的电阻率数值。
2模型假设
（1）数据真实可靠；
（2）电阻率不随各类外部因素发生改变；
（3）三维物体具有各向同性。
3模型建立与求解
3.1问题分析
要求我们根据文本中所给的4851组数据，使用两种适当的插值计算方法计算出空间中某点处所对应的电阻率值并说明该方法插值后所得到的电阻率数据极值大小与位置不变。由于题中所给的条件限制，因极值的大小和位置不能发生改变，若直接采用普通的Lagrange、Newton和Hermite插值时，当所选取的节点个数较多，插值多项式的次数会相应提高，会出现龙格现象。因此针对这种情况，本文给出了两种方法求解。
3.2基于Delaunay三角剖分的三维插值计算方法
Delaunay三角剖分过程中采用Bowyer-Watson算法，其基本步骤为：（1）构造一个大三角形，包含所有散点，放入三角形链表；（2）将点集中的散点依次插入，在三角形链表中找出外接圆包含插入点的三角形（称为该点的影响三角形），删除影响三角形的公共边，将插入点都影响三角形的全部顶点连接起来，完成一个点在Delaunay三角形链表中的插入；（3）根据优化准则对局部新形成的三角形优化，将形成的三角形放入Delaunay三角形链表；（4）循环执行上述第2步，直到所有散点插入完毕。
以上计算过程可通过Matlab中的griddata函数所实现，其计算公式如下所示：
以下为方便叙述，使用“方法一”代指基于Delaunay三角剖分的三维插值算法。
3.3基于质心原理的三维插值计算方法
本方法首先选取所给数据中距离待求点（'，'，'）最近的八个节点及其函数值，如表1所示：
其中八个节点和是方向上与'的前后间距，和是方向上与'的前后间距，和是方向上与'的前后间距。由质心原理，设点（'，'，'）到各节点的权重因子分别为、、、、、、和。
以下为方便叙述，使用“方法二”代指基于质心原理的三维插值计算方法。
4模型评价与推广
4.1模型评价
4.1.1优点
（1）在方法一中，根据分段插值的原理，采用基于Delaunay三角剖分的三维插值算法进行求解，先将原始数据进行插值处理，再将循环所得的点所对应的坐标（或要求点的坐标）写入Matlab中，运用Matlab内置的griddata函数对各点进行插值插值运算。根据问题二的复杂性评估和定量分析易看出：方法一平均值与标准差的大小更接近原始数据且模型简单易懂，Matlab中有内置函数griddata便捷了此法的运算实现。问题四的相关数据更是能够清晰地反映方法一与原始数据的近似度较高；
（2）在方法二中，基于质心原理的三维插值算法，读入循环所所得的点所对应的坐标（或要求点的坐标），在待求点附近寻找距离其最近的8个节点，代入公式进行计算（通过Codeblocks编程实现）。根据问题二的复杂性评估和定量分析易看出：方法二运行时间较短，为方法一的近10倍，而且从数据结果可看出，方法二的精确度可达小数点后12位，较高。
4.1.2不足
（1）方法一运用Matlab实现算法，耗时较长，占用计算机内存较大，对计算机配置有一定要求，并且数据结果的精确度低于方法二；
（2）方法二运用Codeblocks实现算法，时间短、精确度高，但是与原始数据的相似性劣于方法一。
4.2模型的推广
采用图像插值技术来提高数字图像的分辨率和清晰度，从软件方面进行改进具有十分重要的实用价值，以及用于遥感图像的几何校正、医学城乡以及电影、电视和媒体广告等影像特技处理中。
应用插值加密的方法处理图像示例图如图2：上图为放大图片后的结果，易看出马赛克色块且清晰度较低；下图为插值加密处理后的图像，图片清晰度明显提高。
参考文献
[1]杨盛林，刘昱，刘玉峰.惯性平台热场分析及热设计的改进[J].中国惯性技术学报，2005.13（1）：5-9.
[2]王瑾.基于强对流数值模拟的贵州冰雹识别及临近预报方法研究[D].中国科学院博士论文，2009.
[3]肖艳娇.新一代天气雷达三维组网技术及应用研究[D].南京信息工程大学博士论文，2007.