CN110110025A - 基于特征向量空间滤值的区域人口密度模拟方法 - Google Patents

Abstract

一种基于特征向量空间滤值的区域人口密度模拟方法，包括获取区域矢量文件和统计数据，选择特征向量空间滤值法，以区域夜间灯光平均亮度作为自变量，选择辅助自变量；根据区域矢量文件对遥感夜间灯光影像进行处理，计算区域总亮度和平均亮度；建立邻接关系，得到空间邻接矩阵并进行中心化，计算矩阵特征值和特征向量；提取合适的特征向量作为夜光亮度的空间影响因子，添加到自变量中，求解回归系数，得到人口密度的特征向量空间滤值回归模型，根据模型实现区域人口密度模拟。本发明能够有效消除空间异质性和空间自相关性对人口密度分布的影响，采用自动化手段替代人工统计，节约人力物力，对于城市化智能监测、环境污染检测等应用具有重要意义。

Description

基于特征向量空间滤值的区域人口密度模拟方法

技术领域

本发明属于空间统计分析服务应用技术领域，以夜间灯光影像作为材料，特别涉及一种基于特征向量空间滤值的区域(城市级)人口密度模拟方法。

背景技术

夜间灯光指从太空观测夜间无云的地球时，人类聚居区和经济带发出的夺目的光芒，其来源主要是照明设施的普及。由于夜光遥感影像的独特魅力，谷歌地球已经将灯光影像作为影像层之一。科学工作者可以对这些影像进行数据挖掘从而发现社会或者自然规律。相比于普通遥感卫星影像，夜光遥感影像能更多地反映人类活动，这一特点使其在社会科学领域和自然科学领域都得到了广泛应用。

大量的统计研究表明夜间灯光与国民生产总值(GDP)或区域生产总值(GRP)存在较高的相关性。夜间灯光数据由于其覆盖范围广、时间跨度大、经济易获取等优点，成为了解经济发展动态、监测人类活动的有效数据源。夜光遥感技术被广泛地应用于经济参数估算、区域发展研究、城市化监测、光污染等诸多研究领域，能够客观地反映经济的变化趋势，便于在空间大尺度下进行时空动态分布研究，对分析中国中部地区经济空间格局具有重要意义。如美国国防气象卫星计划Defense Meteorological Satellite Program(DMSP)是美国国防部极轨卫星项目，DMSP上的线性扫描业务系统Operational Linescan System(OLS)最初是专门为云层监测设计的振荡扫描辐设计，它的两个波段可适用于白天和夜间进行观测。目前，DMSP/OLS夜间灯光影像主要用于城镇扩展研究，经济因子估算以及其它环境、灾害、渔业、能源等领域。通常来说，卫星传感器获取的主要是地表的太阳辐射反射信号，而DMSP/OLS传感器独辟蹊径，采集的是夜间灯光、火光等产生的辐射信号。DMSP/OLS传感器在夜间工作，能探测到城市灯光甚至小规模居民地、车流等发出的低强度灯光，并使之区别于黑暗的乡村背景。因此，以DMSP/OLS为代表的夜间灯光影像可作为人类活动的表征，成为了人类活动监测研究的良好数据源。除了DMSP-OLS数据之外，suomi national polarorbiting partnership(NPP)卫星的visible infrared imaging radiometer suite(VIIRS)传感器，作为DMSP卫星的继任者，有着740米的空间分辨率，星下分辨率更是高达400米，从2011年至今，其数据量已经极为丰富，每日影像均可免费下载。

夜间灯光亮度反演的必要性，首先是数据的获取性，DMSP-OLS以及NPP-VIIRS数据获取并不复杂，其中种类繁多，年合成，月合成，是否去云去火点等，。且DMSP-OLS数据的日合成数据需要付费下载，未对大众免费开放，这样使人们获取2013年之前的夜光数据上造成困难。其次是数据的完整性和连续性，由于DMSP-OLS数据在2013年已经不再更新，而2011年开始更新的NPP-VIIRS图像与前任有诸多不同，不但提升了分辨率，也去除了饱和效应。然而这在一定程度上给数据的连续性分析带来了困难，比如如果想研究改革开放四十年来中国夜光亮度的变化情况，则需要将时间序列分成两段，因为两段时间的研究载体是不同的。第三点则是数据的完整性以及数据质量问题，常用的DMSP-OLS图像分辨率较低，只有2700米，且存在过饱和效应，NPP-VIIRS数据中，每天的数据是未去火点的，即便是月合成数据，异常值也非常多。第四则是数据处理的复杂性，NPP-VIIRS数据的处理过程繁琐，需要进行异常值处理，空缺值处理等一系列步骤。

现有发明中有很多利用夜间灯光研究例如ELVIDGE C D等利用DMSP-OLS对每周21个国家的夜间发光面积和GDP进行回归分析，发现回归的决定系数达到0.97.此后，类似的研究在欧盟，中国，美国等国分别展开，分析发现夜光总量与GDP/GRP的R²在这些区域均可达到0.8到0.9。布朗大学的研究显示可以通过夜间灯光亮度来修正国民生产总值的增长率。CHEN利用1992年至2008年的夜光遥感影像数据提出了了全球GDP格网产品。考虑到夜间灯光在不同程度上与人口，经济等均存在较高的相关性，因此夜间灯光为经济指标的空间化提供了新的途径。例如人口的估算，LO C P利用中国1997DMSP/OLS夜光影像的点面积，夜光总量，平均发光亮度，夜光比例等多个指标与人口，非农业人口在县级和城市单元上进行回归分析发现夜光数据能在县级尺度上很好地模拟非农业人口。利用夜光遥感数据和一些辅助数据如土地利用，遥感植被数据可得到更为精细地人口网格数据。夜光影像能够反映照明设施地密度和使用程度，因此能够为电力消费空间化提供依据。

空间滤值方法最早是由Getis和Griffith提出，该方法的核心思想是，把模型中的变量分解成空间影响和非空间影响两部分，即将变量的空间影响部分提取出来并将其滤去，剩下无空间相关性地可以利用常用的回归方法进行分析。Getis提出的空间滤值方法利用局部Gi统计将自变量通过公式进行变化，实现残差中空间影响部分的滤去，这种方法能滤去空间自相关影响，但是对变量的要求是必须满足为正值。Griffith提出的特征函数空间滤值方法通过选取特征向量加入到自变量中构建滤值算子来代替模型残差中的自相关部分，使得剩余的残差部分只受随机误差影响，从而消除空间自相关的影响。滤值算子相当于残差的自相关部分，需要包含地理单元间的空间关系。Patuelli利用空间滤值方法研究德国失业现象，发现空间滤值的加入使得回归模型对于失业现象的预测准确率提高，从实证研究的角度验证了空间滤值方法的有效性。Chun提出了一种更快更有效的生成特征向量子集的方法，使得空间滤值方法的效率大大提高。但是，目前尚未有研究提出用于区域人口密度模拟。

由上文所述，夜间灯光数据有着良好的可获得性，而且与部分人文经济参量的相关性大，因此具有通过夜间灯光而对区域人文经济参量进行估计的潜力。传统利用夜间灯光反演人口密度的研究中，大多未考虑空间相关性，造成残差中的空间相关性较高。

发明内容

为了解决人口统计复杂，流程长，周期长，无法快速获得等缺点，以及偏远区域统计困难，统计数据不完全等问题，本发明以夜间灯光为材料，利用特征函数空间滤值法对区域人口密度进行模拟。

本发明采用的技术方案是一种基于特征向量空间滤值的区域人口密度模拟方法，模拟过程包括以下步骤：

步骤1，数据获取和模型选择，包括获取区域矢量文件和统计数据，选择特征向量空间滤值法，以区域夜间灯光平均亮度作为自变量，按照以下原则选择辅助自变量，

一是变量与人口密度之间存在显著相关性；

二是加入所有自变量之后模型不存在严重的共线性问题；

步骤2，下载遥感夜间灯光影像，根据步骤1获得的区域矢量文件对遥感夜间灯光影像进行处理，计算区域总亮度和平均亮度；

步骤3，针对步骤1获得的区域矢量文件建立邻接关系，得到相应的空间邻接矩阵W₀，并对空间邻接矩阵W₀进行中心化得到矩阵W₁；

步骤4，计算矩阵W₁的特征值和特征向量；

步骤5，提取合适的特征向量作为夜光亮度的空间影响因子，

步骤6，将提取的所有特征向量添加到自变量中，使用最小二乘法求解回归系数，得到人口密度的特征向量空间滤值回归模型，

y＝β₀+β₁MEAN+β₁X_aux+β_kE_k+ε

其中，y代表地级市的人口密度，X _aux表示除区域平均亮度MEAN以外其他自变量的集合，E_k表示最终选择的特征向量集合，β₀、β₁和β_k为系数，ε表示残差，ε服从正态分布；

步骤7，根据人口密度的特征向量空间滤值回归模型实现区域人口密度模拟。

而且，步骤2中，根据步骤1获得地级市矢量文件对遥感夜间灯光图像进行处理，包括裁剪、投影转换、几何校正和异常值矫正处理。

而且，步骤5中，提取合适的特征向量作为夜光亮度的空间影响因子，包括首先筛选出所对应莫兰指数除以最大莫兰指数大于等于0.25的特征向量，然后使用前向选择法从中挑选合适的特征向量作为空间影响因子。

而且，对步骤6建立的人口密度的特征向量空间滤值回归模型进行精度评价。

而且，包括计算所得人口密度的特征向量空间滤值回归模型的拟合优度、调整后拟合优度、均方根误差、残差莫兰指数、和均方根误差作为评价指标，以验证所提出的模型的精度，与之对比的模型有普通线性回归，空间滞后模型，空间误差模型。

而且，选择辅助自变量的实现方式为，

1)设有辅助自变量备选集{X₁，X₂，……，X_n}，计算每个备选变量与夜间灯光亮度的Pearson指数及显著性检验，剔除未通过显著性检验的变量；

2)对剩余自变量{X₁，X₂，……，X_m}进行共线性诊断，剔除方差膨胀因子VIF值大于10的自变量，得到最终所选的辅助自变量{X₁，X₂，……，X_k}。

本发明所提供的基于特征向量空间滤值法的人口密度回归能够有效消除空间异质性和空间自相关性对人口密度分布的影响，提高估算模型的精度，建模过程和模型结构简单，可以有效模拟出地级市级别城市人口密度情况，采用自动化手段替代人工统计，可以节约大量人力物力，在一定程度上可以定量分析对人口分布的影响因素，对于后续的城市规划，城市化智能监测、环境污染检测等应用具有重要意义。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例的夜光影像处理流程图。

图3为本发明实施例的自变量筛选流程图。

具体实施方法

为了便于本领域普通技术人员理解以及实施本发明，下面结合附图及实施例对本发明作进一步阐述，应当理解，此处所述的实施例仅用于说明和解释本发明，并不用于限定本发明。

本发明要解决的问题是：人口统计流程繁复，在进行大范围统计时虽能获得精确数据，但费时费力，周期长，人力高远，人口稀疏地区无法准确统计，且受人口流动影响大，难以获得实时人口分布信息。传统拟合人口密度的方法如线性回归方法，无法有效消除模型残差的空间相关性，故本发明实施例利用特征向量空间滤值法，且利用夜间灯光影像反映实时信息，为人口分布提供参考。

空间滤值法虽然能有效降低拟合结果中残差的空间相关性，以提高拟合精度，但是也存在着计算量大，特征向量计算，筛选繁琐等问题。因此本发明首先对研究区域进行规划筛选，控制地级市数量，从而避免计算量过大。利用特征向量空间滤值法的优势有三点：1.由于人口统计流程复杂繁琐，受地区，时间影响大，例如偏远地区统计困难，且部分时间人员迁徙流动大，而基于夜间灯光影像可以有效反映实时人口密度情况，为人口数据提供参考。在于当确定研究区域集合的空间邻接关系后，根据该邻接关系确定的邻接矩阵可以重复利用，因为对于特定的矩阵而言，其特征向量是确定的，因此计算过程类似，可重复利用。3.相对于空间误差模型，空间滞后模型等考虑到空间相关性的模型，空间滤值法的精度更高。

基于特征向量空间滤值法由GDP，区域平均夜间灯光亮度，电力消耗量等对人口密度建模，实施例提供的一种基于特征向量空间滤值的区域人口密度模拟方法流程参见附图1，包括以下步骤：

步骤1：数据获取和模型选择，包括获取区域矢量文件和统计数据，选择特征向量空间滤值法。

具体实施时，数据获取可采用的方式为：收集整理城市统计年鉴，空值数据查询补全，整理配准地级市矢量文件。

实施例中，步骤1实现如下，

步骤1.1：从权威数据源获得必要数据，模型选择特征向量空间滤值法，从中国国家统计局(http://www.stats.gov.cn/)获取省市，主要地级市统计年鉴数据。

步骤1.2：统计数据预处理，由于地级市个数从2010年的283个逐渐增长到2017年的294个，故地区数据配准还需要参考城市统计年鉴的县级数据，将升为地级市的区域单独统计配准。所选地级市整体必须相互邻接，遇到飞地合并到当地，遇到岛的情况如舟山市则建立与之最近的城市邻接关系，较小面积的海岛则删除。筛选数据时遇到缺失值则查询该地级市所在官网重新查找统计年鉴，以补全数据。区域矢量文件需要与年鉴数据相互配准。以区域夜间灯光平均亮度作为自变量，从年鉴中选择GDP，耗电量作为辅助自变量。辅助自变量的选择应符合两个原则：一是变量与人口密度之间存在显著相关性，以皮尔逊指数(Pearson指数)衡量，公式如下：

其中，和分别表示X、Y相应的平均值，其中X表示因变量Y对应的自变量(如某个城市的平均亮度)，i用于表示第i个地理对象(如第5个地级市)，n是所选城市总数，即X_i表示第i个地理对象的自变量，Y_i表示第i个地理对象的因变量，PC表示Pearson指数。若PC值不等于0且通过了显著性检验(PC小于等于0.1)，则该变量与灯光亮度相关，可以保留下来。二是加入所有自变量之后模型不存在严重的共线性问题，需要构建初步的线性回归模型进行共线性诊断，剔除方差膨胀因子(VIF)大于10的变量，该步骤可在R语言中使用VIF函数进行判断([1]张宇山.多元线性回归分析的实例研究[J].科技信息,2009(09):54-56.)。Pearson系数计算与共线性诊断都可在R语言中进行。参见图3，实施例选择辅助自变量的实现方式为以下流程，

1)设辅助自变量备选集{X₁，X₂，……，X_n}，计算每个备选变量与夜间灯光亮度的Pearson指数计算及显著性检验，剔除未通过显著性检验的变量(p值大于0.1)；

2)对剩余自变量{X₁，X₂，……，X_m}进行共线性诊断，剔除方差膨胀因子(VIF)值大于10的自变量，得到最终所选的辅助自变量{X₁，X₂，……，X_k}。

步骤2：下载遥感夜间灯光影像，根据步骤1获得地级市矢量文件对遥感夜间灯光图像进行裁剪，投影转换，几何校正，异常值矫正等处理，计算区域总亮度和平均亮度。

实施例中，夜光影像的下载和处理参见附图2。夜间灯光影像选取NPP/VIIRS年合成影像(https://governmentshutdown.noaa.gov/)。亦可选择年底的月合成影像。

步骤2.1：投影转换和裁剪。影像自身采用WGS84坐标系，需要转换成CGCS2000坐标系，若矢量文件本身为WGS84坐标系则不用转换。进行转换后，根据矢量文件进行裁剪，选出研究区域，裁剪滞后栅格整体亮度值可能变深，是由于原栅格和结果栅格的拉伸方式和统计值不同造成，如果亮度值发生变化，将原栅格统计数据另存为XML文件，再加载到结果栅格中即可。

步骤2.2：亮度值异常值矫正处理，由于所选年合成影像数据进行了亮度校准，故最低亮度值均为0，但异常高亮度较多，实施例中将北京，上海，广州，深圳等中国经济最发达的地区中的最高亮度值作为亮度阈值。最后将空置数据NO DATA统一赋值为0。以上步骤均可在ArcGIS Desktop中的栅格计算器中进行。

步骤3：针对步骤1获得的区域矢量文件建立邻接关系，得到相应的空间邻接矩阵W₀，并对空间邻接矩阵W₀进行中心化得到矩阵W₁。

实施例中，步骤3构建空间邻接关系，根据矢量文件的空间邻接关系构建邻接矩阵，包括采用rook连接，由多边形的邻接关系构建二元邻接矩阵W₀，即多边形i和j相邻，则元素W₀(i，j)等于1，否则等于0，所构建邻接矩阵记为W₀。

步骤4：计算矩阵W₁的特征值和特征向量，将所有特征向量记录为{E₁,E₂,……，E_n}。E₁表示第一个特征向量，E₂表示第二个特征向量，以此类推到第n个特征向量E_n。

实施例中，中心化邻接矩阵W₀得到矩阵W₁,计算矩阵W₁的特征值和特征向量，矩阵中心化的公式如下：

其中I为n维单位矩阵，11^T是一个n×n的矩阵，矩阵内所有元素都等于1，n是城市的数量。再使用数学分解的方法，求解W₁的特征值和特征向量E＝{E₁，E₂，E₃，……，E_n}，其中W₀为步骤3中的空间邻接矩阵。

完成特征值计算后，进行特征向量的预选择，要求其特征值不为0，同时一般遵循规律为：

其中，λ_i表示第i个特征向量所对应的特征值，λ_max表示最大的特征值。该过程可使用Matlab等软件中自带函数Eig()进行计算，也可用R语言中的spmoran包实现。

步骤5：提取合适的特征向量作为夜光亮度的空间影响因子，

实施例对步骤4得到的特征向量首先筛选出所对应莫兰指数除以最大莫兰指数大于等于0.25的特征向量，然后使用stepwise的前向选择法从中挑选合适的特征向量作为空间影响因子，具体步骤如下：

步骤5.1：对X、Y中心化得X_cent和Y_cent，X表示自变量，Y表示因变量，对其进行中心化后得到相应变量X_cent和Y_cent，求解回归残差e和e的莫兰指数Moran’s I，公式如下：

X_cent＝(I-P)X,Y_cent＝(I-P)Y

其中，I为n×n的单位矩阵；P为一个n×n的矩阵，且所有元素的值都为1/n，n是地级市的数量；

步骤5.2：对e的Moran’s I进行显著性检验。显著性检验的方法具体如下：

a.计算残差e和标准化的莫兰指数Moran’s I₀；

b.随机排列残差得e_rnd,计算随机排列后的残差的莫兰指数Moran’sI_rnd；

c.步骤b重复进行999次,计算其PC值p2，p2＝(num+1)/(999+1)。

其中，num为999次随机排列中，Moran’sI_rnd大于Moran’s I₀的个数。若p2值小于阈值(一般取0.05或0.01)，结果显著，说明回归残差存在空间自相关性，不符合线性模型的假设，执行步骤5.3。否则，结果不显著，则执行步骤6。

步骤5.3：X的原始值包括区域平均亮度，GDP等，循环遍历所有特征向量，每次从E中选择一个特征向量E_i加入到自变量X中，即

X＝X+E_i

逐个计算回归残差Moran’s I，当所有特征向量遍历一遍之后，所有残差Moran’sI记位i₁,i₂,……，i_n，组成向量I＝{i₁,i₂,……，i_n}。选取最小残差Moran’s I所对应的特征向量，将其作为一个新的自变量保留下来，并从原特征向量组E中剔除，然后重复步骤5.2，进入下一次迭代。

步骤6：将提取的所有特征向量添加到自变量中，使用最小二乘法求解回归系数，得到人口密度的特征向量空间滤值回归模型。

y＝β₀+β₁MEAN+β₁X_aux+β_kE_k+ε

其中：y代表地级市的人口密度，X _aux表示除了MEAN(区域平均亮度)，其他自变量的集合(如GDP等)，E_k表示步骤5最终选择的特征向量的集合，通过极大似然法或限制性的极大似然法可以对上述式子进行求解，得到β₀，β₁，β_k各系数的取值，ε表示残差，上述值即为拟合结果，ε服从正态分布，即消除了残差的空间相关性，提高了拟合的精确程度。

步骤7：根据人口密度的特征向量空间滤值回归模型实现区域人口密度模拟。模拟结果中，拟合结果决定系数R₂可达到0.9以上，误差在百分之十以内。

具体实施时，以上流程可采用计算机软件技术实现自动运行流程。实现该方法的装置也应当在本发明的保护范围内。

为证明本发明方案的有效性，保证模拟效果，可进行模型精度评价与模型对比选择。

步骤1)：计算所得模型的拟合优度R²、调整后拟合优度R²(Adj.R²)、均方根误差(RMSE)、平均绝对误差(MAE)以及残差的莫兰指数Moran’s I、均方根误差MSE等作为评价指标，以验证所提出的基于特征向量空间滤值法的精确性。与之对比的模型有普通线性回归，空间滞后模型，空间误差模型。

为便于实施参考起见，提供各指数计算方式如下：

其中n是地级市个数，y_i是第i个城市的城市人口密度的真实值，表示第i个城市人口密度的拟合值。

其中m是自变量的个数；R²和Adj.R²的取值范围是[0,1],值越大说明模型精度越高。

式中参数含义同上，RMSE和MAE越小说明模型精度越高其中是由模型得到的夜间灯光亮度的残差，是平均值，是城市i和j之间的反距离空间权重。的取值范围是[-1,1]，值越接近于0，残差空间自相关性越弱，模型越可靠。

步骤2)：根据精度评价选择最优模型。进行了步骤1)的模型精度评估之后，进行比较拟合优度(R²)、调整后拟合优度(Adj.R²)、均方根误差(RMSE)、平均绝对误差(MAE)、残差Moran’s I以及均方根误差RMSE等评价指标。可以获得其中一个最优的模型作为最终模型。

应当理解的是，上述针对本发明中较佳实施例的表述较为详细，但不能因此认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在本发明权利要求保护范围内，可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围以所附权利要求书为准。

Claims (6)

1.一种基于特征向量空间滤值的区域人口密度模拟方法，其特征在于，模拟过程包括以下步骤：

步骤1，数据获取和模型选择，包括获取区域矢量文件和统计数据，选择特征向量空间滤值法，以区域夜间灯光平均亮度作为自变量，按照以下原则选择辅助自变量，

一是变量与人口密度之间存在显著相关性；

二是加入所有自变量之后模型不存在严重的共线性问题；

步骤2，下载遥感夜间灯光影像，根据步骤1获得的区域矢量文件对遥感夜间灯光影像进行处理，计算区域总亮度和平均亮度；

步骤3，针对步骤1获得的区域矢量文件建立邻接关系，得到相应的空间邻接矩阵W₀，并对空间邻接矩阵W₀进行中心化得到矩阵W₁；

步骤4，计算矩阵W₁的特征值和特征向量；

步骤5，提取合适的特征向量作为夜光亮度的空间影响因子，

步骤6，将提取的所有特征向量添加到自变量中，使用最小二乘法求解回归系数，得到人口密度的特征向量空间滤值回归模型，

y＝β₀+β₁MEAN+β₁X_aux+β_kE_k+ε

其中，y代表地级市的人口密度，X_aux表示除区域平均亮度MEAN以外其他自变量的集合，E_k表示最终选择的特征向量集合，β₀、β₁和β_k为系数，ε表示残差，ε服从正态分布；

步骤7，根据人口密度的特征向量空间滤值回归模型实现区域人口密度模拟。

2.根据权利要求1所述基于特征向量空间滤值的区域人口密度模拟方法，其特征在于：步骤2中，根据步骤1获得地级市矢量文件对遥感夜间灯光图像进行处理，包括裁剪、投影转换、几何校正和异常值矫正处理。

3.根据权利要求1所述基于特征向量空间滤值的区域人口密度模拟方法，其特征在于：步骤5中，提取合适的特征向量作为夜光亮度的空间影响因子，包括首先筛选出所对应莫兰指数除以最大莫兰指数大于等于0.25的特征向量，然后使用前向选择法从中挑选合适的特征向量作为空间影响因子。

4.根据权利要求1所述基于特征向量空间滤值的区域人口密度模拟方法，其特征在于：对步骤6建立的人口密度的特征向量空间滤值回归模型进行精度评价。

5.根据权利要求4所述基于特征向量空间滤值的区域人口密度模拟方法，其特征在于：包括计算所得人口密度的特征向量空间滤值回归模型的拟合优度、调整后拟合优度、均方根误差、残差莫兰指数、和均方根误差作为评价指标，以验证所提出的模型的精度，与之对比的模型有普通线性回归，空间滞后模型，空间误差模型。

6.根据权利要求1或2或3或4或5所述基于特征向量空间滤值的区域人口密度模拟方法，其特征在于：选择辅助自变量的实现方式为，

1)设有辅助自变量备选集{X₁，X₂，……，X_n}，计算每个备选变量与夜间灯光亮度的Pearson指数及显著性检验，剔除未通过显著性检验的变量；

2)对剩余自变量{X₁，X₂，……，X_m}进行共线性诊断，剔除方差膨胀因子VIF值大于10的自变量，得到最终所选的辅助自变量{X₁，X₂，……，X_k}。