基于热红外和夜间灯光遥感的GDP空间化研究——以福建省为例

张雪红 吴雨阳 王永俊 甄晓菊 孙艺 薛青宇 刘开颜

0 引言

国内生产总值(GDP)是衡量区域经济发展水平的一个重要的综合性统计指标．传统以各级行政单元统计的GDP资料无法显示统计区域内部的GDP差异,解决此问题的有效方法之一是实现 GDP 的空间化[1]．

早期GDP空间化方法主要利用地理信息系统(GIS)技术并结合土地利用类型或者社会统计数据的人均GDP数据进行研究[2-3]．随着夜间灯光遥感技术的发展,国内外专家发现夜间灯光遥感影像可以客观、实时地记录城市地物自身发出的灯光强度,并由此尝试利用夜间灯光数据进行城市土地覆盖、人口估算、经济活动等城市化关键要素的快速评估与空间综合分析[4]．Elvidge等[5]指出了夜间灯光数据与GDP指标之间存在较强的线性关系．Tilottama等[6]通过分析夜间灯光影像数据和校正后官方经济数据之间的空间关系建立了回归模型,从此拉开了通过建模反演GDP空间化的序幕．在此之后,韩向娣等[7]首次提出了结合DMSP/OLS与土地利用数据建立GDP空间化模型的方法,从而提高了反演模型的精度．不同产业所占GDP比重能有效反映区域内产业结构特征和产业布局,若实现单一产业的精细化空间分布,将对区域产业结构优化、城乡规划乃至精准扶贫等提供重要的决策依据．而夜间灯光遥感数据难以区分GDP中第二和第三产业[7],以往研究中也鲜见报道．

此外,GDP空间化研究中夜间灯光遥感数据主要使用DMSP/OLS、NPP-VIIRS等,但其中低空间分辨率、灯光过饱和等因素在很大程度上制约着遥感模型的精度．虽然部分研究针对以上问题提出了一些解决方案,如王俊华[8]在对四川省GDP空间化的研究中结合DMSP/OLS和NPP/VIIRS两种数据连接使用的方法分析了不同时相的相对误差,李峰等[9]采用土地利用类型数据与DMSP/OLS数据相结合的方法消除后者本身存在的灯光过饱和以及溢出现象所导致的噪声误差,但这些手段和方法无法完全消除以上制约因素所引起的不确定性．

相关研究表明,能源消耗及碳排放与城市化有较强的相关关系[10-13],城市空间形态能够影响因城市热岛产生的碳排放[14]．碳排放会导致在人类涉及的活动区域地表温度相对自然区域温度更高,因此地表温度信息在一定程度上可以反映以碳排放以及能源消耗为显著特征的第二产业空间分布[15]．为了评价结合热红外遥感数据进行不同产业空间化的可行性,本文以福建省为研究区,基于“珞珈一号”科学实验卫星夜间灯光遥感数据和Landsat8热红外遥感数据,并结合土地利用信息,尝试实现对福建省第二、第三产业增加值的空间化．并进一步以厦门市为例,从局部分别对两种方法的第二、第三产业增加值空间化结果进行对比,以精细化评价本文提出的第二、第三产业增加值空间化方法的可行性及有效性．

1 研究区概况与研究方法

1.1 研究区概况

福建省位于我国东南沿海(115°50′～120°40′E,23°30′～28°22′N),其森林覆盖率为65.95%,居全国首位．福建省区域特色经济发达,涉及制造、建筑、加工、服务等行业．2018年福建省实现地区生产总值35 804.04亿元,比上年增长8.3%．其中,第二产业增加值17 232.36亿元,增长8.5%，第三产业增加值16 191.86亿元,增长8.8%．第二和第三产业增加值比重分别为48.1%、45.2%．福建省第二产业比重与全国相比偏高,并长期处于第二产业占比高于第三产业的状态[16]．福建省厦门市、泉州市以及福州市的第二、第三产业较为发达,产业增加值明显．

1.2 数据来源

1.2.1 遥感数据

1)夜间灯光数据

夜间灯光信息由“珞珈一号”遥感影像数据提供,其影像幅宽达到260 km,空间分辨率为130 m[11-12]．本文采用2018年拍摄的福建省的“珞珈一号”01星遥感影像,数据来源为高分湖北中心提供的全国灯光数据图(http:∥www.hbeos.org.cn),所获取数据已经进行过辐射校正、大气校正以及几何纠正等预处理工作．

2)Landsat8影像数据

Landsat8卫星上携带两个传感器,分别是陆地成像仪(Operational Land Imager,OLI)和热红外传感器(Thermal Infrared Sensor,TIRS)．OLI有9个波段,可用于后续的土地利用信息提取与归一化差值植被指数(NDVI)计算．TIRS主要用于收集地球两个热区地带的热量流失．由于福建地处东南沿海,Landsat8影像云覆盖比例较大,经过筛选最终选取2017—2019年同一季节的12景Landsat8影像用于本文研究(通常短期内土地利用类型、工矿用地等变化较小,因此选取2018年前后相邻两年同季节Landsat8影像对本文研究结果影响不大),其影像成像时间和条带号如表1所示．数据下载于http:∥ids.ceode.ac.cn/．

表1 Landsat8遥感数据

1.2.2 土地利用数据

土地利用数据为中国科学院地理科学与资源研究所发布的1∶10万土地利用遥感监测数据．本文所使用的土地利用数据是在2015年土地利用遥感监测数据的基础上,基于表1中所选取的Landsat8/OLI遥感影像通过人工目视解译更新而成的．由于本文主要研究福建省第二、第三产业经济空间化,故选取土地利用类型数据中的“城镇用地”、“农村居民点”和“工矿用地”的二级土地利用类型,并与其对应的夜间灯光遥感数据和反演的地表温度数据利用Arcgis分别进行空间叠加分析．

1.2.3 GDP统计数据

本文使用的社会经济资料主要来自于《福建统计年鉴(2019)》,它收录了2018年福建省全省及各地区、各部门经济和社会发展各方面的统计数据．

1.3 数据处理及研究方法

在前人的GDP空间化研究中,由于夜间灯光数据对第二、第三产业本身的不可区分性,通常将两者综合起来进行分析,且大多采用单一的夜间灯光遥感数据或同时结合土地利用数据进行建模．考虑到热红外遥感数据可以有效反映第二产业分布信息,本文在土地利用数据和夜间灯光遥感数据的基础上增加热红外遥感数据,以此来评价基于遥感技术对第二、第三产业分开进行空间化的可行性．为了后续建模的需要,将夜间灯光遥感数据和热红外遥感数据重采样为与土地利用类型对应的30 m空间分辨率．此外,福建省共有64个县级行政区,由于数据缺失或异常,剔除影像中存在异常值的样本数据,最终选择福建省的47个县级行政区用于本文研究,并随机选择36和11个县级行政区分别作为训练样本和验证样本．

以下为剔除异常样本的原因,具体分为3种情况:

1)县区所在的地形起伏明显,造成“珞珈一号”夜光遥感影像的几何误差较大,难以进行校正．具体县区为:罗源县、闽清县、连城县、松溪县、福安市、德化县、华安县和平和县．

2)县区被覆盖的热红外遥感数据云量较大．具体县区为:龙海市、漳浦县．

3)县区被不同时相的热红外遥感数据覆盖,难以反演所在地区的准确温度．具体县区为:武平县、浦城县、安溪县、永春县、长泰县、龙海市、平谭县和福清市．

1.3.1 夜间灯光遥感数据处理

首先对“珞珈一号”数据进行辐射定标,01星产品辐射亮度转换公式如下:

L=X3/2·10-10,

(1)

式中,L为定标后的辐射亮度,X表示“珞珈一号”影像单像元的原始值，即DN值．

然后对辐射定标后的数据进行几何校正等预处理．为了数据的有效性,将影像的辐射亮度值线性拉伸到0～64间．

1.3.2 Landsat8影像预处理及NDVI计算

首先结合元数据文件中的定标系数对表1中所选取的每景Landsat8 OLI影像数据进行辐射定标,获得辐射亮度影像;然后基于ENVI软件的FLASSH大气校正模块对每景影像进行大气校正;再将每景大气校正影像进行镶嵌,使得影像完全覆盖福建省,并用福建省行政矢量图裁切;最后基于式(2)计算NDVI(其量值记为INDVI),获得福建省NDVI影像,为地表温度反演提供输入数据．

(2)

式中,ρNIR为近红外波段,ρR为红光波段．

1.3.3 地表温度反演与温度归一化

地表温度计算公式为

L10=gD+b,

(3)

(4)

(5)

其中L10为Landsat8 TIRS传感器第10波段像元在传感器处的辐射值,D为像元灰度值,g为波段增益值,b为波段偏置值．对于TIRS传感器的第10波段而言,g为3.342 W/(m2·sr·μm),b为0.1 W/(m2·sr·μm),T为传感器处的亮温值,K1=774.89 W/(m2·sr·μm),K2=1 321.08 K,TLST为反演的地表温度,λ为热红外波段中心波长,ρ=1.438×10-2m·K,ε为地物比辐射率[17],其值根据Sobrino等[18]的模型通过NDVI进行估算得到．

有关研究表明高温区域主要分布在建筑用地密集、功能单一以及周边植被覆盖少的地区[19],因此地表温度与城市工业化呈正相关性．但由于温度本身的自然属性,即使在第二、第三产业不发达地区的实际温差依然存在,故本文采取以县为单位进行归一化处理,以此来提高温度与工业化的相关性,这在一定程度上也可以减小时相差异所带来的影响．归一化公式如下:

(6)

式中,TLST为所在像元的地表温度,TLST,min和TLST,max分别代表像元所在县的最低地表温度和最高地表温度,每个像元的温度值取在0～1之间,0代表此地增加值为最低．

1.3.4 模型建模及校正

本文采用逐步回归法构建模型:

Y1=β1X1+β2X2+…+βiXi,

(7)

式中,Y1为因变量,即相应区(县)级行政单元GDP的回归值(预测值),Xi(i=1,2,…,n)为回归模型自变量,即不同土地利用类型的灯光参数指标和温度指标,βi为Xi的偏回归系数．

由于直接采用预测模型模拟的第二、第三产业增加值误差较大,还需要根据第二、第三产业增加值的实际统计增加值对模拟增加值进行线性纠正,最后生成纠正后的第二、第三增加值密度图,使得误差只分布在该县内部[12],这样有利于进一步提高空间化的准确性．线性纠正公式为

PT=Pi(Pt/PS),

(8)

式中,PT为纠正后的单个像元的分产业增加值,Pi为像元i模拟的分产业增加值,Pt为该区县统计的实际分产业增加值,PS为该区县模拟的分产业增加值．

2 结果与分析

2.1 GDP 空间化建模

2.1.1 空间化模型

为了评价热红外遥感数据在第二、第三产业增加值空间化中的潜力,本文基于夜间灯光遥感数据、土地利用类型数据和热红外遥感数据(下文简称方法2),并以传统的结合夜间灯光数据以及土地利用类型数据(下文简称方法1)为对照分别进行空间化建模,利用SPSS软件进行逐步多元回归分析,获得如表2所示的建模结果．表2中,P2表示第二产业增加值,P3表示第三产业增加值．从表2中可以看出:使用方法1对福建省2018年第二、第三产业增加值建模的逐步回归模型决定系数R2偏小,分别为0.743、0.776,其主要原因是夜间灯光数据无法准确反映第二产业分布信息;方法2因引入热红外遥感数据,使得其对福建省第二、第三产业增加值的多元回归模型决定系数R2分别达到0.966和0.870,较方法1有了显著提高．

表2 方法2与方法1回归模型比较

1)第二产业增加值空间化模型

经过逐步回归分析,方法2第二产业增加值模型入选的自变量为城镇用地的归一化温度、工矿用地的归一化温度以及工矿用地的灯光强度．建立的模型为

P2=158.741T51+73.378T53+88.15L53,

(9)

式中,P2为模拟的第二产业增加值,T51为城镇用地的归一化温度,T53为工矿用地的归一化温度,L53为工矿用地的灯光强度．城镇用地、工矿用地同原材料加工和碳排放有一定的关联,而碳排放跟局部地区温度有显著关系,因此城镇用地和工矿用地的温度与灯光强度变量最终被入选进行回归建模．

2)第三产业增加值空间化模型

在方法2中的第三产业增加值模型入选的自变量为城镇用地的归一化温度、城镇用地的灯光强度以及农村居民点的灯光强度．构建的模型为

P3=178.352T51+36.114L51+42.997L52,

(10)

式中,P3为模拟的第三产业增加值,L51为城镇用地的灯光强度,L52为农村用地的灯光强度．由于第三产业增加值主导因素很多,除了城镇用地的灯光数据和温度数据可以反映第三产业增加外,据调查福建省农村金融、信息等现代农村服务业正在逐年稳步发展,农村居民点也存在第三产业增加值,故农村居民点的灯光数据也被入选作为自变量．

2.1.2 模型验证

基于选取的11个县级行政区作为验证样本对采用方法2建立的第二、第三产业模型分别进行验证,各县相对误差均在0～40%之间,平均相对误差大致在20%左右(表3)．

表3中部分地区相对误差较大,其可能原因如下:

1)热红外数据时相差异．上杭县、清流县、沙县的第二产业增加值的拟合值高于真实值,且相对误差较大,其原因可能是热红外数据时间差异造成的．尽管对温度数据进行0～1的归一化之后减小了不同时相带来的误差,使得县与县之间温度变化相对较小,但是因单景影像温度整体偏高,导致一个县内放大或缩小其热分布差异的情况存在可能,且难以避免．

2)第二产业的构成差异．若第二产业以轻工业为主,而其能源消耗及碳排放不显著,从而导致热红外遥感难以反映这一类型的第二产业．例如泉州市的第二产业以轻工业为主,该地区的温度分布较为均衡,基本无高值存在,故第二产业增加值拟合误差相对较大,且拟合值低于真实值．

3)复杂地形导致影像数据几何校正精度下降．由于福建省多山地丘陵,地形坡度因子使得影像几何校正产品仍存在不同程度的几何误差．受“珞珈一号”夜间灯光数据本身的局限性影响,辐射校正后部分山地地区的误差仍然存在,导致某些农村地区与其对应的灯光数据难以充分匹配,例如柘荣县、永安市、永泰县,从而造成了第三产业增加值拟合的误差相对较大,且拟合值低于真实值．

为了与方法1进行对比,基于选取的11个县区级行政区作为验证样本对采用方法1建立的第二、第三产业增加值模型分别进行验证,各县区相对误差较大(表4).

从表4可以看出,各县区第二、第三产业增加值模型相对误差较大,平均值分别为72.60%、60.10%．究其原因主要在于只用单一的夜间灯光数据不足以区分第二、第三产业的空间分布情况,所以前人大多将第二、第三产业增加值看作一个整体,与结合土地利用类型下的夜光遥感数据进行建模,从而忽略了单产业空间化建模,为单产业空间化分布研究造成了困难．还有一点造成方法1模型验证不理想的原因为福建省中西部山地多、地形起伏大的特殊地形,导致“珞珈一号”夜光影像数据几何校正精度下降,不能完全与实际土地利用进行匹配,这一原因也同样存在于方法2 的研究中,为本文研究中的不足之处．但总体来看,热红外遥感数据结合夜光遥感数据能够较好地分别与第二、第三产业增加值进行回归分析,建立的模型更加准确,与实际情况较为相符．

2.2 空间化结果分析

基于2.1中建立的2018年福建省第二、第三产业增加值模型,并经过校正后分别生成2018年福建省第二、第三产业增加值分布(图1a和图1b)．根据空间化结果发现福建省第二、第三产业增加值分布具有如下特点:

表3 方法2模型精度验证

表4 方法1模型精度验证

1)第二产业增加值高值区主要位于福建省的东部以及东南部,如泉州、福州和厦门等主要城市;第二产业低值区主要位于如南平大部分地区和宁德部分地区的福建省北部,该区域第二产业相对落后．

2)第三产业增加值分布特点与第二产业相似,其高值区集中分布在东部以及东南沿海地区,而西北部以及一些多山地区第三产业则欠发达．

3)受福建省多山地的地理分布格局的影响,增加值空间分布多山地区增加值较低,平原沿海地区经济分布高;内部局部小区域也有高值,如三明市辖区．此外,莆田市虽在福建省的东部,其第二产业和第三产业经济均较为落后．

2.3 空间化优势分析

厦门是福建副省级城市,为我国东南重要的中心城市、港口城市和优秀旅游城市,也是我国的经济特区、综合配套改革试验区、自由贸易试验区、两岸新产业和现代服务业合作示范区．因此,厦门市的第二、第三产业均较发达．为了进一步评价热红外遥感数据在GDP空间化中的潜力,下文以厦门市为例,对本文提取的新空间化方法与传统空间化方法进行对比．

2.3.1 第二产业增加值空间化对比

图2a为未考虑热红外数据的方法1的厦门市2018年第二产业增加值空间化结果,从中可以看出第二产业增加值空间分布较为零散,而且分级效果不明显．图2b为本文提出的方法2生成的第二产业增加值空间分布,通过与厦门市主要的工业园区分布图对比发现:方法2空间化结果与工业园区分布较为吻合,且增加值分级显著．

图1 2018年福建省第二、第三产业增加值空间分布Fig.1 Distribution map of added values of the secondary (a) and tertiary industry (b) of Fujian province in 2018

图2 2018年厦门市第二产业增加值不同空间化方法对比Fig.2 Comparison of added values of the secondary industry in Xiamen in 2018 between method 1 (a) and method 2 (b)

为了能从局部详细验证以上结论,进一步选取厦门市特新工业园区、厦门市方品工业园区以及厦门市鼓浪屿等典型局部区域进行两种空间化方法对比．

图3 厦门市典型区第二产业增加值空间化局部对比Fig.3 Comparison of added values of the secondary industry in typical districts of Xiamen between method 1 (a), method 2 (b),and actual distribution (c)

1)厦门市特新工业园区属于典型的第二产业分布区,方法2反演的第二产业增加值较高(图3b1),与实际情况相吻合(图3c1),而方法1空间化结果则与实际差异较大,甚至无明显的高值区(图3a1)．

2)厦门市方品工业园区也属于典型的第二产业分布区,方法2与方法1空间化结果分别见图3b2和图3a2,其差异与厦门市特新工业园区类似．

3)据调查,厦门市鼓浪屿旅游业发达,属于典型的第三产业分布区,而第二产业增加值较小,方法2反演的第二产业增加值较低(图3b3),符合实际情况,而方法1反演的第二产业增加值较高(图3a3),与实际情况差异较大．

4)方法1在海沧区第二产业增加值存在大范围高值区,根据高分辨率影像考察发现该地区主要为商业区以及居民居住用地,较强的灯光强度使得方法1反演结果与实际不相符,而方法2反演结果则更为合理．

因此,本文提出的方法2相对方法1的合理性在于第二产业把“工矿用地高于其他区域的地表温度的影响”纳入增加值预测,从而导致地表温度显著高于其他周边区域．因此,方法2通过结合热红外遥感数据考虑了地表温度信息可以显著提高第二产业增加值的预测与空间化效果．

2.3.2 第三产业增加值空间化对比

图4为基于方法2与方法1对厦门市2018年第三产业增值值空间化结果．方法1(图4a)和方法2(图4b)均能从不同程度反映厦门市第三产业经济的空间分布情况．厦门市第三产业主要集中分布于厦门岛的湖里区及思明区、集美区的杏林镇及其区政府所在地、同安区区政府所在地以及翔安区的马巷镇．但方法1的第三产业增加值分布较为破碎,且相比之下,方法2空间化结果更佳．

厦门市是福建省经济最发达的省辖市之一,尤其是厦门岛(图5c1),岛内的思明区是厦门市的政治经济中心,厦门市的主要旅游景点(如鼓浪屿、曾厝垵等)及商业街(思明区的中山路步行街、世贸商城、湖里区的SM城市广场)均在辖区内，旅游业非常发达,同时也带动第三产业迅速发展,因此该区域除了南部思明区的云顶岩和北部高崎国际机场外,岛内其他区域第三产业增加值较高．方法2空间化结果(图5b1)中增加值高值区域比方法1(图5a1)面积更大且更均匀,因此前者空间化结果相对更为合理．此外,厦门市集美区(图5c2)分布着商业区(包含购物中心、银行等服务业),为典型的第三产业分布区,方法2反演的第三产业增加值较高(图5b2),与实际情况相吻合(图5c2),而方法1(图5a2)空间化结果中高值区小且破碎．因此,本文提出的方法2对第三产业增加值空间化效果更为合理．

图4 2018年厦门市第三产业增加值空间分布对比Fig.4 Comparison of added values of the tertiary industry in Xiamen in 2018 between method 1 (a) and method 2 (b)

图5 厦门市典型区第三产增加值空间化局部对比Fig.5 Comparison of added values of the tertiary industry in typical districts of Xiamen between method 1 (a), method 2 (b),and Landsat8 image (c)

3 结论

针对GDP空间化研究中难以对第二、第三产业增加值分别进行空间化的问题,本文以福建省为研究区,在土地利用数据、夜间灯光遥感数据的基础上,通过增加热红外遥感数据考虑地表温度信息,以实现第二、第三产业增加值的空间化．研究结论如下:

1)与利用夜间灯光遥感数据结合土地利用数据进行GDP空间化相比,增加热红外遥感数据的本文方法构建的第二、第三产业增加值空间化模型(R2分别为0.966、0.870)均优于夜间灯光遥感数据+土地利用数据方法(R2分别为0.743、0.776)．本文方法的平均相对误差约20%．进一步以厦门市为研究对象,发现本文方法由于热红外遥感数据考虑了地表温度信息,可以大幅提高第二、第三产业增加值空间化效果,其空间分布格局与实际更为吻合．

2)空间分布上从总体上看,福建省第二产业增加值与第三产业增加值分布特点相似,高值区主要集中在东部及东部沿海地区．受福建省多山地的地理分布格局的影响,增加值空间分布呈多山地区产值增加值较低,而平原沿海地区较高的特点．

3)第二产业可分为重工业和轻工业．对于工业园区,其工业生产带来的碳排放等会影响区域上空的热力因子,与温度数据有较强的相关性;而对于以轻工业为主的地区,其正常生产与温度信息相关性不大,从而导致空间化误差可能会偏大．

由于Landsat热红外遥感数据自身的局限性,无法同时获取多幅无云的2018年冬季的福建省热红外遥感数据,尽管本文对反演的温度进行归一化处理以进行修正,但是仍无法解决单一县级单元内同时包含多幅不同时相的热红外遥感数据所造成的误差,后续可以运用其他的数据或算法来解决这一问题．同时,本文研究时未能考虑到地形因素对模型的影响,而本文研究区多山地,地形起伏大,地形坡度因子成为一个不可忽略的影响因子,致使部分处于山地丘陵地区的中西部地区难以与其对应的灯光数据进行准确的几何匹配,使得地形复杂地区反演结果误差偏大．在后续研究中还需进一步评估其他地形类型区域对本文方法的适用性．