论文精读：喀斯特山地流域耕地流转的时空演变与地形梯度效应——以贵州南北盘江流域为例

喀斯特山地流域耕地流转的时空演变与地形梯度效应——以贵州南北盘江流域为例一、研究背景喀斯特山区的耕地困境贵州南北盘江流域位于珠江上游是典型的喀斯特山地流域。这里山高坡陡碳酸盐岩厚度占沉积层85%以上坡度大于25°的面积约占45%。独特的地质背景造就了“地表破碎、地下河发育、土层浅薄”的生态环境也使耕地资源变得尤为珍贵。然而2000—2020年间随着城镇化加速、生态退耕政策实施以及农业结构调整流域内的耕地面积经历了“先增后减再增”的剧烈转型总体呈现减少趋势。耕地不仅数量减少还在空间上表现出明显的地形梯度效应低海拔、缓坡的优质耕地被建设用地占用高海拔、陡坡的坡耕地则逐渐撂荒或退耕还林还草。理解这种耕地流转的时空规律及其与地形的关系对于喀斯特地区的耕地保护和生态修复具有重要意义。本文基于2000—2020年的土地利用数据和DEM数据综合运用核密度分析、地形位指数、转移矩阵和弦图可视化等方法系统揭示了贵州南北盘江流域耕地流转的时空特征及地形梯度效应并附上完整的Python实现代码。二、研究区概况与数据准备2.1 研究区简介南北盘江流域贵州部分位于贵州省西南部地势自西北向东南递减海拔233~2856 m相对高差大。气候属中亚热带季风湿润区年均温16.2℃年降水1200~1650 mm。土地利用类型包括耕地、林地、草地、水域、建设用地和未利用地其中耕地主要分布在河谷盆地、山间坝子和缓坡地带。2.2 数据来源土地利用数据中国科学院资源环境科学数据中心http://www.resdc.cn2000—2020年共5期2000、2005、2010、2015、2020空间分辨率30 m。按《土地利用现状分类》重分类为6类。DEM数据地理空间数据云http://www.gsccloud.cnASTER GDEM v3.030 m分辨率。2.3 数据预处理在ArcGIS或Python中完成投影转换统一为WGS_1984_UTM_zone_48N土地利用栅格重采样保持30 m从DEM提取坡度、地形起伏度、地形位三、研究方法与代码实现3.1 核密度分析KDE核密度估计用于刻画耕地空间分布的聚集程度。计算公式[f(x) \frac{1}{nh} \sum_{i1}^{n} K\left(\frac{x - x_i}{h}\right)]我们使用scikit-learn或scipy实现。importnumpyasnpimportrasteriofromscipy.spatialimportcKDTreefromscipy.statsimportgaussian_kdeimportmatplotlib.pyplotaspltdefkernel_density_from_raster(raster_path,band1,sample_rate0.1): 从土地利用栅格中提取耕地类型代码1的坐标计算核密度 withrasterio.open(raster_path)assrc:datasrc.read(band)transformsrc.transform rows,colsnp.where(data1)# 假设耕地编码为1xs,ysrasterio.transform.xy(transform,rows,cols)coordsnp.array(list(zip(xs,ys)))# 如果点太多随机采样以加速计算iflen(coords)50000:idxnp.random.choice(len(coords),int(len(coords)*sample_rate),replaceFalse)coordscoords[idx]# 高斯核密度估计kdegaussian_kde(coords.T)# 生成网格xmin,ymin,xmax,ymaxcoords[:,0].min(),coords[:,1].min(),coords[:,0].max(),coords[:,1].max()xi,yinp.mgrid[xmin:xmax:200j,ymin:ymax:200j]zikde(np.vstack([xi.flatten(),yi.flatten()]))zizi.reshape(xi.shape)returnxi,yi,zi# 示例计算2000年耕地核密度并绘图xi,yi,zikernel_density_from_raster(landuse_2000.tif)plt.figure(figsize(10,8))plt.contourf(xi,yi,zi,levels20,cmapYlOrRd)plt.colorbar(labelKernel Density)plt.title(Cultivated Land Kernel Density (2000))plt.show()3.2 地形因子提取1坡度、地形起伏度importrasteriofromrasterio.featuresimportslopeimportnumpyasnpdefcalculate_slope_and_relief(dem_path,window_size3):withrasterio.open(dem_path)assrc:demsrc.read(1).astype(np.float32)transformsrc.transform cell_sizetransform[0]# 坡度度slope_degslope(dem,cell_sizecell_size)# 地形起伏度窗口内最大高程-最小高程fromscipy.ndimageimportmaximum_filter,minimum_filter reliefmaximum_filter(dem,sizewindow_size)-minimum_filter(dem,sizewindow_size)returnslope_deg,relief2地形位指数Terrain Position Index, TPI论文中公式为[T \ln\left[\left(\frac{E}{E_0}1\right)\times\left(\frac{S}{S_0}1\right)\right]]其中 (E, S) 为某点的高程和坡度(E_0, S_0) 为区域平均值。defterrain_position_index(dem,slope_deg):计算地形位指数E0np.mean(dem)S0np.mean(slope_deg)Tnp.log((dem/E01)*(slope_deg/S01))returnT3.3 耕地利用转移矩阵转移矩阵表示不同时期各地类之间的转换面积。使用numpy和pandas实现。importpandasaspdimportnumpyasnpdeftransfer_matrix(raster_before,raster_after,class_names):计算两个时期土地利用转移矩阵像元计数或面积withrasterio.open(raster_before)assrc1,rasterio.open(raster_after)assrc2:beforesrc1.read(1)aftersrc2.read(1)# 确保形状一致ifbefore.shape!after.shape:afterafter[:before.shape[0],:before.shape[1]]nlen(class_names)transnp.zeros((n,n),dtypenp.float64)foriinrange(n):forjinrange(n):mask(beforei1)(afterj1)trans[i,j]np.sum(mask)*(30*30)/1e6# 转换为km²dfpd.DataFrame(trans,indexclass_names,columnsclass_names)returndf# 示例2000→2005年转移矩阵classes[耕地,林地,草地,水域,建设用地,未利用地]trans_2000_2005transfer_matrix(landuse_2000.tif,landuse_2005.tif,classes)print(trans_2000_2005)3.4 弦图可视化Chord Diagram用plotly或holoviews绘制弦图直观展示地类间的流转。importplotly.graph_objectsasgodefchord_diagram(transfer_df,titleLand Use Transition):sources[]targets[]values[]labelstransfer_df.index.tolist()fori,srcinenumerate(labels):forj,tgtinenumerate(labels):ifi!jandtransfer_df.iloc[i,j]0:sources.append(i)targets.append(j)values.append(transfer_df.iloc[i,j])figgo.Figure(data[go.Sankey(nodedict(pad15,thickness20,linedict(colorblack,width0.5),labellabels),linkdict(sourcesources,targettargets,valuevalues))])fig.update_layout(titletitle,font_size12)fig.show()# 绘图chord_diagram(trans_2000_2005,2000-2005 Cultivated Land Transition)3.5 地形梯度分级与统计按照论文分级标准海拔233~900, 900~1200, 1200~1500, 1500~2000, 2000 m坡度0~2°, 2~6°, 6~15°, 15~25°, 25°地形起伏度0~30, 30~70, 70~200, 200~500 m地形位自然间断法分为15级再合并为5个等级defzonal_statistics(raster,zone_raster,zone_classes):按分区统计各类面积importrasteriofromrasterio.featuresimportgeometry_maskimportnumpyasnp# 简化实现利用numpy的bincountwithrasterio.open(raster)assrc1,rasterio.open(zone_raster)assrc2:landsrc1.read(1)zonesrc2.read(1)valid(land0)(zone0)land_valsland[valid]zone_valszone[valid]# 组合编码combinedland_vals*1000zone_vals countsnp.bincount(combined)# 解析...# 返回DataFramepass由于实际代码较冗长这里给出核心思路将地形因子分级后的栅格作为分区统计每个分区内耕地的像元数乘以像元面积得到面积。四、主要结果与讨论4.1 耕地时空分布演变特征1核密度结果2000—2020年间南北盘江流域耕地空间分布呈现“中部多、西部较多、东南部少”的格局。高密度区逐渐向流域中部集中低密度区面积由2000年的880 km²先增至2010年的1040 km²后降至2020年的760 km²反映了退耕还林和城镇化双重影响下耕地斑块的重新集聚。2不同地形条件耕地面积变化海拔耕地主要分布在1200~1500 mⅢ级2000—2020年该范围耕地面积基本稳定而1500~2000 m范围耕地显著减少。坡度6°15°Ⅲ级分布最广其次是15°25°Ⅳ级。所有坡度带耕地均减少陡坡25°减少幅度相对较小。地形起伏度0~30 mⅠ级集中了约70%的耕地但该区域耕地退化最严重。地形位中低地形位1~3级耕地面积占比最高且减少量占总减少量的43.55%表明低地形位是耕地流失的主战场。4.2 耕地流转的弦图分析通过弦图论文图4–7发现主要流转路径耕地 ↔ 林地 ↔ 草地三者之间频繁转换体现了退耕还林还草政策的持续影响。海拔梯度三级海拔1200~1500 m上耕地转出最多主要流向林地和草地中低海拔区域耕地向建设用地流转趋势明显。坡度梯度三级、四级坡度6°~25°耕地流转最活跃陡坡耕地主要转为林草地。地形位梯度三级地形位中等上耕地动态最剧烈表明该区域是人类活动与生态修复交锋的焦点。4.3 地形梯度效应与驱动机制论文提出了一个综合驱动模型图8归纳为三条路径低地形位缓坡、低海拔产业结构调整 → 耕地向高效经济作物或建设用地流转。驱动因素市场需求、城镇化、交通便利。中地形位中等坡度、中等海拔政策推动 劳动力迁移 → 退耕还林、生态修复。驱动因素退耕还林政策、农村人口外流。高地形位陡坡、高海拔资源禀赋限制 → 撂荒 → 自然恢复为林草地。驱动因素耕作成本高、生态脆弱、人口迁出。量化数据显示低地形位区域的PD斑块密度下降速率比高地形位区域高12个百分点LSI下降速率高23个百分点说明人类活动强烈的低地形位对尺度变化更敏感。五、结论与政策建议5.1 主要结论面积与格局变化2000—2020年南北盘江流域耕地总面积减少253.47 km²空间上向流域中部集聚高密度耕地比例提高。地形分布特征耕地主要集中在中低海拔1200~1500 m、缓坡6°15°、低起伏030 m及中低地形位区域。流转规律主要流转方向为“耕→草→林→耕”体现政策驱动的生态转换。建设用地侵占耕地主要发生在低地形位。地形梯度效应低地形位区域耕地变化最剧烈陡坡高海拔区域以撂荒和自然恢复为主。5.2 政策建议低地形位区域鼓励经果林种植、高标准农田建设同时严控建设用地占用优质耕地。中地形位区域巩固退耕还林成果实施生态补偿推广混农林业。高地形位区域禁止新开垦坡耕地加强生态修复监测发展生态旅游或碳汇项目。六、完整代码整合一键运行版以下是一个整合脚本执行从数据读取到结果可视化的全部流程需自行准备tif文件。importosimportnumpyasnpimportpandasaspdimportrasteriofromscipy.statsimportgaussian_kdeimportmatplotlib.pyplotaspltimportplotly.graph_objectsasgo# 1. 核密度 defplot_kde(raster_path,out_png):withrasterio.open(raster_path)assrc:datasrc.read(1)transformsrc.transform rows,colsnp.where(data1)xs,ysrasterio.transform.xy(transform,rows,cols)coordsnp.array(list(zip(xs,ys)))iflen(coords)50000:idxnp.random.choice(len(coords),50000,replaceFalse)coordscoords[idx]kdegaussian_kde(coords.T)xmin,ymin,xmax,ymaxcoords[:,0].min(),coords[:,1].min(),coords[:,0].max(),coords[:,1].max()xi,yinp.mgrid[xmin:xmax:200j,ymin:ymax:200j]zikde(np.vstack([xi.flatten(),yi.flatten()])).reshape(xi.shape)plt.figure(figsize(10,8))plt.contourf(xi,yi,zi,levels20,cmapYlOrRd)plt.colorbar(labelDensity)plt.title(os.path.basename(raster_path))plt.savefig(out_png,dpi300)plt.close()# 2. 地形位计算 defcalc_topographic_position(dem_path,slope_deg):withrasterio.open(dem_path)assrc:demsrc.read(1).astype(np.float32)E0,S0np.mean(dem),np.mean(slope_deg)Tnp.log((dem/E01)*(slope_deg/S01))returnT# 3. 转移矩阵 deftrans_matrix(before_path,after_path,class_names,out_csv):withrasterio.open(before_path)assrc1,rasterio.open(after_path)assrc2:beforesrc1.read(1)aftersrc2.read(1)ifbefore.shape!after.shape:afterafter[:before.shape[0],:before.shape[1]]nlen(class_names)transnp.zeros((n,n))foriinrange(n):forjinrange(n):trans[i,j]np.sum((beforei1)(afterj1))*900/1e6dfpd.DataFrame(trans,indexclass_names,columnsclass_names)df.to_csv(out_csv)returndf# 4. 弦图 defdraw_chord(trans_df,title,out_html):sources,targets,values[],[],[]labelstrans_df.index.tolist()fori,srcinenumerate(labels):forj,tgtinenumerate(labels):ifi!jandtrans_df.iloc[i,j]0:sources.append(i);targets.append(j);values.append(trans_df.iloc[i,j])figgo.Figure(data[go.Sankey(nodedict(pad15,thickness20,labellabels),linkdict(sourcesources,targettargets,valuevalues))])fig.update_layout(titletitle)fig.write_html(out_html)# 主程序 if__name____main__:# 假设文件列表years[2000,2005,2010,2015,2020]class_names[耕地,林地,草地,水域,建设用地,未利用地]# 计算每个年份的核密度图foryinyears:plot_kde(flanduse_{y}.tif,fkde_{y}.png)# 计算逐年转移矩阵并绘制弦图foriinrange(len(years)-1):bef,aftyears[i],years[i1]transtrans_matrix(flanduse_{bef}.tif,flanduse_{aft}.tif,class_names,ftrans_{bef}_{aft}.csv)draw_chord(trans,f{bef}-{aft}Cultivated Land Transition,fchord_{bef}_{aft}.html)# 地形因子分析需要DEMdem_pathdem_30m.tifslope_deg,reliefcalculate_slope_and_relief(dem_path)# 需定义该函数Tcalc_topographic_position(dem_path,slope_deg)# 保存地形位栅格withrasterio.open(dem_path)assrc:profilesrc.profile profile.update(dtyperasterio.float32)withrasterio.open(topographic_position.tif,w,**profile)asdst:dst.write(T.astype(np.float32),1)print(分析完成)七、总结本文基于贵州南北盘江流域20年的土地利用数据结合地形梯度分析和弦图可视化系统揭示了喀斯特山地流域耕地流转的时空规律。研究发现耕地总体收缩且向中低地形位集中流转主要发生在耕地、林地和草地之间地形梯度效应显著不同地形位的主导驱动因素各异。文中所附代码可实现核密度估计、地形位计算、转移矩阵和弦图绘制为类似研究提供了可复现的技术方案。关键词喀斯特山地耕地流转地形梯度核密度转移矩阵弦图字数统计约4300字含代码及注释。