1、引言

滨海湿地是自然作用和人类活动剧烈地带，受淤涨侵蚀和围垦开发利用的影响，滨海湿地变化频繁且类型多样，及时监测湿地资源及其生态环境的变化具有重要的现实意义[1,2]。遥感变化检测技术通过数理模型分析得到不同时相上地表发生的变化，能够及时对滨海湿地资源及其生态环境进行动态监测，为减少人财物等方面的投入具有重要的现实意义。

对于中低分辨率遥感影像，像元级变化检测方法是常采用的变化检测方法，具有无需人工干预、简单易行等优点，但像元级方法存在对影像配准、辐射校正误差敏感，检测结果椒盐现象严重等问题[3]。随着遥感影像空间分辨率的提高，像元级变化检测方法已无法满足遥感应用的需求[4]。与像元级方法不同，对象级变化检测方法以影像对象为分析单元，具有对影像配准、辐射校正要求低，改善椒盐现象等优点，但对象级方法受影像分割尺度等参数影响较大、实现过程繁琐。可见，像元级方法和对象级方法各有优势又存在不足，因此目前将像元级与对象级变化检测方法各自的优势进行结合，是高分辨率遥感影像变化检测研究的热点[5]。国内外学者已开展了一些研究[6,7,8],Cao等[6]首先采用改进的多尺度水平集方法对差异图像进行分割得到像元级变化检测结果，其次在面向对象分割的基础上，计算对象内变化像元的比例来确定发生变化与未变化的训练样本，最后采用支持向量机进行分类得到最终的变化检测结果。Xiao等[7]首先采用自适应阈值方法对差异影像进行分割得到像元级变化检测结果，其次对像元级结果进行形态学闭、形态学开等操作得到发生变化的对象，最后利用形态学建筑物指数[9]识别建筑区域的变化对象。冯文卿等[8]首先利用水平集方法对光谱差异和纹理差异融合影像进行分割，得到像元级变化检测结果，其次对多时相叠加影像的主成分进行超像素分割，得到不同尺度大小的超像素区域，最后通过对超像素区域内的像元级变化检测结果建立决策级融合规则得到最终的变化检测结果。张志强等[5]首先采用随机森林分类器对影像特征集进行分类得到像元级建筑物变化检测结果，然后对后时相影像进行分割，通过计算分割对象内发生变化的建筑物像元比例识别建筑物变化区。

以上结合像元级与对象级的变化检测方法均是像元级方法和对象级方法串行的方式，即像元级方法获得的结果是对象级方法的前提，其中像元级方法多数是采用分割或分类方法直接对差异影像进行分析处理，这种处理方式在研究区未发生变化的像元所占比重较大时，较多的未变化像元势必影响变化像元的检测过程，且变化检测过程耗时。为了解决该问题，分析差异影像发现，差异影像上变化区域具有一定的显著性，因而可以在像元级变化检测过程中引入视觉显著性检测技术，利用视觉显著性检测算法初步定位出可能发生变化的区域，然后再利用分割或分类算法将初步提取的可能发生变化区域分为变化与未变化区域[10,11]，这样降低了较多的未变化像元对提取变化像元的影响，提高了变化检测结果精度。然而，在像元级变化检测阶段如何顾及高分辨率遥感影像的空间邻域信息来提高初始变化检测效果，仍值得深入研究，另外，在融合像元级和对象级的变化检测阶段，如何自适应地选择训练样本，也有待继续研究。

基于上述分析，本文提出一种显著图引导的结合像元级与对象级变化检测方法。该方法首先提取影像特征，得到特征差异影像；其次利用视觉显著性检测算法获得显著图，采用结合模糊C均值和马尔可夫随机场方法对显著区域进行分割获得初始像元级变化检测结果；最后进行面向对象的影像分割，通过构建对象的不确定性指数自适应选择训练样本，采用随机森林分类器进行分类得到最终变化检测结果。以江苏盐城滨海湿地为典型试验区，利用高分辨率的资源三号（Ziyuan-3,ZY-3影像数据验证该方法的有效性。

2、结合像元级与对象级的变化检测方法

显著图引导的结合像元级与对象级变化检测方法总体技术流程如图1所示。首先，对变化前T1时相和变化后T2时相的影像数据进行预处理，根据试验区湿地场景的特点，提取湿地特征；其次，采用视觉显著性检测算法提取显著区域，这些显著区域在特征差异图上与周围相邻区域对比明显，为可能发生变化的区域，利用结合模糊C均值(Fuzzy C Mean,FCM)和马尔可夫随机场(Markov Random Field,MRF)方法对显著区域进行分割，得到初始的像元级变化检测结果；最后，采用面向对象的分割方法进行影像分割，根据影像对象内像元的变化情况构建不确定性指数自适应选择训练样本，利用随机森林分类器进行分类得到最终的变化检测结果。

图1显著图引导的结合像元级与对象级变化检测方法

2.1特征提取

分析光学影像上湿地的场景特点得出，盐沼植被、农田是与植被相关的地物类型，内陆水体、海域和养殖塘是与水体相关的地物类型，光滩和建筑用地在光学影像上表现出高反射特性，基于此，提取光学影像的亮度、归一化差异植被指数（Normalized Difference Vegetation Index,NDVI），归一化差异水体指数（Normalized Difference Water Index,NDWI)3个特征描述湿地场景。

不同特征的取值范围不同，为了使不同特征在后续处理中具有相同的地位，需进行归一化。假设Fi(i=1,2)为提取的Ti(i=1,2)时刻影像特征，根据式(1)进行归一化，采用差值法对归一化特征Fin进行运算，生成特征差值图DF，表达如式(2)所示。

式中：F2n为归一化处理后的变化后影像特征；F1n为归一化处理后的变化前影像特征。

2.2显著图引导的像元级变化检测

2.2.1视觉显著性检测算法

特征差异图上发生变化的像元值往往大于其周围未发生变化的像元值，那么变化像元由于视觉上的强对比度能引起计算机视觉注意，于是变化检测与视觉显著性检测能联系起来。利用计算机视觉感知对图像进行显著性检测，能够提高图像分割、目标检测的准确性[12]。目前显著性检测算法有：基于图的(Graph-Based Visual Saliency,GBVS)方法[13]，基于上下文感知的方法[14]，最大对称环绕(Maximum Symmetric Surround,MSS)算法[15]等。当显著性区域占研究区整个区域比例超一半以上，或研究区域背景复杂时，MSS算法能够突出目标的显著性，获得高质量的显著图。

本文采用MSS算法对特征差异图进行显著性检测，得到显著图S(xi)。利用Otsu法对显著图进行分割，提取显著性区域，表达如式(3)所示。当S(xi)大于阈值τ时，表明像元xi位于显著区域，即为可能发生变化的区域，显著度SE用1表示；否则像元xi位于非显著区域，即为未发生变化区域，显著度SE用0表示。显著度为0的区域不再参与后续的变化检测过程。

2.2.2 FCM-MRF法的变化检测

MRF模型能够刻画图像的空间信息相关性和上下文关系，在变化检测过程中，不仅利用了当前像元的图像特征信息，而且还考虑了当前像元的空间上下文信息，因而能够很好地抑制噪声、配准误差等对变化检测结果的影响。特别是对于高分辨率影像，利用MRF模型建立高分辨率影像的空间邻域关系，一定程度上能够削弱像元级变化检测结果的椒盐效应。

传统MRF模型利用“0”和“1”表示像元xi与其邻域N(i)的空间关系，1表示中心像元xi的邻域像元xj发生变化，0表示邻域像元xj未发生变化，通过计算邻域像元xj中发生变化和未变化的数量确定当前中心像元xi的变化类型。当邻域像元xi发生变化的个数大于未发生变化的个数时，当前中心像元xi确定为发生变化像元，否则中心像元xi确定为未发生变化像元。由于地表覆盖的复杂性和影像分辨率的限制，图像上通常存在大量的混合像元，此时利用“0”和“1”不再能准确表达空间邻域关系，而且还容易造成空间邻域信息的过度利用。针对该问题，郝明[16]提出利用空间引力模型改进的Potts模型来定义像元xi的空间邻域关系，公式表达如式(4)—式(5)所示。

式中：xj,j∈N (i){j=1,2,⋯,8}表示中心像元xi的3×3邻域窗口内的像元；z(pi)和z(pj)分别表示中心像元xi与邻域像元xj的类别标记；Dij表示中心像元xi与邻域像元xj的欧式距离；pi和pj分别表示中心像元xi和邻域像元xj属于类别k,k∈{未变化,变化}的隶属度，是利用FCM算法计算得到，因而上述改进的MRF方法记作FCM-MRF方法。本文利用FCM-MRF方法仅对显著性检测算法得到的显著区域进行处理，得到初始的变化检测结果。

2.3像元级与对象级结合的变化检测

2.3.1面向对象的影像分割

采用多尺度分割算法进行面向对象的影像分割，多尺度分割算法的参数有尺度、颜色因子和紧致度，其中尺度在颜色因子、紧致度的参与下控制分割结果的异质性，尺度参数的取值直接影响面向对象影像分析的结果，为此本文采用不同分割尺度下的局域方差变化率（Rates of Change of Local variance,ROC-LV）估计最佳分割尺度[17,18]。

2.3.2训练样本自适应选择

训练样本选择是在像元级变化检测结果和面向对象的分割结果基础上开展的，根据对象Ri中变化和未变化像元个数占整个对象Ri像元个数的比值，计算分割对象的不确定性指数T[11]，如式(6)所示。

式中：nc、nu和n分别为对象Ri中像元级变化检测得到的变化像元个数、未变化像元个数和对象总像元个数。设置阈值Tm，利用多数投票法判定对象Ri的属性li，如式(7)所示。

式中：li=1,2,3表示对象Ri的属性为未变化对象、不确定对象和变化对象。利用多数投票法，当分割对象内变化像元的比例大于设置的阈值时，则该分割对象归为变化对象；当分割对象内非变化像元的比例大于设置的阈值时，该分割对象归为未变化对象。否则分割对象为不确定对象。这里将通过阈值Tm确定的变化对象和未变化对象选为训练样本，不确定性阈值Tm取值范围为[0.50,0.95]。当不确定指数T>Tm分割对象被选为变化样本；当不确定指数T<-Tm时分割对象被选为未变化样本；当-Tm≤T≤Tm时分割对象为不确定是否发生变化的对象。之后利用选择的变化和未变化训练样本，采用随机森林分类方法对不确定是否发生变化的分割对象进行分类。

2.4精度验证

为了评价变化检测效果，将变化检测结果与真实变化参考图进行对比，计算受试者工作特征(Receiver Operating Characteristic,ROC)曲线下的面积(Area Under Curve,AUC)、总体精度(Over Accuracy,OA)，精度比(Accuracy Rate,AR)和虚检率(False Detection Rate,FDR)指标值。AUC越大，表明变化检测结果越好，OA指标反映变化检测结果的整体检测精度，AR指标反映变化检测结果的漏检情况，FDR指标反映变化检测结果的虚检情况。OA、AR越大，变化检测结果越好，FDR越小，变化检测结果越好。OA、AR和FDR指标计算公式如下：

式中：Sr为变化检测结果中正确检测为变化的像元个数；Sc为真实变化参考图中变化的像元个数；Sa为变化检测结果中检测为变化的像元个数；Sn为变化检测结果中正确检测为未变化的像元个数；S为像元个数总数。

由本文总体流程图可知，本文方法受到分割尺度和不确定性阈值的影响，因此为了验证本文方法的有效性，首先分析不同分割尺度和不确定性阈值对变化检测结果的影响，然后将本文方法与像元级变化检测方法、对象级变化检测方法以及现有的两种显著图引导的变化检测方法（SG-PCAK方法[10]和SG-RCVA-RF方法[11]）进行对比分析。对比方法的实现过程描述如下：

(1)像元级变化检测方法：该方法为图1中显著图引导的像元级变化检测方法。

(2)对象级变化检测方法：该方法首先采用分割尺度55，颜色因子0.4，紧致度0.5对亮度、NDVI和NDWI 3个特征影像进行分割；然后采用差值法对分割对象的特征均值进行运算得到特征差异图，采用MSS算法对特征差异图进行显著性检测提取出显著性区域；最后采用FCM-MRF方法对显著性区域进行分割得到对象级变化检测结果。

(3)SG-PCAK方法：该方法首先进行变化矢量分析得到变化强度图；然后采用基于上下文感知的方法提取显著性区域，仅保留显著区域的变化强度值，得到新变化强度图；最后采用基于主成分分析和K均值聚类的方法对新变化强度图进行分割，得到变化检测结果。

(4)SG-RCVA-RF方法：该方法首先采用变化矢量分析得到变化强度图，采用GBVS方法提取显著性区域，通过阈值法对显著区域进行分割得到像元级变化检测结果；然后采用分割尺度55，颜色因子0.4，紧致度0.5进行分割得到影像对象，最后利用像元级变化检测结果对影像对象进行分析，得到最终变化检测结果。

3、试验区概况与数据源

3.1试验区概况

以城镇化、工业化背景下的江苏盐城滨海湿地为典型研究区。江苏盐城滨海湿地位于江苏省盐城市的沿海地带，地跨响水县、滨海县、射阳县、大丰市和东台市五县/市，位于119°40′E—121°E,32°40′N—34°30′N，为我国典型的淤泥质滩涂湿地，其海岸线长约582 km，约占整个江苏海岸线长度的56%，北窄南宽且呈条带状分布，射阳河口以北为侵蚀型海岸，斗龙港口以南为淤涨型海岸，射阳河口至斗龙港之间为过渡型海岸。研究区地势平坦，绝大部分地区高程小于5 m，地跨暖温带和北亚热带，年平均气温为13.7～14.6℃，光照充足，雨水丰沛，无霜期210～224 d，年平均降水量1000 mm左右。2002年盐城湿地被列入《国际重要湿地名录》[19]，其包括两大国际重要湿地，分别为国家级珍禽自然保护区和大丰麋鹿国家级自然保护区。本文以盐城湿地的东南部——新开港区域为试验区进行实验，影像区域如图2所示。

图2盐城滨海湿地ZY-3卫星影像

3.2数据源

由图2可知，试验区包括海域、农田、养殖塘、盐沼植被、内陆水体、光滩、建筑用地等地物，湿地类型复杂，为了精细地提取试验区湿地的变化，在自然资源卫星影像云服务平台[20]获取了2013年和2018年ZY-3卫星数据，影像的基本信息如表1所示。比照2013年和2018年影像可以看出，试验区发生的变化是：海域转为盐沼植被（即盐沼植被向海域扩散，如图2区域A所示）、海域转为养殖塘（即海域围垦为养殖塘，如图2区域B所示）、农田转为养殖塘（如图2区域C所示）、盐沼植被转为农田（如图2区域D所示）和盐沼植被转为养殖塘（如图2区域E所示）。

表1资源三号卫星影像基本信息

对ZY-3数据进行预处理，首先是图像配准，在2景影像上分别选取一定数量的同名点进行自动配准，配准误差控制在1个像元之内，该步骤在ENVI软件下完成。其次是图像裁剪，将2景影像进行叠加，裁剪公共区域，生成试验区影像。最后利用ENVI FLAASH工具对光学影像进行大气校正。

4、 结果及分析

4.1分割尺度的影响

不同分割尺度对变化检测结果的影响分析如图3所示。图3(a)中尺度55，尺度180，尺度360，尺度475为ROC-LV方法估计的分割尺度，尺度55的AUC值最大，为0.86174，由图3(b)可知，尺度55的总体精度最高，精度比最高，虚检率最小。由此可见，尺度55的变化检测结果最好，其分割结果基本与湿地地块边界吻合，能够将发生变化的地块与未发生变化的地块各自分割成独立的地块，因此在结合像元级和对象级的变化分析过程中，变化地块与未变化地块能够分开分析。

4.2不确定性阈值的影响

采用尺度55进行面向对象的影像分割，得到影像对象，根据式(6)和式(7)自适应选择训练样本，不确定性阈值Tm对变化检测结果的影响分析如图4所示。由图4可知，Tm为0.7时，AUC、总体精度和精度比均最大。Tm取值逐渐增大时，AUC值先下降后上升，在Tm为0.7时达到最大值后又逐渐下降。Tm逐渐增大时，变化检测总体精度保持平稳，在Tm=0.7达到最大值93.51%，随后逐渐下降。精度比曲线在[0.5,0.6]范围内下降，在[0.6,0.7]范围内上升，在[0.7,0.95]范围内逐渐下降。由图4(a)和图4(b)分析可知，阈值Tm不宜取过低或过高，取值低时，将变化像元比重小的对象归为变化对象，造成虚检率过高；取值高时，将变化像元比重高的对象归为变化对象，造成漏检率过高，Tm为0.7时变化检测结果最优。

图3分割尺度对变化检测结果的影响

图4不确定性阈值Tm对变化检测结果的影响

4.3不同变化检测方法的对比

不同变化检测方法的变化检测结果如图5所示，图5(f)为真实变化参考图，黑色表示未变化区域，白色表示变化区域，在对研究区进行野外调研的基础上，参照高分辨率影像资料进行人工勾绘得到。

图5不同变化检测方法的变化检测结果

图5(a)和(b)分别为像元级变化检测方法和对象级变化检测方法得到的变化检测结果，像元级方法由于考虑了像元空间邻域关系对像元进行分割，所以与SG-PCAK方法相比海域的虚检像元少，然而像元级方法在农田区域将未变化的农田误检为变化，所以仅利用像元级方法提取湿地变化区域，会产生较多的虚检。对象级方法相较于像元级方法，虚检现象有所降低，椒盐效应有所减弱，但漏检现象比像元级方法有所明显。原因是在面向对象的分割过程中，若将发生变化的像元分割至未变化像元个数比重较大的对象中时，该对象被归为未变化对象，则对象内那部分发生变化的像元划归为未变化，导致了漏检现象。

图5(c)和(d)分别为SG-PCAK方法和SG-RC-VA-RF方法得到的变化检测结果。SG-PCAK方法直接利用光谱波段进行变化检测，造成盐沼植被转变为农田的变化（图2区域D）未能在差异影像上表现出来，导致SG-PCAK方法的变化检测效果欠佳，同时海域虚检现象明显，将由海浪波动引起的伪变化虚检为变化。SG-RCVA-RF方法相较于像元级、对象级和SG-PCAK方法，检测效果有明显提升，虚检现象降低，然而该方法利用光谱波段进行变化检测，并未考虑湿地的场景特性，表现为海域围垦开发为养殖塘的变化（图2区域B）未能检测出来。图5(e)为提出方法的变化检测结果，相较上述4种方法，最接近真实变化参考图，将海域围垦为养殖塘（图2区域B）、盐沼植被开发为农田（图2区域D）等变化检测出来，检测效果最优。

不同变化检测方法的定量对比结果如图6。由图6(a)可知，提出方法的AUC值（0.861 74）最大，像元级方法的AUC值（0.811 23）次之，SG-PCAK方法的AUC值（0.752 66）最低。由图6(b)可知，提出方法的变化检测总体精度、精度比最高，总体精度为93.51%，高于SG-PCAK方法5.95%，同时虚检率最低，与SG-PCAK方法相比，虚检率降低了35.96%，精度比提高了29.24%。与SG-RCVA-RF方法相比，精度比提高了14.23%，与像元级方法和对象级方法相比，虚检率分别降低了29.04%和22.78%。AUC、OA、AR和FDR 4个指标分析结果均表明，提出方法的湿地变化检测结果优于像元级、对象级、SG-PCAK和SG-RCVA-RF方法。

5、结论与讨论

5.1结论

本文以江苏盐城滨海湿地为典型试验区，针对像元级变化检测对配准误差敏感，椒盐现象严重；而对象级变化检测受分割参数影响大，过程繁琐问题，提出了显著图引导的结合像元级与对象级变化检测方法。实验探讨了分割尺度和不确定性阈值对该方法变化检测结果的影响，通过与像元级、对象级变化检测方法相比，虚检率分别降低了29.04%和22.78%，通过与现有显著图引导的变化检测方法SG-PCAK、SG-RCVA-RF相比，精度比分别提高了29.24%和14.23%，验证了该方法的有效性。

本文的主要贡献体现在：

(1)在滨海湿地的变化检测应用中，充分考虑了滨海湿地的场景特点，提取了能够有效表征滨海湿地植被、水体等变化的特征，实现了在输入特征层面湿地变化表达的完整性。

图6不同变化检测方法的对比

(2)在像元级变化检测过程中，引入了视觉显著性检测方法，采用“先初步定位可能变化区，再精确提取变化”的思路，有效提高了变化检测的运算效率，另外，在精确提取变化像元的过程中，顾及了高分辨率影像的空间邻域关系，缓解了像元级变化检测结果的椒盐现象。

(3)在结合像元级与对象级的变化检测过程中，将像元级方法与对象级方法有效串行，全自动化、自适应地生成训练样本，避免了人工选择训练样本的费事费力过程。

5.2讨论

为了提高变化检测精度，国内外学者采用串行或并行方式将像元级方法和对象级方法进行集成[8]。并行方式是指像元级变化检测和对象级变化检测同时进行，然后将两种方法得到的结果进行综合，串行方式是一种方法是另一种方法的前提，最终目的都是得到更好的变化检测结果。其中串行方式的结合像元级与对象级方法总体思路是在得到初始像元级变化检测结果的基础上，通过分析分割对象内像元的变化情况来判断对象的变化属性。文献[7]针对建筑物的变化检测，通过形态学建筑物指数判断变化区域是否是建筑物的变化。文献[8]通过多值逻辑函数判断超像素分割对象的变化属性。文献[21]针对洪水的监测，在像元级变化检测结果的基础上，采用模糊聚类和最近距离聚类方法得到最终的洪水区域。

本文以滨海湿地变化的检测为研究对象，按照串行方式将像元级与对象级方法集成，在像元级变化检测阶段和集成像元级与对象级方法阶段提出新方法，在像元级变化检测阶段，引入显著性检测算法提取特征差异图上的显著区域，用于快速获得初始变化检测结果，在集成像元级与对象级方法阶段，构建不确定指数自适应获取用于分类的训练样本，通过采用随机分类器进行分类得到最终的变化检测结果。本文研究成果在理论层面解决了目前结合像元级与对象级变化检测方法仍存在的问题，在应用层面提高了湿地变化监测的精度，整个技术流程实现了全自动化，无需人工设置阈值和选择训练样本，降低了变化检测过程中的主观因素影响，而且可操作性强，本文方法可以推广至多种数据源的其它湿地类型以及国土遥感监测等领域中。但是需要说明的是：

(1)本文在对湿地的变化进行检测时，针对湿地场景的特点，考虑了湿地类型是与植被、水体相关的，提取了3个湿地特征。提取的特征不同，也会影响变化检测效果，因而为了实现更高精度的变化检测，特征提取环节需根据研究区域的特点，挖掘合适的特征进行分析。

(2)在面向对象的影像分割环节，尺度参数直接影响生成的对象，本文利用ROC-LV方法估计理想的分割尺度，未来可以考虑利用新颖的最优尺度参数估计方法确定面向对象的分割尺度，以提高变化检测精度。

参考文献:

[1]宫宁,牛振国,齐伟,等.中国湿地变化的驱动力分析[j].遥感学报,2016,20(2):172-183.

[2]左平,李云,赵书河,等. 1976年以来江苏盐城滨海湿地景观变化及驱动力分析[j].海洋学报(中文版),2012,34(1):101-108.

[3]王成军,毛政元,徐伟铭,等.超像素与主动学习相结合的遥感影像变化检测方法[j].地球信息科学学报,2018,20(2):235-245.

[4]李杨,江南,侍昊,等.lANDSAT-8影像的lda模型变化检测[j].地球信息科学学报,2015,17(3):353-360.

[5]张志强,张新长,辛秦川,等.结合像元级和目标级的高分辨率遥感影像建筑物变化检测[j].测绘学报,2018,47(1):102-112.

[8]冯文卿,眭海刚,涂继辉,等.联合像元级和对象级分析的遥感影像变化检测[j].测绘学报,2017,46(9):1147-1155,1164.

[12]何鹏飞.基于不确定性分析的遥感影像变化检测方法研究[d].徐州:中国矿业大学,2017.

[16]郝明.基于空间信息准确性增强的遥感影像变化检测方法研究[d].徐州:中国矿业大学,2015.

[19]郝敬锋,刘红玉,李玉凤,等.基于转移矩阵模型的江苏海滨湿地资源时空演变特征及驱动机制分析[j].自然资源学报,2010,25(11):1918-1929.

[20]自然资源部国土卫星遥感应用中心.自然资源卫星影像云服务平台[db/ol]

[21]冷英,李宁.一种改进的变化检测方法及其在洪水监测中的应用[j].雷达学报,2017,6(2):204-212.

吴瑞娟,何秀凤,王静.结合像元级与对象级的滨海湿地变化检测方法[j].地球信息科学学报,2020,22(10):2078-2087.

基金:国家自然科学基金项目(41830110､41871203);内江师范学院科研资助项目(2019yz02).