¹曲阜师范大学地理与旅游学院，山东 日照

²现代古生物学和地层学国家重点实验室(中国科学院南京地质古生物研究所)，江苏 南京

收稿日期：2021年12月10日；录用日期：2022年1月31日；发布日期：2022年2月9日

摘要

基于黄河三角洲1998年、2008年和2018年的Landsat TM/OLI遥感影像数据，在ENVI5.1和ArcGIS10.2.2的支持下分析了黄河三角洲植被时空演化特征。研究结果表明：1998年到2018年这20年间本研究区的林地植被覆盖面积一直稳定在第三位，1998年植被覆盖面积第一位是芦苇及灌木植被，第二位是为耕地植被，2008年和2018年植被覆盖面积第一位变为耕地植被，第二位变为芦苇及灌木植被。从1998年到2018年，植被覆盖度面积变化最大的是极低植被覆盖度，面积减少了298.79 km²，其次是高植被覆盖度，面积增加了289.92 km²。该研究区NDVI总体上呈中部高，周边低，黄河河道、黄河故道高，沿海低的分布格局。较高等级的植被覆盖度主要分布在研究区内部，少量分布在河口地区，较低等级的植被覆盖度主要分布在北部、东南部的沿海滩涂上，以及水体、盐田及养殖池处，与NDVI分布趋于一致，较高等级的植被覆盖度由中部向西南、东南方向扩散，较低等级的植被覆盖度由东北部的盐田及养殖池向西扩散。通过转移矩阵计算得出，黄河三角洲1998年到2018年，芦苇及灌木植被、林地植被和耕地植被三者在空间上相互转移，其中芦苇及灌木植被和耕地植被之间的相互转移较为明显。

关键词

植被，时空演化，NDVI，植被覆盖度，黄河三角洲

Analysis of Temporal and Spatial Evolution Characteristics of Vegetation in the Yellow River Delta Based on Remote Sensing Image

Qingyue Xie¹, Xiaoping Xie^1,2

¹School of Geography and Tourism, Qufu Normal University, Rizhao Shandong

²State Key Laboratory of Paleontology and Stratigraphy (Nanjing Institute of Geology and Palaeontology, Chinese Academy of Sciences), Nanjing Jiangsu

Received: Dec. 10^th, 2021; accepted: Jan. 31^st, 2022; published: Feb. 9^th, 2022

ABSTRACT

Based on the Landsat TM/OLI remote sensing image data of the Yellow River Delta in 1998, 2008 and 2018, the temporal and spatial evolution characteristics of vegetation in the Yellow River Delta are analyzed with the support of ENVI5.1 and ArcGIS10.2.2. The results show that during the 20 years from 1998 to 2018, the forest land vegetation coverage in the study area has been stable in the third place. In 1998, the first vegetation coverage was reed and shrub vegetation, and the second was cultivated land vegetation. In 2008 and 2018, the first vegetation coverage became cultivated land vegetation, and the second became reed and shrub vegetation. From 1998 to 2018, the largest change in the area of vegetation coverage was the extremely low vegetation coverage, with an area decrease of 298.79 km², followed by the high vegetation coverage, with an area increase of 289.92 km². The NDVI in the study area is generally high in the middle, low around, high in the Yellow River channel and the old yellow river channel, and low along the coast. The higher level vegetation coverage is mainly distributed in the study area, and a small amount is distributed in the estuary area. The lower level vegetation coverage is mainly distributed on the coastal beaches in the north and southeast, as well as water bodies, salt fields and aquaculture ponds, which is consistent with the distribution of NDVI. The higher level vegetation coverage diffuses from the middle to the southwest and southeast. The lower vegetation coverage spreads westward from the salt field and aquaculture pond in the northeast. According to the calculation of transfer matrix, from 1998 to 2018, reed and shrub vegetation, forest vegetation and cultivated land vegetation in the Yellow River Delta transferred to each other in space, among which the mutual transfer between reed and shrub vegetation and cultivated land vegetation was more obvious.

Keywords:Vegetation, Temporal and Spatial Evolution, NDVI, Vegetation Coverage, Yellow River Delta

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

植被是覆盖地表的植物群落的总称，同时也是重要的自然资源，与自然环境中各因子之间存在着紧密的联系 [1]。植被作为环境保护者和生态圈正常健康运转的基础和条件，其作用可以有效地促进全球水循环、碳平衡以及保持生物的多样性，同时也对于提高国家和区域环境的质量以及改善人们的生活环境发挥着积极的作用 [2]。随着遥感技术的逐渐成熟以及遥感的宏观观测、动态监测、应用领域广泛等的特点，遥感技术越来越多的被应用于植被方面的研究。在上世纪80年代，国外学者讨论了遥感对土地利用研究的应用及遥感影像对土地覆被研究的贡献和潜力 [3]。国内学者刘耀林等人应用遥感技术与实地调查相结合的方法对武汉市武昌区城市用地进行了较系统地分类 [4]。许学工等人利用遥感和植被指数等方法研究了黄河三角洲地区数年至20余年的土地覆被质量变化 [5]。侯学会等人利用遥感影像分析了1992~2010年黄河三角洲自然保护区土地利用/覆被变化情况 [6]。植被的时空变化不仅体现了生态环境质量本身的变化，而且也体现了自然和人类对生态环境的作用过程 [7]。

本研究区位于黄河三角洲，地理位置特殊，生态环境较脆弱，受自然和人类活动因素的影响较大，植被的时空演化对该研究区的经济和生态环境意义重大。随着胜利油田开发及油田经济的飞速发展，城市化进程加快以及黄河河道的多次迁移，使得该地区植被呈现明显的时空变化规律。因此，本文以黄河三角洲植被时空演化为研究内容，利用2008年、2013年和2018年的Landsat TM/OLI遥感影像数据，对该研究区进行土地覆被分类、计算转移矩阵、动态度、归一化植被指数和植被覆盖度来研究黄河三角洲植被的时空演化特征，较前人研究更加完善。这有利于了解近20年来植被的空间和时间变化，就缓解该区域土地盐碱化和生态环境脆弱等问题提出相应的解决措施；有利于针对未来生态环境发展趋势制定合理的可持续发展战略，促进当地的生态环境和社会经济协调发展。

2. 研究区概述

本文研究区域位于山东省东营市东北部，渤海湾南岸与莱州湾西岸之间。地理坐标为东经118˚33'~119˚和北纬37˚35'~38˚12'之间。黄河三角洲是黄河入海所携带的泥沙在河口沉积形成的冲积平原，是黄河与渤海陆海相互作用的过渡地区。它以垦利县宁海为轴点，北起套尔河口，南至淄脉河口 [8]，是中国最年轻的土地之一 [9]。

地形以平原为主，地势由西南向东北倾斜，由于黄河在历史上的多次改道，地面略有起伏，多见坡地、河滩高地等微地貌景观，时常发生旱、涝、碱灾害。黄河三角洲沿海滩涂分布广阔，多为沙质和泥质，是重要的海产基地，包括黄河刀鱼、黄河鲤鱼、毛虾等。东营市境内已建设成了国家级自然保护区，以保护湿地生态系统、珍稀濒危鸟类为主体。

3. 研究方法及数据准备

3.1. 动态度模型

利用该模型可以分析研究区内各土地覆被类型在不同时期的变化速度和剧烈程度 [10]，其表达公式为：

$K=\frac{U_b-U_a}{U_a}\times \frac{1}{T}\times 100\%$ (1)

式中：K为研究时段内某一土地覆被类型动态度； $U_a$ 、 $U_b$ 分别为研究期初及研究期末某种土地覆被类型的面积；T为研究时间的范围。

3.2. 整体动态度模型

利用该模型可以分析整个研究区在一定时间范围内土地覆被类型变化的快慢情况 [11]，其表达公式为：

$LC=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}\Delta L{U}_{i-j}}}{2LU}\times \frac{1}{T}\times 100\%$ (2)

式中： $LU$ 为研究区总面积； $\Delta L{U}_{i-j}$ 表示研究时间范围内第i类土地覆被类型转变为非i类土地覆被类型面积的绝对值；T为研究时间的范围。

3.3. 转移矩阵模型

利用该模型可以对某地区各土地覆被类型内部转移进行分析 [12]，能够更好地认识该研究区内某种土地覆被类型的流失方向和补充来源。该模型为：

$A_{ij}=\left[\begin{array}{ccccc}{A}_{11}& A_{12}& A_{13}& \cdots & A_{1n}\\ A_{21}& A_{22}& A_{23}& \cdots & A_{2n}\\ A_{31}& A_{32}& A_{33}& \cdots & A_{3n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ A_{n1}& A_{n2}& A_{n3}& \cdots & A_{nn}\end{array}\right]$ (3)

3.4. 归一化植被指数

NDVI (归一化植被指数)是反映植被生长及空间分布密度的最有效指标 [13]，归一化植被指数是指近红外波段的反射值与红光波段的反射值之差比上两者之和 [1]，其取值范围为−1~1，负值表示地面覆盖为云、水、雪等，0值表示地面覆盖为岩石、裸地，正值表示地面有植被覆盖，且数值越大代表植被覆盖度越高。其表达公式为：

$NDVI=\left(NIR-R\right)/\left(NIR+R\right)$ (4)

式中：NDVI为归一化植被指数，NIR为近红外波段，R为红波段。

由于不同遥感影像的波段组成不同，因此在进行NDVI计算时要将波段与公式对应。所以，基于Landsat5 TM影像，公式(4)可转化为

$NDVI=\left(Band4-Band3\right)/\left(Band4+Band3\right)$ (5)

TM影像的Band4为近红外波段，Band3为红波段。

基于Landsat8 OLI影像，公式(4)可转化为

$NDVI=\left(Band5-Band4\right)/\left(Band5+Band4\right)$ (6)

OLI影像的Band5为近红外波段，Band4为红波段。

3.5. 植被覆盖度

植被覆盖度是衡量地表植被覆盖的一个重要指标 [14]，能够反映植被占地面面积的比例，其取值范围为0~1，数值越大，表示植被覆盖度越高。本研究采用像元二分模型，利用归一化植被指数提取植被覆盖度，其表达公式为：

$VFC=\left(NDVI-NDV{I}_{soil}\right)/\left(NDV{I}_{veg}-NDV{I}_{soil}\right)$ (7)

式中：VFC是植被覆盖度，NDVI_soil是无植被覆盖区域的NDVI值，NDVI_veg是纯植被覆盖的NDVI值。

在实际操作中受各种因素的影响，可将NDVI_soil和NDVI_veg近似等于NDVI_min和NDVI_max。所以，公式(7)可以转化为公式(8)。

$NDVI=\left(NDVI-NDV{I}_{min}\right)/\left(NDV{I}_{max}-NDV{I}_{min}\right)$ (8)

同时NDVI_min和NDVI_max的取值应在一定的置信区间，该研究区的NDVI_min和NDVI_max如表1所示。

表1. 研究区NDVI_min、NDVI_max统计表

3.6. 数据准备

本研究选取1998年9月10日、2008年10月7日Landsat5 TM和2018年10月3日Landsat8 OLI的遥感影像，空间分辨率为30 m，时间间隔为10年。选取遥感影像时尽量选择云量较低的影像，这样有利于遥感影像的处理进而达到更好的成像效果(表2) [15]。

表2. 研究区遥感影像

在ENVI5.1和ArcGIS10.2.2的支持下，对1998年、2008年和2018年这3期的遥感影像进行大气校正、几何校正、辐射定标、裁剪等预处理。根据先前的研究成果 [16] [17] [18] [19] [20]、相近年份黄河三角洲土地覆被类型图和黄河三角洲特点，本研究确定该研究区的土地覆被类型体系为：滩涂、水体、林地植被、盐田及养殖池、芦苇及灌木植被、耕地植被、城镇及建设用地和裸地。使用监督分类和传统目视解译及现场考察相结合的方法进行土地覆被类型分类。采用监督分类中的最大似然法进行分类，然后采用聚类处理(Clump Classes)对分类结果进行小斑块处理 [18]。运用ArcGIS10.2.2制作黄河三角洲1998年、2008年和2018年的土地覆被类型图(图1)。通过混淆矩阵方法进行精度评价，三幅土地覆被类型图的精度检验均在80%以上。运用上述公式和ENVI5.1软件，得到黄河三角洲归一化植被指数分布图(图2)。在归一化植被指数的基础上，利用像元二分模型计算植被覆盖度，并参照相关研究 [1] [2] 将植被覆盖度划分为5个等级(表3)，生成1998年~2018年植被覆盖度分布图(图3)。

图1. 1998~2018年土地覆被类型图

图2. 1998~2018 NDVI分布图

表3. 植被覆盖度等级划分

图3. 1998~2018植被覆盖度分布图

本研究区土地覆被类型中的滩涂主要指沿海大潮高潮位与低潮位之间的潮间带及河流湖泊低水位至常水位间的滩地；水体主要指该研究区内的黄河、水库和湖泊。盐田及养殖场主要指盐田和水产品养殖场；城镇及建设用地主要指城镇居民用地、农村居民用地、道路和油田开采用地；耕地植被主要是指水浇田植被和旱地植被；林地植被主要指人工刺槐林和天然林；芦苇及灌木植被主要包括芦苇、柽柳和翅碱蓬等覆盖面积较大的植被；裸地主要指未利用的土地 [19]。

4. 研究方法及数据准备

4.1. 黄河三角洲植被时间演化特征分析

根据黄河三角洲的3景土地覆被类型图、植被覆盖分布图，分别计算出每年土地覆被类型的面积(图4)、比例(表4)、面积变化、变化率和动态度(表6)，植被覆盖度不同等级的面积(图5)、面积变化和动态度(表5)。

图4. 黄河三角洲各年土地覆被类型面积(单位：km²)

表4. 黄河三角洲各年土地覆被类型比例

根据图4和表4可以看出，1998年芦苇及灌木植被覆盖面积最大，所占面积达到1287.80 km²，占研究区土地覆被类型比例37.72%；其次是耕地植被，面积为689.76 km²，占研究区土地覆被类型比例22.69%；面积最小的林地植被，占地面积为54.89 km²。2008年耕地植被面积大量增加，增加到1166.50 km²，占研究区土地覆被类型比例39.24%；其次是面积为552.61 km²的芦苇及灌木植被，占研究区土地覆被类型比例18.59%；位于第三位的植被是林地植被，占研究区土地覆被类型比例1.47%。2018年植被覆盖面积最大的是耕地植被，面积为1117.66 km²，占研究区土地覆被类型比例36.21%；第二位的是芦苇及灌木植被，占研究区土地覆被类型比例12.87%；第三位是林地植被，占研究区土地覆被类型比例2.82%。总体来看，1998年到2018年这20年间本研究区的林地植被覆盖面积一直稳定在第三位，1998年植被覆盖面积第一位是芦苇及灌木植被，第二位是为耕地植被，2008年和2018年植被覆盖面积第一位变为耕地植被，第二位变为芦苇及灌木植被。

图5. 1998~2018各等级植被覆盖度面积(单位：km²)

由图5和表5可知，1998年高植被覆盖度所占面积最大，为896.85 km²，其次是极低植被覆盖度，为873.39 km²，面积最小的是低植被覆盖度，为392.56 km²；2008年面积最大的是1347.45 km²的高植被覆盖度，其次是低植被覆盖度，最小的是极低植被覆盖度；2018年高植被覆盖度所占面积仍然最大，为1186.77 km²，其次是极低植被覆盖度，第三位的是中高植被覆盖度，最后依次相对较少的为中植被覆盖度、低植被覆盖度。1998年到2008年，植被覆盖度面积变化最大的是极低植被覆盖度，以6.38%的速度减少了557.24 km²，其次是高植被覆盖度，变化率为50.24%，动态度为5.02%，面积增加了450.60 km²；2008年到2018年，植被覆盖度面积变化最大的是极低植被覆盖度，以8.17%的速度增加了258.44 km²，其次是低植被覆盖度，动态度为−4.17%；从1998年到2018年，植被覆盖度面积变化最大的是极低植被覆盖度，变化率为−34.21%，动态度为−1.71%，面积减少了298.79 km²，其次是高植被覆盖度，以1.62%的速度增加了289.92 km²。

表5. 1998~2018年黄河三角洲各植被覆盖度面积变化及动态度

表6. 1998~2018年黄河三角洲土地覆被类型面积变化及动态度

根据图1和表6可以看出，1998年到2018年黄河三角洲的植被覆盖面积存在变化。1998年到2008年，研究区的整体动态度为2.01%，2008年到2018年，研究区的整体动态度为1.66%，1998年到2018年，研究区的整体动态度为1.34%。从耕地植被面积变化来看，1998年到2008年，变化率为69.12%，呈增加趋势，动态度为6.19%，说明耕地植被以6.19%的速度增加了476.74 km²；2008年到2018年，耕地植被的变化率为−4.19%，呈下降趋势，但下降的幅度较小，动态度为−0.42%，耕地植被以0.42%的速度下降了48.84 km²；1998年到2018年，耕地植被以3.10%的速度增加了427.9 km²；所以在此二十年间，耕地植被的面积在上升–下降的趋势中增加。从芦苇及灌木植被面积变化来看，1998年到2008年，面积变化率为−51.81%，呈下降趋势，面积以5.18%的速度减少了594.23 km²；2008年到2018年，芦苇及灌木植被的变化率为−28.09%，呈下降趋势，动态度为−2.81%；1998年到2018年，芦苇及灌木植被以3.27%的速度减少了749.44 km²；所以在这二十年中，芦苇及灌木植被的面积在下降–下降的趋势中减少。从林地植被面积变化来看，1998年到2008年，变化率为−20.33%，呈下降趋势；2008年到2018年，林地植被的变化率为98.93%，上升趋势较大；1998年到2018年，林地植被以2.92%的速度增加了32.1 km²；因此在这二十年间，林地植被面积在下降–上升的趋势中增加。

4.2. 黄河三角洲植被空间演化特征分析

根据图1(A)~(C)这3景黄河三角洲土地覆被类型图，我们可以看出，耕地植被主要分布于黄河三角洲内部的大部分地区，靠近黄河两岸及建设用地附近。芦苇及灌木植被1998年分布面积最广，三角洲内部大量分布，2008年和2018年三角洲内部的芦苇及灌木植被大量减少，主要分布在老河口附近以及现黄河入海口两岸，呈带状分布，同时在黄河三角洲内部的城镇及建设用地和耕地植被中零散分布。林地植被主要是刺槐林，现代黄河三角洲刺槐林是纯人工林，是现代黄河三角洲森林的主体。因为林地植被为人工林，所以在靠近河口区域的自然保护区人类活动较少因而分布面积较少 [19]，但1998年到2018年这二十年间，林地植被有增加的趋势，主要分布于自然保护区黄河口管理站黄河北岸中心路以西和黄河故道两侧，与前人的研究结果一致 [16]。由图2可知，该研究区NDVI总体上呈中部高，周边低，黄河河道、黄河故道高，沿海低的分布格局，这与前人的研究结果一致 [15]。NDVI的负值主要分布在水体、沿海的盐田及养殖池，在0~1区间的主要分布在林地、耕地、芦苇及灌木植被。由图3可以看出，较高等级的植被覆盖度主要分布在研究区内部，少量分布在河口地区，较低等级的植被覆盖度主要分布在北部、东南部的沿海滩涂上，与路广等人的研究结果一致 [1]，以及水体、盐田及养殖池处，与NDVI分布趋于一致，较高等级的植被覆盖度由中部向西南、东南方向扩散，较低等级的植被覆盖度由东北部的盐田及养殖池向西扩散。

以1998年的土地覆被类型图为前一时相，以2008年的土地覆被类型图为后一时相 [21]，在ENVI5.1的支持下，利用Change Detection Statistics功能计算1998年到2008年黄河三角洲土地覆被类型转移矩阵(表7)。

表7. 黄河三角洲1998~2008土地覆被类型转移矩阵(单位：km²)

结合以上分析和表7可知，1998年到2008年，芦苇及灌木植被大量减少，主要转出为耕地植被，有423.46 km²的芦苇及灌木植被转变为了耕地植被，其次，有101.60 km²的芦苇及灌木植被转变为了盐田及养殖池，芦苇及灌木植被的转入主要来自水体、滩涂、耕地植被和裸地，但转入较少，所以这也是在这十年间该植被大量减少的主要原因。林地植被的转出主要为耕地植被，有15.95 km²转变为耕地植被。耕地植被的主要转出方向为城镇及建设用地、芦苇及灌木植被，分别转出了13.82 km²和29.74 km²，其转入主要来自芦苇及灌木植被，有423.46 km²的芦苇及灌木植被转入为耕地植被，耕地植被转出面积较少，因此，在该十年耕地植被面积大量增加。总体来说，这十年间芦苇及灌木植被覆盖面积大量减少，耕地植被覆盖面积大量增加，同时两者之间大量转移。

以2008年的土地覆被类型图为前一时相，以2018年的土地覆被类型图为后一时相，在ENVI5.1的支持下，利用Change Detection Statistics功能计算2008年到2018年黄河三角洲土地覆被类型转移矩阵(表8)。

表8. 黄河三角洲2008~2018土地覆被类型转移矩阵(单位：km²)

结合以上分析和表8可以得知，2008年到2018年，芦苇及灌木植被依然呈减少趋势，主要转出方向与1998年到2008年基本一致，主要向耕地植被、盐田及养殖池转出，分别转出了121.01 km²和110.36 km²，其转入主要来自滩涂、水体、耕地植被和裸地。林地植被的转出较少，有37.87 km²的耕地植被向林地植被转入。耕地植被主要由芦苇及灌木植被、裸地转入，转出方向是主要是盐田及养殖池、城镇及建设用地、芦苇及灌木植被和林地植被。

从整体上来看，黄河三角洲1998年到2018年，芦苇及灌木植被、林地植被和耕地植被三者在空间上相互转移，其中芦苇及灌木植被和耕地植被之间的相互转移较为明显。随着经济的发展，该研究区的城镇及建设用地和盐田及养殖池面积逐年增加，大量芦苇及灌木植被向盐田及养殖池、耕地植被转移，造成芦苇及灌木植被大量减少，耕地植被面积大幅增加。

5. 结论

本文对20年的3景黄河三角洲遥感影像数据进行了解译分析，提取了该地区的土地覆被类型、NDVI、植被覆盖度，对黄河三角洲植被时空演化特征进行了分析，结论如下：

1) 黄河三角洲植被时间演化特征：1998年到2018年这20年间本研究区的林地植被覆盖面积一直稳定在第三位，1998年植被覆盖面积第一位是芦苇及灌木植被，第二位是为耕地植被，2008年和2018年植被覆盖面积第一位变为耕地植被，第二位变为芦苇及灌木植被。从1998年到2018年，植被覆盖度面积变化最大的是极低植被覆盖度，面积减少了298.79 km²，其次是高植被覆盖度，以1.62%的速度增加了289.92 km²。1998年到2018年黄河三角洲的植被覆盖面积存在变化，耕地植被的面积在上升–下降的趋势中增加，芦苇及灌木植被的面积在下降–下降的趋势中减少，林地植被面积在下降–上升的趋势中增加。

2) 黄河三角洲植被空间演化特征：耕地植被主要分布于黄河三角洲内部的大部分地区，靠近黄河两岸及建设用地附近。芦苇及灌木植被1998年分布面积最广，三角洲内部大量分布，2008年和2018年三角洲内部的芦苇及灌木植被大量减少，主要分布在老河口附近以及现黄河入海口两岸，呈带状分布，同时在黄河三角洲内部的城镇及建设用地和耕地植被中零散分布。林地植被有增加的趋势，主要分布于自然保护区黄河口管理站黄河北岸中心路以西和黄河故道两侧。该研究区NDVI总体上呈中部高，周边低，黄河河道、黄河故道高，沿海低的分布格局。较高等级的植被覆盖度主要分布在研究区内部，少量分布在河口地区，较低等级的植被覆盖度主要分布在北部、东南部的沿海滩涂上，以及水体、盐田及养殖池处，与NDVI分布趋于一致，较高等级的植被覆盖度由中部向西南、东南方向扩散，较低等级的植被覆盖度由东北部的盐田及养殖池向西扩散。通过转移矩阵计算得出，黄河三角洲1998年到2018年，芦苇及灌木植被、林地植被和耕地植被三者在空间上相互转移，其中芦苇及灌木植被和耕地植被之间的相互转移较为明显。

文章引用

谢清月,谢小平. 基于遥感影像的黄河三角洲植被时空演化特征分析
Analysis of Temporal and Spatial Evolution Characteristics of Vegetation in the Yellow River Delta Based on Remote Sensing Image[J]. 地理科学研究, 2022, 11(01): 73-85. https://doi.org/10.12677/GSER.2022.111009

参考文献