朱传林， 王学良， 柴 健， 杨仲江， 余田野

(1.湖北省防雷中心,湖北 武汉 430074; 2.南京信息工程大学 大气物理学院,江苏 南京 210044)

0 引 言

ArcGIS Engine是ESRI公司系列软件中的桌面GIS开发组件，ArcGIS Engine提供了空间分析的相关方法，在空间插值方面支持常用的反距离加权法插值、样条函数插值、克里金插值等插值方法。目前，一些学者已利用ArcEngine对降雨量的空间插值及等值线的绘制做了分析[1-4]，但国内利用ArcEngine自动绘制闪电密度等值线尚未见相关报道。闪电密度分布是闪电活动规律中最基本的参数之一，表征一定区域闪电频度的重要指标，快速、精确地绘制闪电密度等值线，可以有针对性地对闪电密度较大的地区加强雷电防护，具有指导意义。

与闪电强度等值线比较，闪电强度等值线的绘制相对简单，直接用ArcEngine组件进行空间插值，生成连续的闪电强度栅格图，再将闪电强度值相等的点连成强度等值线，而闪电密度数据是经过原始数据加工而来，如何科学地统计出闪电密度数据也是本文要解决的关键问题之一。闪电密度数据统计出来之后，再借鉴闪电强度等值线的绘制方法绘制闪电密度等值线。

在富客户端平台(Eclipse Rich Client Platform,RCP)之上，首先利用降维思想将闪电强度数据转化为闪电密度数据[5-6]，采用ArcEngine控件，基于已有数据进行空间插值，获得全局空间范围内各个点位的闪电密度数据，进而通过插值、等值线提取等方法，搭建了闪电密度等值线的自动绘制系统。

1 降维思想

闪电定位系统的原始数据包含的最基本的四要素有：时间、经度、纬度、强度。时间主要用于筛选出符合条件的闪电数据，在绘制二维等值线的过程中该要素是不用考虑的。因此，如何科学地利用闪电强度数据(包含经度、纬度、强度，如表1所示)转化为闪电密度数据(经度、纬度、密度)是绘制闪电密度等值线的关键环节。

表1 处理前的部分闪电强度数据

通常用网格的方法统计闪电密度，为了便于辨认及定位网格，在划分网格时，需要对网格各段的经线方向、纬线方向赋予特定代号(称之为经度代号、纬度代号)。某一条闪电属于哪一个网格需要与经纬度网格边界值做四次比较，定位当前闪电记录归属的网格后，赋予该条闪电网格的经度代号、纬度代号，然后，经度代号、纬度代号利用拼接的方式合并为一个代号，同时要保证其唯一性。经过循环程序,用以上方法判断以后，所有的原始数据都相应增加了三个属性，即经度代号、纬度代号、唯一代号。最后，统计唯一性代号的个数，再除以网格面积及闪电年数，便可得到密度参数，如表2所示。

从以上分析可以看出，降维思想统计密度数据时主要分为两个过程，①降维统计过程，从四次比较中提取出二维代码，二维代码组合为一维代码，进而对一维代码统计的过程；②密度参数赋值过程，查找一维代码的个数，除以相应网格的平方公里数及闪电数据的年限后，再赋予闪电原始数据密度参数值。

表2 降维程序处理后的部分闪电密度数据

注：C列即为经度代号，D列为纬度代号，E列为一维代码，F列为一维代码的个数，G列为密度参数

2 自动绘制闪电密度等值线的流程

闪电原始数据通常是txt文本格式，为了方便统计查询，需要将闪电原始数据经Java语言处理之后导入到Oracle数据库，然后利用SQL语句从数据库中查询出某区域的闪电强度数据，该数据包含经度、纬度、电流强度。利用降维思想将闪电强度数据转化为密度数据，这里的转化过程有两种方法：①利用Java语言实现降维思想的密度统计；②利用ArcEngine的Fishnet组件划分网格统计，再使用frequency函数实现一维代码的统计，两种方法实质上一脉相承的。

将统计的闪电密度数据加载成shp格式点图层，然后利用插值算法生成闪电密度栅格图，结合shp格式地图，用掩膜组件生成所需区域的闪电密度等值线。如图1所示为绘制闪电密度度等值线的技术流程图。其中：①利用Java语言对txt格式的闪电数据进行预处理，生成可以自动导入数据库的sql语句；②基于Eclipse RCP平台，采用Oracle数据库按省份建立了全国的闪电资料数据库，在数据的导入与输出的过程中要选择相应的省份，然后利用所研究区域的经纬度初步筛选出该地区的闪电数据；③利用降维思想统计密度数据时，根据统计区域合理的划分网格，不宜过大也不宜过小；④将统计出的闪电密度数据转换成ArcGIS支持的dbf文件，有经度、纬度、闪电密度3个字段；⑤将基数据加载到闪电密度数据图层；⑥选用插值方法，对关联后的闪电密度数据图层的闪电密度字段进行插值，生成闪电密度等值线[6-9]。

图1 闪电密度等值线自动绘制流程

3 插值算法

3.1 反距离加权插值法

反距离加权插值法(Inverse Distance Weighting,IDW)以数据点与网格点的距离平方倒数作为权值，距离越大权值越小，反之则越大，具体算法如下[11-13]：

设平面上分布一系列离散点，已知某坐标值为Xi，Yi，Zi(i=1,2,…，n)，依据周围离散点的值，通过距离加权值求z点值，则

IDW通过对邻近区域的每个采样点值平均运算获得内插单元值。IDW是一个均分过程，这一方法要求离散点均匀分布，并且密集程度足以满足在分析中反映局部变化。

反距离加权插值法优点是算法简单，容易实现，其插值的结果介于原始数据的最大值与最小值之间。其不足是没有考虑数据场在空间的分布，有时会因为原始数据的分布不均而导致估值产生偏差。所以，用这种插值结果绘出的等值线，平滑美观，但与实际可能有一定出入。

3.2 克立格插值法

克立格插值法有两个步骤：①对空间数据场进行结构分析，即，首先了解空间数据场的性质再提出变差函数模型；②在变差函数模型之上进行克立格计算。依据空间数据场是否存在漂移，将克立格插值法分为普通克立格插值法和泛克立格插值法。其中，普通克立格插值法常称为局部最优线性无偏估计，其算法如下[13-14]:

(3)

其中，λi为待定加权系数。该算法与其它各种内插法不同，克立格内插法是根据无偏估计和方差最小两项要求来确定上式中的加权系数λi的，故称为最优内插法。

当为无偏估计时，应满足,

(4)

(5)

克立格插值法考虑的因素较多，因此算法比较复杂，而且耗时较多；克立格插值过程中可能会产生负值，有时需要对负值区域进行人为干预。

3.3 样条函数插值法

样条函数插值是采用特殊分段多项式进行插值的形式，可以使用低阶多项式样条实现较小的插值误差，可以修改少数数据点配准而不必重新计算整条曲线。该插值算法的缺点是样条内插的误差不能直接估算，同时在实践中要解决的问题是样条块的定义以及如何在三维空间中将这些“块”拼成复杂曲面，又不引入原始曲面中所没有的异常现象等问题。

4 自动绘制闪电密度等值线的关键步骤

4.1 Java语言实现降维统计闪电密度数据

首先，根据区域的最大、最小经纬度，获取经纬度的跨度，其代码如下：

float longitudeSpan = maxlongitude-minlongitude;

float latitudeSpan= maxlatitude-minlatitude;

然后，利用经纬度的跨度以及网格宽度，得到经度、纬度方向的分段总数，其关键代码如下：

int longitudetotalNUM= new

BigDecimal(longitudeSpan/(unit*getLatitudeChange())+"").intValue();

int latitudetotalNUM= new

BigDecimal(latitudeSpan/(unit*getLongitudeChange())+"").intValue();

其中，getLatitudeChange()、getLongitudeChange()两个方法是用来获取单位公里数对应经纬度的变化量。

接着，使用循环语句来判断，当前某一条闪电对应的经度代号及纬度代号，并将其结合成一维代码，为保证一维代码的唯一性，需要对纬度代号的位数进行适当判断，采取经纬代号乘以10的n次方的形式串接，这里的n是指纬度代号的位数，部分代码如下：

if(latitudetotalNUM>=10&&latitudetotalNUM<100) {onecode =i*100+j; }

最后，将一维代码的个数与网格公里数、闪电数据年限做除法运算，得到密度参数值，并把最后的数据写入到文件，其部分代码如下：

double density=getDensity(countnumber,unit,year);

output.println(longitude+","+latitude","+ longitudecode+","+ latitudecode +","+onecode+","+density);

4.2 ArcEngine组件实现降维统计密度数据

(1) 闪电原始数据加载生成shp文件。将统计出的原始闪电数据加载生成shp格式点图层，在点数据的加载过程中主要用IFeatureCursor、IFeatureBuffer、IFeature、BufferedReader等接口，首先利用IFeature Cursor、IFeatureBuffer接口获取添加图层缓冲区，然后用BufferedReader从文件中读取数据，向IFields实例化对象中添加数据。其部分关键代码如下：

IFeatureCursor featureCursor = featureClass.IFeatureClass_insert(true);

IFeatureBuffer featureBuffer = featureClass.createFeatureBuffer();

IFields fields = featureBuffer.getFields();

featureCursor.insertFeature(featureBuffer);

featureCursor.flush();

(2) 降维统计密度数据。首先用CreateFishnet接口根据研究的区域范围划分公里网格，该图层与闪电原始数据图层关联后便获得了一维代号，进而利用Frequency接口统计一维代号的个数，然后CalculateField接口实现密度参数的计算，并添加密度参数字段，最后将计算的密度参数结果赋予该字段。其部分关键代码如下：

createFishnet.setOriginCoord(originpoint);

frequency.setFrequencyFields("FID_onecode");

calculate.setExpression("[FREQUENCY]/"+(nuit*nuit*year)+"");

4.3 密度等值线的自动绘制

(1) 闪电强度数据的插值。以反距离加权平均插值为例。利用IInterpolationOp接口的IDW方法实现反距离加权插值，关键代码如下，其中参数inFeatures为具有密度参数字段的shp文件的路径；Output为输出路径；densityValue为闪电密度值，CellSize参数有图片的大小来决定，power为幂参数，要求为正实数(默认值为2)。

interpolationIDW.setInPointFeatures(inFeatures);

interpolationIDW.setOutRaster(Output);

interpolationIDW.setZField("densityValue");

interpolationIDW.setCellSize(new Double(LightningData2Image.cellsize));

interpolationIDW.setPower(2);

(2) 等值线的提取。使用Contour接口从插值后的栅格文件提取等值线，关键代码如下。其中OutPolyinFeatures为插值后的闪电密度栅格图，ContourInterval为密度等值线间隔，BaseContour为绘制等值线的起点。

contour.setInRaster(inFeatures);

contour.setContourInterval(10);

contour.setOutPolylineFeatures(Output);

contourtool.setBaseContour(0);

(3) 等值线的出图。等值线的图形输出部分使用到的接口有PageLayoutControl、IMapGrid、IGridLabel等，IMapGrid用来添加经纬度网格，IGridLabel用来添加文字等标签。部分关键代码如下，其中的参数densityRect用来控制输出图片的大小，pageLayoutControl用于获取当前视图，export对象实现闪电密度等值线图形的输出。

IEnvelope pEnvelope = new Envelope();

pEnvelope.putCoords(densityRect.left densityRect.bottom, densityRect.right densityRect.top);

pageLayoutControl.getActiveView().output(export.startExporting(), (int) export.getResolution(), densityRect,pageLayoutControl.getActiveView().getExtent(), pageCancle);

export.finishExporting();

5 实验结果与分析

闪电密度是雷电灾害风险评估中重要的量化指标之一，为了更好地实现闪电密度等值线的自动绘制，本文采用Eclipse RCP技术搭建了自动绘制闪电密度等值线的软件平台，选用了跨平台的Java语言编写。该软件选用Oracle数据库建立了全国的闪电资料数据库，txt格式的原始闪电数据可以自动入库，同时需要注意的是需要选择与数据相应的省份以及闪电定位仪型号(如图2所示)，以保证数据的正常导入。考虑到闪电数据入库是一个相当耗时的过程，加入了Job后台任务并发处理机制，数据在入库的过程中不会出现“假死机”现象，同时也不影响用户的其他操作。当数据入库结束时，要用对话框提醒用户，因为在非UI线程内不能直接进行UI处理，采用了异步处理机制，即使用Display.getDefault().asyncExec(new Runnable())语句进行UI弹出对话框，给使用者良好的体验效果[15-16]。

图2 闪电数据入库

本实验以武汉地区为例绘制武汉地区的闪电密度等值线，首先将湖北省的闪电原始数据导入库，使用闪电数据库查询功能可以查询数据的起止时间以及根据起止时间查询数据量等，如图3所示。在某些特殊情况下，如：由于误操作导致数据重复导入了，可以利用闪电数据库的管理功能将该省份的数据删除，然后再重新入库。如图4所示。

图3 闪电数据的查询

绘制武汉地区的闪电密度等值线时，系统根据所选择的起止时间查询出该地区的所有闪电数据(主要包括闪电时间、经纬度、强度等要素)并将其写入到指定路径的txt文件，通过预先设定的降维统计密度数据的方法，系统会自动选择“Java语言”或是“ArcEngine组件”进行密度数据统计，将得到的包含密度参数的shp文件通过反距离加权插值方法自动绘制出闪电密度栅格图(如图5所示)，再用Contour组件从密度栅格文件中提取出闪电密度等值线，图6为其局部放大平滑后的图，然后利用武汉地区的地图掩膜出该区域的等值线，通过掩膜等工具就可以把研究区域以外的数据剔除，只保留地图之内的等值线数据。图形输出时需要通过pageLayoutControl获取当前活动的视图，输出的图形支持多种格式，如jpg、tiff、emf、gif等，值得一提的是在图形输出时设置了一些手动选项，用户可以根据个人需要进行定制图名、图例、单位、指北针、比例尺，及调整网格大小等。

图4 闪电数据的管理

图5 闪电密度栅格图

图6 闪电密度等值线

6 结 语

本文运用Oracle数据库统计出武汉地区的闪电数据，利用降维思想将闪电原始数据转化为闪电密度数据，结合ArcGIS Engine提供的空间插值方法分析了闪电密度等值线的自动绘制步骤并给出了关键步骤的主要程序代码，设计实现了闪电密度等值线的自动绘制模块，该模块在武汉地区地理信息系统中得到应用并取得较好的效果。

在实际的应用中，如果没有shp格式地图，等值线区域将是一个矩形区域。为了更好地结合地图绘制等值线，可以结合Google Earth做进一步的开发，把Google Earth下载的影像图进行矢量化，进而可定制特定区域的闪电密度等值线。

[1] 黄丙湖，孙根云.雨量线自动提取的分析与实现[J].计算机工程与设计,2010,31(15):3499-3502.

HUANG Bing-hu,SUN Gen-yun. Analysis and realization of automatic extracting isohyets[J]. Computer Engineering and Design,2010,31(15)：3499-3502.

[2] 彭思岭，邓 敏，黄丙湖，等. 基于ArcEngine的等雨量线绘制[J].地理信息世界,2010(2):79-83.

PENG Si-ling, DENG Min, HUANG Bing-hu,etal. Drawing Rain Isolines Based on ArcEngine[J].Geomatics World,2010(2)：79-83.

[3] 傅希德，唐 俊，张晓盼，等.基于Surfer和ArcGIS Engine的雨量等值线自动生成法[J].水电能源科学,2008, 26(6):9-10.

FU Xi-de, TANG Jun, ZHANG Xiao-pan,etal. Approach on Automatic Drawing of Rainfall Isoline Based on Surfer and ArcGIS Engine[J]. Water Resources and Power,2008,26(6)：9-10.

[4] 宋丽琼,田 原,邬 伦,等.日降水量的空间插值方法与应用对比分析—以深圳市为例[J].地球信息科学,2008,10(5):566-572.

SONG Li-qiong, TIAN yuan, WU lun,etal. On Comparison of Spatial Interpolation Methods of Daily Rainfall Data:A Case Study of Shenzhen[J]. Geo-information Science,2008,10(5):566-572.

[5] 朱传林，王学良，杨仲江，等．降维思想在统计雷击大地密度中的应用[J]．气象科技，2012，40(5)：839-842．

ZHU Chuan-lin, WANG Xue-liang, YANG Zhong-jiang,etal.Application of Dimension Reduction Concept iin Computing Ground Flash Density[J]. Meteorological Science and Technology,2012,40(5):839-842.

[6] 朱传林，范宏飞，段振中，等．基于ArcEngine的闪电强度等值线的自动绘制[J]．湖北大学学报(自然科学版)，2012，34(3)：368-372．

ZHU Chuan-lin, FAN Hong-fei, Duan Zhen-zhong,etal.Lightning Intensity Automatic Contour Drawing Based on ArcEngine[J]. Journal of Hubei University(Natural Science),2012,34(3):368-372.

[7] 汤国安,杨 昕.ArcGIS地理信息系统空间分析实验教程[M].北京：科学出版社,2006.

[8] 刘旭林，赵文芳.气象观测数据等值线自动绘制系统[J].气象,2009, 35(4):102-107.

LIU Xu-lin, ZHAO Wen-fang. Automated Contour System for Meteorological Observation Data[J]. Meteorological Monthly,2009,35(4)：102-107.

[9] 韩 鹏,王 泉，王 鹏，等. 地理信息系统开发—ArcEngine方法[M].武汉：武汉大学出版社,2008.

[10] 方书敏,钱正堂,李远平.甘肃省降水的空间内插方法比较[J].干旱区资源与环境,2005,19(3):47-50.

FANG Shu-min, QIAN Zheng-tang, LI Yuan-ping. Comparison of Four Precipitation Spatial Interpolation Methods in Gansu[J]. Journal of Arid Land Resources & Environment,2005,19(3):47-50.

[11] 孟庆香,刘国彬,杨勤科.黄土高原降水量的空间插值方法研究[J].西北农林科技大学学报(自然科学版),2006,34(3):83-88.

MENG Qing-xiang, LIU Guo-bin,YANG Qing-ke. Research on Spatial Interpolation Methods of Precipitation on Loess Plateau[J].Journal of Northwest Sci-Tech University of Agriculture and Forestry(Natural Science Edition), 2006,34(3):83-88.

[12] 秦 涛，付宗堂. ArcGIS中几种空间内插方法的比较[J].物探化探计算技术,2007,29(1):72-75.

QIN Tao, FU Zong-tang. Comparison of Several Spatial Interpolation Methods in ArcGIS[J]. Computing Techniques for Geophysical and Geochemical Exploration,2007,29(1):72-75.

[13] 李 新,程国栋,卢 玲.空间内插方法比较[J].地球科学进展,2000,15(3):260-265.

LI Xin, CHENG Guo-dong, LU ling. Comparison of Spatial Interpolation Methods[J]. Advance in Earth Sciences, 2000,15(3):260-265.

[14] 朱求安,张万昌,余钧辉.基于GIS的空间插值方法研究[J].江西师范大学学报(自然科学版),2004,28(2):183-188.

ZHU Qi-uan,ZHANG Wan-chang,YU Jun-hui. The Spatial Interpolations in GIS[J]. Journal of Jiangxi Normal University (Natural Sciences Edition), 2004,28(2):183-188.

[15] 陈 冈.Eclipse RCP应用系统开发方法与实践[M]. 北京：电子工业出版社.2007.

[16] 赵满来.可视化Java GUI程序设计—Eclipse VE开发环境[M]. 北京：清华大学出版社.2010.