剖线强度梯度变化的作物软硬变化区精确划分*
2018-09-28张锦水崔有祯
朱 爽,张锦水,崔有祯
(1.北京工业职业技术学院,北京 100042; 2.北京师范大学遥感科学与工程研究院地理科学学部,北京 100875; 3.北京师范大学地理科学学部,北京 100875)
0 引言
我国是农业大国,及时、准确地获取农作物播种面积信息,对于制定国家/区域农业经济发展规划、指导种植业结构调整,提高农业生产管理水平具有重要的意义[1]。遥感技术具有覆盖范围广、探测周期短的特点,是进行农作物种植面积提取的重要技术手段[2-3]。利用多期遥感影像进行作物检测识别,能够根据作物短期内的光谱变化,定量刻画出作物的生长物候特征,以此为基础进行作物识别,有效消除作物光谱相混的问题,提高作物的识别精度[4]。在众多遥感变化检测方法中,软硬变化检测是一种新型的农作物识别方法,该方法是依据遥感影像作为格网模型产生软硬变化区共存现象提出的[5-6]。
软硬变化检测识别实施的关键是有效地识别出软变化区(Soft change region,SCR)、硬变化区(Hard change region,HCR)[7-9]。当前常用的方法是利用像元变化强度计算生成变化强度图,通过阈值来确定地表的变化、未变化区域(Non change region,NCR)[10-11]。显然,阈值的设定是关键。目前,设定阈值的方法多采用人工目视判读[12-13],但该方法受到人为主观因素的影响比较大,且操作困难,难以形成统一的标准进行推广应用。自动分类法[14-15]是选择样本进行直接分类,如SVM分类方法[16-17]或自动聚类方法,该类方法速度快、效率比较高,在研究中应用较多。此外,还有一种阈值确定方法介于两者之间,即交互式确定阈值,如陈晋等[18]提出双窗口变步长阈值设定方法确定变化区域,从变化强度影像上分别选出变化、未变化的样本,通过迭代设定阈值来不断地逼近最高识别精度,确定全局最优的划分阈值。双窗口变步长阈值确定方法只要选择出合适的变化样本,就能够快速确定阈值,且对变化对象的标识精度较高。但是,上述方法不适合中分辨率遥感影像的SCR、HCR与NCR的划分,原因在于该划分是3类地物特征的划分,而传统阈值划分多用于2类地物分割。
文章提出了一种基于变化强度图的剖线梯度变化方法(Profile based Gradient Change Magnitude,PGCM)来确定阈值,通过从提取的作物图斑内部向外部绘制剖线,计算剖线变化强度的梯度,根据梯度变化的突变特性确定阈值进行HCR、SCR和NCR划分,为作物的软硬变化检测提取提供基础。
1 研究区与数据
研究区位于北京市大兴区,覆盖范围为10km×10km(图1a),属永定河冲积平原。该地区农作物种类较多,主要包括小麦、玉米、花生、大豆等,作物种植地块破碎是这一区域典型的农业景观特征,耕地地块尺寸在0.03~8.53hm2范围之内, 0.67hm2以下的耕地地块占90%以上,在遥感影像上表现为典型的纯净、混合像元共存现象。因此,该区域适合开展PGCM方法进行冬小麦软硬变化区划分的实验研究。
该文采用QuickBird多光谱影像,分辨率为2.4m,包含4个波段(蓝光: 0.45~0.52μm,绿光: 0.52~0.60μm,红光: 0.63~0.69μm,近红外: 0.76~0.90μm),获取日期为2006年4月23日,处于冬小麦物候的拔节期(图1b)。采用文献[19]所使用的方法,以拔节期QB影像为基础,对播种期影像(2005年10月,图1c)进行模拟。在构建模拟影像过程中,充分考虑到各种地物的物候特征和遥感影像的成像特点,将拔节期的休耕地地块按照10%比例模拟为植被; 将冬小麦光谱模拟为裸地光谱信息; 对于在拔节期与冬小麦相混的林地、草地等自然植被,考虑到不同地物在时相上地物反射率的变化,需要在保留各自地物光谱的基础上,随机混入10%的其他类型地物光谱,作为播种期的地物光谱(S±10%S,S为林地、草地的光谱信息)。播种期模拟影像和拔节期真实影像的结合,可有效表达出各种地物的物候特征变化,尤其是准确地刻画出冬小麦从裸地到植被这一独特的光谱特征变化,真实表达冬小麦的生长情况。此外,模拟数据与真实数据之间不存在配准误差的问题,消除了这些因素对变化检测影响。因此,该套研究数据能够有效支撑软硬变化检测的软硬变化区划分的研究。
(a)研究区 (b)4月23日QB影像 (c)10月QB模拟影像图1 研究区及数据
2 基于剖线梯度变化方法(PGCM)
对拔节期影像(图1b)与播种期影像(图1c)逐波段进行差值计算,得到差值影像图。图2是差值影像的一个子区图(以30 m分辨率为例),能够有效地表达冬小麦的空间分布。图2a中红色光谱代表一种作物,作为要提取的目标对象。图2b为对应的变化强度图,高亮的区域变化程度剧烈,为HCR; 黑色区域为NCR; 灰色过渡区域为SCR,直线AB为从小麦地块内部到外部的剖线A→B。图2c为剖线A→B上10个像元的变换强度,显示出冬小麦的变化强度从地块内部向外呈现逐步降低的趋势,这与冬小麦地块内部的HCR向SCR逐步过渡的趋势相吻合(图2c)。
(a) 差值影像(组合RGB=红光、 (b)变化强度图 (c)剖线上像元变化强度 绿光、蓝光差值波段)图2 差值影像子图
图3定量地表达出从地块内部向外部过渡的过程中,整体的变化强度呈逐步降低的趋势(图3a)和变化强度的梯度变化(即剖线上的下一个像元变化强度与上一个像元变化强度之差,图3b)。从剖线上来看(图2c和图3),硬变化区域中多为纯净冬小麦像元,变化强度处于高值,因此在这一范围内变化强度的差值变化不大(1~3个点),但随着距离的增加,由纯净冬小麦像元进入混合区域的冬小麦像元,第3点至第4点的差值幅度突然增大,且为负值,记为T1,该位置可定义为SCR像元的上限; 随着距离的进一步增加,变化幅度继续降低,在第6个像元的时候,差值幅度又一次突然增大(即第5个点与第6个点之间的差值),这正好位于从软变化区向未变化区过渡的状态,记为T2,该位置可定义为SCR像元的下限。
(a)剖线上的变化强度 (b)剖线上变化强度的梯度变化图3 剖线上变化强度和梯度变化
根据上述分析,在此提出基于剖线梯度变化方法(PGCM)进行硬、软和未变化区域划分方法。具体思路见式(1)(2):
ΔCG=CMpj-CMpi,i=1, 2…,n-1,j=i+1
(1)
(2)
其中,CMpi、CMpj分别代表了剖线上的第i个点与下一个点j(j=i+1)之间的变化强度值;BVi(T1)与BVi(T2)分别代表像元i在T1和T2时期k波段的上的光谱值,波段k取值为[1, 2, 3,…,n];n代表总的波段数;ΔCG代表i,j相邻两点之间变化强度的差值(图3b)。
从图3b可以分析出,从HCR→SCR、SCR→NCR两个阶段的ΔCG均发生了突变。从农业景观特征角度分析,耕地多以规整地块的方式排布,农田周边以休耕地为主。对于从作物地块跨越到周边其他地物区域,遥感影像上产生混合像元,这些混合像元内部一般含有一定比例的作物和其他地物。例如,从作物地块(像元内作物种植面积为100%)向休耕地过渡,休耕地像元内部可能含有70%的作物和30%的其它地物,这一突变是能够通过ΔCG表现出来的。为了客观、快速、准确地确定T1、T2,在此提出一个假设:T1定义为当剖线上第1个出现差值下降且幅度最大时所对应的变化强度取值,T2定义为第2个出现最大差值下降且幅度最大时所对应的像元变化强度值。其中,p1代表剖线上第1个落差最大的点对应的像元位置,p2代表剖线上第2个落差最大的点对应的像元位置。由剖线A→B生成的SCR混合像元为p1→p2-1对应的像元集合Si。如果有n条剖线,则会生成n个SCR(n个Si)像元集合,定义为SCRP,即:
图4 试验流程
T1=max(SCRP)
(3)
T2=min(SCRP)
(4)
通过变化强度来确定3个区域(I):HCR、SCR和NCR,见式(5)。
(5)
由式4可以看出,基于PGCM方法来划分HCR、SCR和NCR,3个区域的关键步骤是确定T1、T2。其中,HCR、SCR是进一步进行作物识别的基础。
考虑到作物地块内部仍然存在光谱的不确定性,在提取出SCR后,采用3×3窗口滤波,当窗口内有1/3以上的SCR像元时,保留中心像元,定义为SCR像元; 否则将该像元归属到窗口内占主导的HCR或NCR。如果HCR和NCR个数相同,则先归属为HCR,再进一步通过变化方向确定该像元的类别归属。
3 实验设计与结果分析
该研究针对单维度多时相构成的变化强度图来绘制空间剖线,计算软硬变化区变化强度梯度落差的阈值; 基于阈值来确定软硬变化区域; 针对软变化区进行窗口滤波,消除环境、传感器等因素造成的“椒盐”现象,具体实验流程见图4。
在影像图上选择34个耕地地块,从地块中心出发向外绘制剖线,并通过PGCM方法自动确定出混合像元。由于该次实验是针对冬小麦进行提取,只需要在冬小麦地块上选取剖线,其他地物不予考虑。剖线要从典型的耕地地块内部向外部绘制,穿越耕地地块后进入由混合像元构成的过渡带,再进入非耕地。由于作物地块种植一般为东西和南北方向,因此剖线方向也要求为东西、南北方向,以避免任意方向剖线穿越地块带来的误差。当地块尺寸比较大时,剖线在地块内、外需要各自跨越至少5个像元。之后将剖线输入到相对应的变化强度影像中,进行 HCR、SCR、NCR划分(图5)。
表1 不同分辨率下变化强度阈值
尺度(m)T1T2517353101705020168463016545401634350160416015839
为验证本方法的有效性,将2.4m分辨率影像下选择出的地块样本聚合到5m、10m、20m、30m、40m、50m和60m 分辨率,分别得到阈值集合,进行不同分辨率下作物软硬区的划分。图6表明了不同分辨率下研究区子区变化强度图和其对应的HCR、SCR和NCR划分结果。从图6可以看出,在不同分辨率尺度下,PGCM均能识别出SCR像元,变化强度阈值见表1。在较高分辨率尺度时,SCR像元主要分布在两个位置:地块边缘和地块内部。例如, 5m分辨率的模拟影像,SCR像元多产生在冬小麦地块的边缘; 另外,当影像分辨率较高时,尽管冬小麦像元多以纯净像元的形式存在,但长势、环境因素或种植情况的不同会导致冬小麦像元之间存在一定的光谱差异,因此还有一部分SCR像元存在于冬小麦地块内部; 除这两种情况之外,部分冬小麦像元的变化强度比较低,与地块边界混合像元的变化强度接近,也被划为SCR像元。已有研究证明,当图像分辨率增高时,地物光谱的异质性会导致地物识别精度的降低[20]。
随着影像分辨率尺度逐步增大,地块内部的光谱异质性逐步降低,光谱异质性对于HCR、SCR之间的划分影响减小。但是,由于研究区作物地块面积较小,地块破碎,随着分辨率的逐步降低,地块中的纯净冬小麦像元数量迅速减少,混合像元现象严重,到了30m分辨率时,大部分冬小麦地块为混合像元或长势不好的冬小麦像元,而HCR冬小麦像元则很少。从识别结果来看,在任何尺度上,PGCM均能检测出相应的SCR像元,针对这些像元如果采用硬变化检测方法进行作物类型提取将会产生误差,因而需要使用软变化检测方法。
(a)研究区数字化解译结果 (b)研究区子区数字化解译结果图5 研究区及子区数字化解译结果注:以QB数据为底图,结合野外数据,解译数字化获得冬小麦真值(黄色框表示); 蓝色框为研究区子区
图6 研究区子区HCR、SCR和NCR划分结果(白色为HCR,红色为SCR,黑色为NCR)
以人工数字化冬小麦为基础,采用聚合方式将整个区域划分为HCR、SCR和NCR 3个区域。表2中的HCR_真值、SCR_真值和NCR_真值代表尺度下的3个区域内的像元数,SCR1 _ PGCM、SCR2 _ PGCM、SCR3_ PGCM分别对应HCR、SCR和NCR 3个区域内PGCM方法提取出的SCR像元个数。可以看出,真值获取的3个区域均含有PGCM方法提取出的SCR像元,但HCR、NCR区域混入的比例相对较低。其中,HCR区混入比例为11%~16%,混入比例随着分辨率下降(5~60m)而降低; NCR区混入比例更低,为3%~4%,且受影像分辨率尺度的影响不大。对于SCR的真值区域,SCR_ PGCM所占比例高达74%~86%,表明PGCM识别的SCR像元的精度较高。整体而言,PGCM对于HCR、SCR和NCR,3个区域SCR像元的提取精度较高,对于在HCR_真值、NCR_真值混入的SCR像元,将进一步通过限制性端元SVM(Constrained Endmember SVM,CESVM)[21]的方式确定端元类型和数量来进行混合像元分解,以消除SCR_PGCM识别误差带来的问题。
表2 PCGM识别出的SCR像元个数
尺度(m)HCR_真值SCR1_PGCMSCR_真值SCR2_PGCMNCR_真值SCR3_PGCM5747 820120 476331 624246 826984 523119 22210150 17822 303157 422119 188707 38424 6472023 5402 97169 86856 163160 5925 539306 49977941 15034 95465 2382 010402 35526727 53624 25033 6091 5515098512019 67416 33019 941906604385014 79412 71912 991552
4 结论与讨论
硬变化区/软变化区的划分是进行软、硬变化检测的基础。该文针对检测过程中变化强度阈值划分的问题,提出了一种基于剖线梯度变化方法(Profile based Gradient Change Magnitude,PGCM)来确定阈值,进行HCR、SCR和NCR 3个区域的划分,以提高软硬变化检测方法的效率和准确性。
从识别结果来看,PGCM方法能够在地块边界探测到软变化像元,根据阈值进一步划分出HCR、SCR像元。在多个尺度下(像元分辨率为5~60 m)HCR区识别误差控制在11%~16%,NCR区混入比例为3%~4%,SCR区识别比例74%~86%,识别精度较高,为作物SHCD变化检测打下了基础,识别结果与冬小麦空间分布结果保持一致性。
研究仍有一些方面有待研究:目前该研究使用的方法无法解决转换方向的问题,会导致其他的变化类型混入到作物SCR之中,这一点需要通过在CESVM变化检测中进一步的确定,排除其他地物的影响。另一方面,同期作物光谱与目标作物相似,对划分结果造成了一定的影响,需要进一步深入研究。