吴自勇，蒋金豹，郭鉴威，王鑫达，季 杨

中国矿业大学(北京)地球科学与测绘工程学院，北京 100083

引 言

天然气是目前应用最为广泛的清洁能源之一。中国的天然气消费量由2001年的274亿m3增长至2017年的2 373亿m3，年均增长14.4%，增幅约为同期世界平均水平的7倍[1]。在此背景下，加快天然气地下储气库的建设成为必然选择，但是储气库在运行过程中因受到地质构造、腐蚀以及施工等多种因素的影响会发生泄漏，泄漏的天然气会通过土壤缓慢渗透到大气中，不仅会对地表土壤以及植被产生影响，而且当泄漏的天然气浓度达到一定程度后还会引发爆炸、人员冻伤、窒息等事故，对生命财产安全造成巨大损失。此外，天然气虽作为一种清洁能源，在未来仍有可能会成为温室效应的重要来源[2]。目前检测地下储气库泄漏的方法存在效率低、技术要求高、效果差的缺陷[3]，因此，迫切需要一种高效准确的检测方法。由于泄漏的天然气会致使地表相应位置的植被产生胁迫反应，因此，可通过探究植被的胁迫反应间接地识别地下储气库泄漏点，提高检测效率，节约检测成本，避免安全事故的发生。

目前，高光谱技术在天然气胁迫下植被的识别研究方面多集中于应用单一光谱信息进行目标识别[4-5]。Smith等[6]研究发现小麦、草地在天然气泄漏胁迫下一阶微分光谱在红边内逐渐减小，并伴有多峰与蓝移现象，并且利用725与702 nm一阶微分比值指数可以识别天然气胁迫下的草地、小麦和大豆； Noomen等[7]发现指数R440/R740能探测仅有的4个天然气泄漏点； Lassalle等[4]发现受油气胁迫61 d后的植被光谱曲线在可见光、近红外区域均有所升高。这些研究中仅利用植被的光谱信息进行识别，虽取得一定的效果，却忽略了空间信息。

成像光谱仪能够同时获取待检测对象的光谱和空间信息，这在天然气微泄漏点检测中具有重要意义，能够高效精确地检测出天然气泄漏点的空间位置。Noomen等[7]在美国Ojai地区发现天然气长期泄漏导致地表草地出现一个半径为20 m的环形裸土区，裸土区外则是茂密的草地，在考虑天然气泄漏胁迫植被症状空间特征的基础上，融合光谱信息成功探测天然气泄漏点。Werff等[8]发现地下天然气泄漏导致地表直径约6～14 m的环状区域内植被出现了明显的胁迫症状，并尝试利用Hough变换提取空间信息，融合光谱分类信息探测天然气泄漏胁迫植被。本文运用高光谱成像技术基于野外模拟实验检测天然气泄漏点。

1 实验部分

1.1 实验设计

实验区位于北京市大兴区(地理位置： 39°39′2.56″ N，116°34′33.10″ E)。实验区长约40 m，宽约20 m，区域内共设置24个边长2.5 m的方形地块，其中12个地块为通气胁迫组，剩余12个地块为控制组，地块之间间隔0.5 m，地块下埋藏有天然气管道，管道泄漏口在地块中心以下0.6 m处。自2018年10月1日起，通过天然气控制系统以1 L·min-1的速度连续不断通气，截止到2018年11月18日。实验期内每天进行巡视，保证植被不受其他因素(如缺肥、缺水)的影响。

1.2 数据采集

使用SOC710VP成像光谱仪在天气状况优良，当地时间10:00—14:00时间段内，获取各地块的高光谱影像。冬小麦影像数据采集时间分别为10月11日、10月24日、11月1日、11月9日和11月18日。模型验证所选的数据分别为大豆、玉米和草地，所选取的大豆影像为8月24日，玉米影像为8月7日，草地影像为9月17日。

1.3 影像处理

1.3.1 预处理

(1)标准反射率校正

SOC710VP采集的影像数据并非地物反射率数据，而是原始DN值，需经过辐射校正处理，将原始DN值转换为反射率数据，本文采用SOC710VP的配套软件SRAnal710对影像数据校正。

(2)光谱平滑

对校正后的影像使用五点移动均值平滑进行平滑处理，平滑公式为[9]

(1)

式(1)中，Rs为每个像元平滑后的光谱，Ri第i波段的反射率值，且i∈[3，k-2]k为高光谱图像的波段数。其他影像操作还包括影像剪裁。

1.3.2 后处理

(1)波段筛选

高光谱图像能够提供众多精细波段，但是信息冗余的问题也会凸显，因此有必要筛选出包含有用信息的特征波段[10]。蒋金豹等对高光谱数据进行ANOVA处理，选择敏感波段构建光谱指数来识别不同品质的建筑涂料。本文运用单因素方差分析(one-way ANOVA)评价胁迫与非胁迫小麦的波段可分性，选取p<0.05置信水平下的有效波段，去除可分性较差的波段[11]。

(2)纹理及颜色特征提取

天然气泄漏胁迫下作物长势及生理特征都会发生变化，这些变化不仅体现在光谱特征上还体现在纹理特征上，因此，可利用纹理特征辅助识别天然气胁迫下的植被[12]。基于统计学的灰度共生矩阵(GLCM)在纹理提取中应用较为广泛，本文采用GLCM计算纹理特征，灰度共生矩阵的距离参数设置为1，分别计算0°，45°，90°和135°方向上的同质性、对比度、非相似性、熵值、能量和相关性。此外，考虑到受天然气胁迫的植被会变黄，本文利用颜色矩计算RGB合成影像的前三阶颜色矩[13]，分别为一阶矩(Mean，μi)、二阶矩(Variance，σi)和三阶矩(Skewness，si)。

(2)

1.4 建模方法

1.4.1 归一化植被指数(NDVI)

NDVI是目前应用最广泛的植被指数之一，在已有研究中，该指数常被用来区分植被和裸土区域[14]。本文应用NDVI设置阈值提取胁迫下的裸露土壤，为防止阈值设置不合理影响后期数据处理，本实验通过多期影像进行结果对比，并最终设置阈值为0.4。

1.4.2 最小二乘支持向量机

支持向量机(SVM)在高光谱影像分类中应用广泛，最小二乘支持向量机(least squares support vector machine，LSSVM)作为一种新型支持向量机方法，用等式约束代替了传统支持向量机中的不等式约束，用误差平方替换了SVM中的松弛变量，通过求解一组等式方程得出最优分类超平面，从而避开了求解计算量相对繁重的二次规划问题，有效提高了求解速度[15]，LSSVM模型的输出可表示为

(3)

式中，ai对应xi的Lagrange乘子，对于核函数k(x)，本文选用高斯径向基(RBF)。

1.4.3 模型建立

本文首先基于NDVI对影像阈值分割，对于植被部分，提取光谱、纹理以及颜色特征，运用LSSVM识别胁迫下的植被区域，并与裸露土壤部分融合，得到识别结果，模型可简单表示如式(4)

(4)

其中，G(x)为胁迫区域，X为光谱、纹理以及颜色特征向量所构成的矩阵，SNDVI<0.4为裸露土壤区域。

为更好提取出天然气泄漏胁迫区域，本文运用图像形态学处理中腐蚀、膨胀和填充等算法对提取结果进行处理，并利用最小二乘拟合方法对识别的胁迫区域边界进行拟合。

2 结果与讨论

2.1 波段筛选与特征提取结果

对预处理后的影像进行ANOVA，在p<0.05水平下，选取F较大的波段，并剔除300 nm附近噪声较大的波段，优选出510，520，570，625，645，680和690 nm七个敏感波段，如图1所示。

由图2可知，胁迫第11天，受胁迫的小麦与为控制组小麦差异不显著。胁迫第24天之后，胁迫组小麦变黄，长势稀疏，胁迫组小麦光谱反射率在可见光范围内(370～760 nm)增大，而近红外范围(760～1 100 nm)降低； 胁迫组影像的同质性、能量和相关性三个纹理特征均低于控制组，对比度、非相似性和熵值均高于控制组； 胁迫组的影像的一阶矩和三阶矩明显高于控制组。随着胁迫的持续进行，受胁迫的小麦干枯、凋落，胁迫组小麦近红外区域光谱反射率仍低于控制组，但是与控制组差距减小。

图1 波段筛选Fig.1 Band Selection

2.2 模型识别

2.2.1 识别结果

针对植被的识别，将上述所提取的光谱、纹理以及颜色特征向量归一化作为LSSVM的输入，运用MATLAB编程，LSSVM中惩罚系数C=100，RBF参数γ=20，系数采用10-fold交叉验证得到； 对于裸露土壤，计算影像的NDVI，将NDVI值小于0.4的区域作为裸露土壤区域。将植被的识别结果与裸露土壤的识别结果融合，并对识别结果分别运用2×2模板进行腐蚀、膨胀处理，提取胁迫边界，得到如图3的结果。

图3 (a) 原始影像； (b) LSSVM识别的天然气胁迫下的小麦； (c) 裸露土壤； (d) 融合结果腐蚀膨胀操作后的结果； (e) 提取的胁迫边界和胁迫区域圆拟合的边界

2.2.2 结果分析

地块持续通气11 d后，胁迫组小麦仅有少量叶片边缘有微黄迹象，胁迫症状不明显，并且小麦生长初期易受到干旱、养分缺失等其他条件的相互作用，致使胁迫区与非胁迫区小麦在光谱、纹理以及颜色特征方面差异不显著，故模型未能较好地识别出胁迫第11天的区域。随胁迫的持续进行，胁迫中心区域小麦逐渐变黄，并且多数并在胁迫后期枯萎，与对照区相比，胁迫中心区小麦叶片萎缩，植株矮小、稀疏，模型均可较好地识别出胁迫区域，由图4的识别结果可以看出，胁迫发生第11天至第40天，胁迫区域逐渐变大，而第40天之后，胁迫区域趋于稳定，通过表1中对4个地块的胁迫区域半径统计发现，胁迫第40天至第49天，胁迫区域有缩小的趋势，通过进一步观测分析可以看出，在胁迫与非胁迫的交界处，随着小麦的生长，未受到胁迫的小麦或受胁迫较轻微的小麦叶片会逐渐伸长至胁迫区域，这种现象映射到冠层的影像中，则会致使胁迫区域的缩小。

2.3 讨论

2.3.1 模型验证

为验证本文所构建的模型能否适用于其他类型植被，本文运用构建的模型对同一试验区的大豆、玉米和草地进行分析，均得到较好地识别结果，如图5所示。

表1 胁迫天数与胁迫半径关系Table 1 Relationship between stress days and stress radius

2.3.2 绿晕现象

通过观察影像可以发现，在距胁迫区0～100像素(0.25 m)的范围内，小麦长势较边缘对照区域旺盛，这一现象与Noomen所发现的“green halo”现象较为一致，Noomen结合前人的研究推断油气泄漏会将水汽带到地表并且在晕圈附近聚集，同时泄漏的气体会使周围的气体分布发生变化，在两者的作用下，将会出现绿晕现象[7]。随着时间的推移，该区域的小麦叶片会向胁迫区伸进，这也进一步解释了胁迫后期胁迫区域有缩小趋势。

图5 模型验证Fig.5 Model Validation

图6 绿晕现象Fig.6 Green halo phenomenon

3 结 论

通过野外模拟天然气微泄漏对地表植被的胁迫实验，分析不同时间序列天然气微泄漏胁迫下冬小麦的高光谱影像数据，基于LSSVM并结合光谱、纹理以及颜色信息，构建天然气微泄漏识别模型，得到如下结论：

(1) 天然气泄漏胁迫下的冬小麦与对照区小麦相比，可见光范围光谱反射率上升，近红外区域下降； 通过方差分析选择510，520，570，625，645，680和690 nm七个可分性较好的特征波段，对特征波段分析可知，受胁迫的小麦影像同质性、能量和相关性三个纹理特征较对照区降低，而对比度、非相似性和熵值三个纹理特征较对照区提高； 胁迫下的冬小麦影像前三阶颜色矩总体上高于控制组。

(2) 本文的模型融合了光谱、纹理和颜色信息，并且在胁迫第24天有效识别出天然气泄漏区域。模型同样适用于大豆、玉米和草地，具有一定的适用性。

(3) 通过对胁迫区域圆拟合可知，冬小麦受天然气胁迫的范围呈现先扩大后缩小的时序特征。

本文构建的模型通过识别植被的变化能够较好的识别天然气微泄漏，但模型仍有不足之处，比如混合像元、阴影以及胁迫区域周围土壤的存在都会对模型的识别效果产生一定的影响。在接下来的研究中，将会结合天然气疑似泄漏区的多期影像、天然气管道埋藏图、地质图等资料解决这些问题，及时高效地识别出天然气泄漏点，为以后的工程应用提供理论基础和技术支持。