基于RS和GIS五常市土地利用变化及模拟研究

龚文峰¹, 王鑫鑫¹, 曹凯华¹, 周 莹², 吴 娟¹

(1.黑龙江大学 水利电力学院，哈尔滨 150080 ;2.黑龙江国防科学技术研究院，哈尔滨 150040)

摘 要： 以五常市2005、2010、2015年LandsatTM影像为主要数据源，在RS和GIS支持下，借助土地利用时空演变模型，探讨研究区域10 a间土地利用格局时空演化规律和特征，基于CA-Markov模型模拟研究区域2020年土地利用空间格局，以期为其土地资源合理利用和可持续发展提供理论依据和技术参考。结果表明：①2005—2015年耕地增加394.91 km²，居民建筑用地增加74.81 km²，林地减少517.49 km²，林地→耕地是主要土地利用转变类型；②耕地保留率高，林地次之；③2015年土地利用模拟结果和实际解译数据十分吻合，Kappa系数为0.900 6。

关键词： 土地利用；空间分析模型；CA-Markov模型；五常市

0 引 言

1995年IGBP(国际地圈—生物圈计划)和IHDP(国际全球环境变化人文因素计划)共同制定了“土地利用变化”研究计划^[1-2]，2005年将其发展为GLP(全球土地计划)^[3]。土地利用变化是全球环境变化的重要组成部分和主要原因之一^[4]，也是当前研究全球变化的核心课题和热点之一^[5-6]。全球土地利用变化不仅需要大尺度、大规模的宏观研究，也需要中小尺度研究加以配合完善。同时，中国是全球变化研究的重要组成部分，但近年来土地利用变化研究主要集中在经济高速发展的城市周边、生态环境脆弱地带及河流湖泊等区域^[7-13]，对中国东北地区县级尺度的土地利用变化研究相对较少。

土地利用变化传统的监测模型有CA模型、CLUE-S模型、GTR模型、SD模型、Markov模型以及多智能体模型^[14-19]等。Markov模型具有长期预测的特点，不具备空间格局变化预测的能力；CA模型只注重元胞周围空间区域的相互影响，具有一定局限性。CA-Markov模型结合了CA模型空间运算和Markov模型长期预测的特点，该模型广泛应用于土地利用变化模拟研究中，但对于东北地区粮食主产区的研究较少，尤其基于典型黑土区的研究更少。

黑龙江省五常市是东北地区典型黑土区域^[20]，以盛产“五常稻花香大米”而出名。然而，近年市场对稻米的需求渐增，在一定程度上刺激了耕地的不合理扩张，在追求经济利益的过程中，给区域农业经济可持续发展带来严重威胁^[21]。基于此，本研究以五常市为研究对象，以2005、2010和2015年Landsat数据为主要数据源，借助于RS和GIS技术，获取研究区土地利用数据，基于土地利用时空演变模型，定量分析土地利用时空演变特征；借助于CA-Markov模型，定量模拟2020年土地利用空间格局，以期为研究区土地利用结构调整、布局优化及区域可持续发展提供理论依据和技术参考。

1 研究区概况

五常市位于黑龙江省南部，其地理坐标为：E126°33′～128°14′和N 44°04′～45°26′，属中纬度温带大陆性季风气候。其西北接松嫩平原，东南靠张广才岭西麓，东北部与阿城市相邻，西南部、南部与吉林省接壤。地理构造处于大兴安岭褶皱山带和长白山、小兴安岭褶皱山带中间的台地上。境内水系发达，河网密布，水资源丰富，有拉林河、阿什河两大水系。五常市夏短冬长，温差较大，年均气温3～4 ℃，全年无霜期130 d左右，平均年降水量625 mm。全市土壤种类有11个土类33个亚类，主要土壤有暗棕色森林土、白浆土、黑土、草甸土、沼泽土和水稻土，其中黑土主要分布在漫岗和坡地上，土层厚，基础肥力高，适合于各种农作物的生长^[20]。

2 数据来源与研究方法

2.1 数据来源与处理

本研究所用遥感数据为五常市2005年和2010年的Landsat7 ETM数据及2015年Landsat8 OLI数据(http://landsat.usgs.gnv/)、30 m分辨率的DEM数据、行政区划矢量数据以及2005—2015年社会发展规划数据。在遥感软件ENVI5.3支持下，完成遥感图像的辐射定标、大气校正、几何校正、投影转换、镶嵌和裁剪等处理。

根据中国土地利用分类体系，结合研究区土地利用现状特点，借助ENVI5.3遥感处理软件，采用监督分类和人机交互相结合的方法完成遥感影像的分类解译，将土地利用类型分为耕地、林地、草地、水域、居民建筑用地及未利用地6大类，利用Google Earth和部分样点数据，应用混淆矩阵验证分类精度，其Kappa系数分别为0.956 2、0.878 0、0.900 6。由Kappa系数分类标准(表1)，一致性程度为几乎一致，表明分类精度高，在ArcGIS的支持下，将分类数据以GRID的格式存储(图1)。

表1 Kappa系数分类标准

Table 1 Kappa coefficient classification standard

图1 五常市2005、2010和2015年土地利用类型图
Fig.1 Land use maps of Wuchang City in 2005、2010、2015

2.2 研究方法

2.2.1 土地利用空间动态演变模型

式中n 为研究区土地利用类型的数目；P _ij 为初期到末期时段内土地类型i 转为类型j 的转移概率。

本文借助ArcGIS10.3的空间分析功能，通过对不同时期土地利用数据叠加分析处理，定量探究研究区土地利用时空演变规律，揭示土地利用类型转变过程和形成机制。公式如下：

(1)

(2)

CCL _i =TRL _i +IRL _i

(3)

基于ArcGIS完成五常市3个时期土地利用类型面积统计(表2)。由表2可见，2005—2010—2015年，五常市土地利用主要以耕地、林地为主，两者之和占总面积比分别为94.8%、92.68%和93.17%，说明研究区域土地利用结构以农林生产为主，符合该区域经济和社会发展的特点。2005—2010年和2010—2015年，耕地和居民建筑用地面积呈现持续增加的趋势，其中，耕地面积增加206.39 km²和188.52 km²，主要原因是受到粮食所带来的经济效益的刺激，农户加大开垦耕地力度；居民建筑用地扩张明显，面积增加70.13 km²和4.68 km²，主要原因是受城市化进程加速的影响，城镇范围逐步扩大；林地面积呈现递减的趋势，面积减少366.29 km²和151.2 km²，林地面积减少主要源于农户为了眼前利益毁林开荒，大片林地遭到破坏；水域呈先增后减趋势，2005—2010年，面积增加19.44 km²，主要源于2006年磨盘山水库建成，淹没了周围的耕地、林地和草地，使得短时间内水域面积快速增加，但2010—2015年受气候变化影响，水域面积略有降低，减少0.76 km²；草地呈先增后减趋势，2005—2010年面积增加78.98 km²，2010—2015年面积减少55.01 km²；未利用地呈先减后增趋势，2005—2010年面积减少8.65 km²，2010—2015年增加13.77 km²。

2.2.2 土地利用转移矩阵

式中X _t+1 为t +1时刻土地利用类型状态；X _t 为t 时刻土地利用类型状态；P _ij 为土地利用变化状态转移概率。

(4)

基于以上对创新领域的发展成效分析，可知美国创新战略具有以下4方面的特点：一是美国致力于经济的可持续发展，通过推进原始创新带来实体经济的长期稳定增长；二是美国高度重视基础研究带来的长期效益，不断加大核心技术研发投入；三是2015年美国创新战略再次强调了政府在推动创新发展中的重要作用，强化了政府创新服务职能；四是美国重视人才培养，大力提倡高校积极参与科技创新活动。

2.2.3 CA-Markov模型

Markov模型计算公式如下：

X _t+1 =P _ij X _t

(5)

土地利用转移矩阵能够描述各土地利用类型之间的转化情况，还能反映土地利用变化的结构特征和各类型间的转移方向。本文基于IDRISI的Markov模块得到2005—2010年和2010—2015年土地利用类型转移矩阵，分析区域土地利用类型的动态变化特征。数学表达式如下：

删除训练集中已有的样本时，采用减量式(或“遗忘”)算法。当样本xc在集合R中，如果它对SVR的解无作用，将其从训练集中移除时影响不大，且不需要调整。另一方面，若xc有一个非零系数，则将这个非零系数的值逐渐减少为0，同时确保训练集中的其它样本依然满足KKT条件。

当前废旧农膜捡拾仍然以人工捡拾为主，1个农民每天的捡拾量为0.2 hm2，可捡拾废旧农膜15～20 kg（折纯），按照1.2元/kg的最高收购价回收，农民通过捡拾废旧农膜的收入为18～24元/d，如果将运输成本考虑进去的话，收入更低。捡拾废旧农膜的收入与打工收入相差太大，农民对捡拾废旧农膜没有参与积极性。

首先，高等教育方面。宿迁市在一本教育上出现断层，目前只有宿迁学院一所民办本科院校为宿迁的经济发展输出人才。虽然其设有电子商务、工商管理以及金融等专业，但其输出人才的数量远不能满足宿迁市经济发展的需要。从GDP角度来看，宿迁市在江苏省属于经济发展排名靠后的城市，难以吸引高学历、高层次的优秀人才流入。

P _t =f (P _t+1 ,N )

(6)

式中P 为元胞状态集合；N 为邻域；t 和t +1为时刻；f 为元胞状态转换规则。

CA-Markov模型具体实现过程如下^[20]：①利用Markov模型获取2005—2010、2010—2015年土地利用转移概率矩阵，并与MCE构建适宜性图像集；②构造CA滤波器，采用5×5滤波器，并将每个邻域元胞对中心元胞的影响权重设为相同；③确定起始时刻和CA循环次数。以2010—2015年土地利用转移矩阵为基础，将2015年作为土地利用空间格局预测的起始时刻，迭代次数设置为5，模拟2020年土地利用空间格局。

模板在不同厚度条件下，外墙保温围护结构的平均传热系数K0与插丝率之间的变化关系如图5所示。从图5可以看出，5条曲线均呈现逐渐递增趋势，即平均传热系数K0均随着插丝率的增大而增大。对于岩棉厚度为150 mm的岩棉复合型保温模板，当插丝率为0.05%时，外墙保温围护结构的平均传热系数K0为0.355 W/(m2·K)，由此可知，岩棉复合型保温模板保温体系的热工性能十分优异。由计算结果可知，岩棉复合型保温模板的保温性能完全可以满足我国现行的墙体最高节能率为65%的节能要求。

3 结果与分析

3.1 土地利用总体特征变化

式中LA _i,t ₁-ULA _i 为在监测期间转出部分面积，即第i 种土地利用类型转化为其他非i 类土地利用类型的面积总和；LA _i,t ₁为监测期初第i 种土地利用类型面积；ULA _i 为监测期间第i 种土地利用类型未变化部分面积；LA _i,t ₂-ULA _i 为监测期间转入部分面积，即在监测期间非i 类土地利用类型转入i 类土地利用类型面积；LA _i,t ₂为监测期末第i 种土地利用类型面积；TRL _i 为转出速率；IRL _i 为转入速率;CCL _i 为第i 种土地利用类型在监测期t ₁～t ₂空间动态变化速率。CCL _i 值高低反映了地类土地利用变化的活跃和稳定程度。

表2 2005,2010,2015年土地利用分类统计

Table 2 Area statistics of lands use and cover change in 2005,2010 and 2015

3.2 土地利用类型动态变化

在ArcGIS的支持下，分别对2005—2010和2010—2015年的土地利用数据进行空间叠加分析，计算不同时段土地利用类型的未变化面积、转出面积、转入面积、转出速率、转入速率和变化速率。

小王在制订目标之后，将操作行为界定为练习难度和求教。我们设计了这样的操作行为：问题的寻找和设计，问题的练习和实践，问题的整理和反思，同伴的合理讨论和师生间的讨论。这五方面就是具体的可操作过程，同时这个五个方面彼此递进，前三个方面是问题的实践，从选择、实践、整理，层层深入，后两个方面重在求教，重视同伴互助和寻找老师的帮忙。这里的具体化和层次性是融合在一起的，最终成为达成目标的种子。

研究区2005—2010年和2010—2015年土地利用动态变化率见表3、表4。由表3、表4可见，2005—2010年，研究区土地利用类型变化速率的排序为：草地>未利用地>居民建筑用地>水域>耕地>林地，2010—2015年变化速率排序为：未利用地>草地>居民建筑用地>水域>耕地>林地。由此可见，未利用地、草地变化速率大，水域、耕地次之，林地最小，表明未利用地和草地较敏感，当受到人类活动的干扰，极易转变为其它地类。2005—2010—2015年，居民建筑用地转入速率远大于同期的转出速率，属高速扩展型，变化速率由26.80%降至19.98%，活跃程度降低，源于城市基础建设已达到一定规模，城市化进程放缓，居民建筑用地面积趋于稳定；耕地转入面积最大，两期转入面积分别为608.48 km²、518.01 km²，占区域所有土地利用类型转入面积总和的50.95%和47.25%，是其他地类转出的主要方向，且转入速率大于其同期的转出速率，属扩展型，但由于初期面积大，耕地变化速率整体较低，由5.98%降至4.72%；林地的转出面积最大，分别为559.48 km²、439.30 km²，是其他地类转入面积的主要来源，且转出速率远大于其同期的转入速率，属于高速衰减型，但变化速率整体仍比较低。研究发现，低海拔、缓坡区域林地是耕地增加主要来源，分别向耕地转出面积为461.45 km²和321.48 km²，因此需加强“退耕还林”，促进区域生态健康发展；水域呈先增后减趋势，变化速率由15.97%降至14.40%，下降了1.57%，原因是磨盘山水库建成，作为重要的饮用水水源地，水域受到一定程度的保护。

元胞自动机(Cellular Automata，CA) 是上世纪 40年代由Ulam首次提出的一种网格动力学模型，其计算公式为：

（4）综合性实验、实训教学的应用与管理。将一系列单一项目按其内在关系有机结合在一起进行的综合性实验、实训教学，是培养和提升学生综合素质和综合能力的主要途径和必然选择，应予以广泛倡导与运用。

表3 研究区2005—2010年土地利用动态变化率

Table 3 Lands use dynamic change rate from 2005 to 2010 in study area

表4 研究区2010—2015年土地利用动态变化率

Table 4 Lands use dynamic change rate from 2010 to 2015 in study area

研究区2005—2010年，2010—2015年土地利用转换概率见表5、表6。由表5、6可见，耕地、林地、水域、居民建筑用地、未利用地、草地保留率分别为88.11%、90.82%，84.99%、86.93%，73.84%、63.57%，51.50%、50.59%，8.99%、11.88%，4.34%、12.50%。10 a间，耕地、林地保留率高，水域和居民建筑用地次之，草地和未利用地低；耕地向林地转移面积比例为3.88%、5.05%，主要源于“退耕还林”政策的推行。耕地向建筑用地转移面积比例为3.99%、2.59%，原因是城市化进程推进，部分耕地被居民建筑用地侵占；林地向耕地转移的面积比例为 12.38%、9.56%，主要原因是受粮食经济效益刺激，农户破坏耕地周边的缓坡林地，种植农作物；居民建筑用地向耕地转移面积比例为45.52%和35.75%，源于政府对一些农村居民点“空心村”的整治，在一定程度上补给耕地资源；水域前期主要向耕地、林地转移，转移面积比例为13.44%和10.40%，后期主要向林地转移，转移面积比例为21.88%；草地和未利用地面积比例小，变化幅度大。草地向耕地、林地转移面积比例分别为41.52%、49.28%和53.20%、33.23%。未利用地向居民建筑用地、耕地转移了46.42%、30.30%和33.05%、81.82%，主要原因是未利用地开发程度低，易受人类活动影响变为城市建设用地或开垦为耕地。

表5 研究区2005—2010年土地利用转移概率表

Table 5 Normalized land use/cover transition probabilities from 2005 to 2010 in study area %

表6 研究区2010—2015年土地利用转移概率表

Table 6 Normalized land use/cover transition probabilities from 2010 to 2015 in study area %

3.3 土地利用变化模拟

基于IDRISI 17.0的CROSSTAB模型，以2010年的数据为基准数据，借助于2005—2010土地利用转移概率，模拟2015年土地利用数据(图2)，计算模拟精度(Kappa=0.900 6)，该结果表明模拟结果与实际现状图匹配度高。基于此，以2015年的土地利用现状为起始数据，基于2010—2015年土地利用转移概率矩阵，运行CA-Markov模型，模拟获取2020年土地利用格局^[20](图2)。

图2 土地利用模拟图
Fig.2 Simulation of land use map

模拟结果(表7)显示：2020年耕地和林地的面积比例分别为53.89%和40.23%，仍是研究区域的主要土地利用类型。与2015年相比，耕地和居民建筑用地依然呈现扩张的趋势。其中耕地面积由2015年的3 777.55 km²增至4 041.62 km²，比《五常市国民经济和社会发展第十三个五年规划纲要》规划的面积增加1 134.36 km²，因此，对于耕地今后的发展是严格控制耕地的不合理扩张，在确保耕地红线的前提下，大力推进“退耕还林还草”等生态措施，促进区域土地资源的合理利用；居民建筑用地面积由2015年的264.45 km²增至2020年的285.32 km²，同规划面积相比，居民建筑用地面积减少66.11 km²，在一定程度上表明居民建筑用地还有一定的扩展空间，今后应继续积极推进城镇基础设施建设，优化发展城镇空间，促进城市化进程中土地利用的合理化；与模拟和规划面积相比，林地、未利用地和水域呈现萎缩趋势，今后在加强对林地保护的同时，更应积极推行“植树造林”和“退耕还林”等生态措施，加大弃耕地和坡耕地的改造力度，避免对林地进一步的破坏，通过人为措施增加林地的面积，促进区域土地合理利用。

亚洲乳水牛的养殖主要集中在印度、巴基斯坦及中国，存栏量分别为4 182万头、1 334万头和559万头，水牛乳产量分别为7 800万t、2 564万t、313万t。自2001年世界水牛大会将水牛的主要用途确定为乳用为主后，亚洲水牛乳产量保持着每年2.5%～3.5%的增长速率，其中印度和巴基斯坦的水牛乳产量增长最为显著，至2016年两国水牛乳产量分别占全国乳类总产量的48%和64%，超过荷斯坦牛乳，成为当地人民的首要乳源。

表7 2020年五常市土地利用分类模拟结果与规划目标对比表

Table 7 Comparison table of the simulation and planning in 2020

4 结 论

本文在RS和GIS支持下，以五常市2005、2010、2015年LandsatTM影像为主要数据源，借助土地利用时空演变模型，探究研究区2005—2015年土地利用格局时空演化规律和特征，并基于CA-Markov模型模拟研究区2020年土地利用空间格局，结果表明：

1)2005—2015年，五常市土地利用结构主要以耕地和林地为主，土地利用时空演变呈现“五增一减”的特点，耕地、居民建筑用地、草地、水域及未利用地呈增长趋势，林地呈下降趋势。

2)2005—2015年主要地类转移为：耕地↔林地、耕地↔居民建筑用地和草地→耕地，其中林地→耕地变化最为显著。耕地保留率分别为88.11%和90.82%，林地保留率分别为84.99%和86.93%。

3)基于CA-Markov模型模拟获取到2020年研究区土地利用格局，耕地、林地、居民建筑用地、未利用地、草地、水域的面积分别为4 041.62、3 016.89、285.32、28.24、54.43、73.1 km²。

参考文献：

[1] Howells M, Hermann S, Welsch M, et al. Integrated analysis of climate change, land-use, energy and water strategies [J]. Nature Climate Change,2013,3(7):621-626.

[2] Wright C K, Wimberly M C. Recent land use change in the Western Corn Belt threatens grasslands and wetlands [J]. Proc Natl Acad Sci U S A, 2013,110 (10):4134-4139.

[3] 王让虎,李晓燕,李艳博,等.东北农牧交错带土地利用/覆被变化及驱动因素分析—以吉林省通榆县为例[J].广东农业科学,2013,40 (21):228-232.

[4] 程洁如.梅州市城市化发展与土地利用研究[J].经济地理,2009,29(10):1665-1670.

[5] 蔡运龙.土地利用/土地覆被变化研究:寻求新的综合途径[J].地理研究,2001,20(6):645-652.

[6] 陈佑启,杨鹏.国际上土地利用/土地覆盖变化研究的新进展[J].经济地理,2001,21(1):95-100.

[7] 张津,朱文博,吴舒尧,等.基于CLUE-S模型的京津冀城市群土地利用变化时空模拟[J].北京大学学报:自然科学版,2018,54 (1):115-124.

[8] 王涛,李贝贝,何亮,等.西安城市扩展时空特征与驱动因素分析[J].测绘科学,2017,42(4):75-79，118.

[9] 张军涛,李颖.近10年来东北农牧交错区土地利用变化研究[J].地理科学进展,2003,22 (6):551-559.

[10] 许小娟,刘会玉,林振山,等.基于CA-Markov模型的江苏沿海土地利用变化情景分析[J].水土保持研究,2017, 24 (1): 213-218.

[11] 龚文峰,王鹏,王双宇,等.基于 GF-2卫星遥感影像的界河水体信息提取方法[J].黑龙江大学工程学报,2018,9(4):1-7.

[12] 李京京.黄土高原地区土地利用/覆被变化及其驱动力分析研究[D].杨凌:西北农林科技大学, 2017.

[13] 魏光辉，姜振盈，张洛晨.新疆和田河流域1990—2015年土地利用/覆被复化及驱动因素研究[J].黑龙江大学工程学报, 2018，9(3)：14-20.

[14] 冯永玖,杨倩倩,崔丽,等.基于空间自回归CA模型的城市土地利用变化模拟与预测[J].地理与地理信息科学, 2016,32(5):37-44.

[15] 李桢,刘淼,薛振山,等.基于CLUE-模型的三江平原景观格局变化及模拟[J].应用生态学报,2018,29(6):1805-1812.

[16] Smith N, Dennis W.The restructuring of geographical scale: coalescence and fragmentation of the Northern Core Region [J]. Economic Geography,1987,63(2):160-162.

[17] 秦贤宏,段学军,李慧,等.基于SD和CA的城镇土地扩展模拟模型—以江苏省南通地区为例[J].地理科学,2009, 29(3):439-444.

[18] 陈颐.基于Markov和Logistic模型的莆田市土地利用变化及林地转出空间模拟[J].中国农业大学学报,2017, 22(2):87-97.

[19] 宋世雄,梁小英,陈海,等.基于多智能体和土地转换模型的耕地撂荒模拟研究-以陕西省米脂县为例[J].自然资源学报,2018,33(3):515-525.

[20] 吴娟.黑土区县级尺度土地利用变化及空间格局模拟研究[D]. 哈尔滨:黑龙江大学,2017.

[21] 王文娟,邓荣鑫,张树文.东北典型黑土区讷谟尔河流域土地利用变化对沟蚀影响研究[J].自然资源学报,2016, 31(5):833-844.

Research on land use change and its simulation in Wuchang City based on RS and GIS

GONG Wen-Feng¹ ,WANG Xin-Xin¹,CAO Kai-Hua¹,ZHOU Ying², WU Juan¹

(1.College of Hydraulic and Electric Power ,Heilongjiang University ,Harbin 150080,China ;2.Heilongjiang Institute of National Deference Science and Technology ,Harbin 150040,China )

Abstract ：The study took the LandsatTM images of 2005, 2010 and 2015 in Wuchang City as the main data sources, with the support of RS and GIS technique, using the spatial and temporal evolution model, the temporal and spatial evolution law and characteristics of land use pattern in the study area were analyzed, on this basis, the CA-Markov model was applied to simulate the spatial pattern of land use in 2020 in this paper ,which provided theoretical basis and technical supports for the rational use and sustainable development in the study area. The results showed that:①During 2005—2015, the cultivated land and residential land increased to 394.91 km² and 74.81 km² ,respectively, the area of the forest land decreaed to 517.49 km². The transformation from forest land to cultivated land was the major land use change of the study area;②Retention rate of cultivated land was higher than that of other land use, followed by forest land;③Land use simulation in 2015 were in good agreement with the actual interpretation data with the kappa coefficient of 0.900 6.

Key words ：land use; spatial analysis model; CA-Markov model; Wuchang City

DOI: 10.13524/j.2095-008x.2019.02.018

中图分类号： F301.24

文献标志码： A

文章编号： 2095-008X(2019)02-0015-08

收稿日期： 2019-01-26;修订日期：2019-02-25

基金项目： 黑龙江省自然科学基金资助项目(QC2017036)

作者简介： 龚文峰(1976-)，男，河南南阳人，教授，博士，硕士研究生导师，研究方向：遥感与GIS应用、水环境遥感监测与土地利用等。E-mail:gwf101@163.com

标签：土地利用论文;  空间分析模型论文;  CA-Markov模型论文;  五常市论文;  黑龙江大学水利电力学院论文;  黑龙江国防科学技术研究院论文;