CN114018833B - 基于高光谱遥感技术估算土壤重金属含量的方法 - Google Patents
基于高光谱遥感技术估算土壤重金属含量的方法 Download PDFInfo
- Publication number
- CN114018833B CN114018833B CN202111309905.3A CN202111309905A CN114018833B CN 114018833 B CN114018833 B CN 114018833B CN 202111309905 A CN202111309905 A CN 202111309905A CN 114018833 B CN114018833 B CN 114018833B
- Authority
- CN
- China
- Prior art keywords
- spectrum
- soil
- heavy metal
- sample
- value
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 239000002689 soil Substances 0.000 title claims abstract description 270
- 238000000034 method Methods 0.000 title claims abstract description 168
- 229910001385 heavy metal Inorganic materials 0.000 title claims abstract description 127
- 238000005516 engineering process Methods 0.000 title claims abstract description 17
- 238000001228 spectrum Methods 0.000 claims abstract description 251
- 230000003595 spectral effect Effects 0.000 claims abstract description 69
- 238000002310 reflectometry Methods 0.000 claims abstract description 55
- 229910017052 cobalt Inorganic materials 0.000 claims abstract description 48
- 239000010941 cobalt Substances 0.000 claims abstract description 48
- GUTLYIVDDKVIGB-UHFFFAOYSA-N cobalt atom Chemical compound [Co] GUTLYIVDDKVIGB-UHFFFAOYSA-N 0.000 claims abstract description 48
- 238000004458 analytical method Methods 0.000 claims abstract description 41
- 238000001914 filtration Methods 0.000 claims abstract description 33
- 238000012795 verification Methods 0.000 claims abstract description 31
- 238000012952 Resampling Methods 0.000 claims abstract description 23
- 229910052785 arsenic Inorganic materials 0.000 claims abstract description 14
- RQNWIZPPADIBDY-UHFFFAOYSA-N arsenic atom Chemical compound [As] RQNWIZPPADIBDY-UHFFFAOYSA-N 0.000 claims abstract description 14
- VYZAMTAEIAYCRO-UHFFFAOYSA-N Chromium Chemical compound [Cr] VYZAMTAEIAYCRO-UHFFFAOYSA-N 0.000 claims abstract description 13
- 229910052804 chromium Inorganic materials 0.000 claims abstract description 13
- 239000011651 chromium Substances 0.000 claims abstract description 13
- RYGMFSIKBFXOCR-UHFFFAOYSA-N Copper Chemical compound [Cu] RYGMFSIKBFXOCR-UHFFFAOYSA-N 0.000 claims abstract description 12
- 229910052802 copper Inorganic materials 0.000 claims abstract description 12
- 239000010949 copper Substances 0.000 claims abstract description 12
- 238000011426 transformation method Methods 0.000 claims abstract description 12
- 238000011156 evaluation Methods 0.000 claims abstract description 8
- 238000000227 grinding Methods 0.000 claims abstract description 8
- 238000010219 correlation analysis Methods 0.000 claims abstract description 7
- 238000007605 air drying Methods 0.000 claims abstract description 6
- 239000000523 sample Substances 0.000 claims description 159
- 230000009466 transformation Effects 0.000 claims description 43
- 238000005259 measurement Methods 0.000 claims description 35
- 238000012937 correction Methods 0.000 claims description 26
- 230000008569 process Effects 0.000 claims description 25
- 238000005070 sampling Methods 0.000 claims description 25
- 238000004364 calculation method Methods 0.000 claims description 18
- 230000000694 effects Effects 0.000 claims description 18
- 238000000611 regression analysis Methods 0.000 claims description 17
- PXHVJJICTQNCMI-UHFFFAOYSA-N Nickel Chemical compound [Ni] PXHVJJICTQNCMI-UHFFFAOYSA-N 0.000 claims description 14
- 238000013519 translation Methods 0.000 claims description 14
- 238000012360 testing method Methods 0.000 claims description 13
- 230000001419 dependent effect Effects 0.000 claims description 12
- 238000009499 grossing Methods 0.000 claims description 12
- 239000004810 polytetrafluoroethylene Substances 0.000 claims description 12
- 229920001343 polytetrafluoroethylene Polymers 0.000 claims description 12
- 238000010521 absorption reaction Methods 0.000 claims description 11
- 238000007789 sealing Methods 0.000 claims description 11
- 238000009826 distribution Methods 0.000 claims description 9
- 238000013210 evaluation model Methods 0.000 claims description 9
- HCHKCACWOHOZIP-UHFFFAOYSA-N Zinc Chemical compound [Zn] HCHKCACWOHOZIP-UHFFFAOYSA-N 0.000 claims description 7
- 229910052793 cadmium Inorganic materials 0.000 claims description 7
- BDOSMKKIYDKNTQ-UHFFFAOYSA-N cadmium atom Chemical compound [Cd] BDOSMKKIYDKNTQ-UHFFFAOYSA-N 0.000 claims description 7
- 229910052759 nickel Inorganic materials 0.000 claims description 7
- 238000002360 preparation method Methods 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 7
- 229910052725 zinc Inorganic materials 0.000 claims description 7
- 239000011701 zinc Substances 0.000 claims description 7
- KRHYYFGTRYWZRS-UHFFFAOYSA-N Fluorane Chemical compound F KRHYYFGTRYWZRS-UHFFFAOYSA-N 0.000 claims description 6
- GRYLNZFGIOXLOG-UHFFFAOYSA-N Nitric acid Chemical compound O[N+]([O-])=O GRYLNZFGIOXLOG-UHFFFAOYSA-N 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 6
- 238000001816 cooling Methods 0.000 claims description 6
- 238000005260 corrosion Methods 0.000 claims description 6
- 238000012067 mathematical method Methods 0.000 claims description 6
- 229910017604 nitric acid Inorganic materials 0.000 claims description 6
- VLTRZXGMWDSKGL-UHFFFAOYSA-N perchloric acid Chemical compound OCl(=O)(=O)=O VLTRZXGMWDSKGL-UHFFFAOYSA-N 0.000 claims description 6
- 239000012498 ultrapure water Substances 0.000 claims description 6
- 229910052720 vanadium Inorganic materials 0.000 claims description 6
- LEONUFNNVUYDNQ-UHFFFAOYSA-N vanadium atom Chemical compound [V] LEONUFNNVUYDNQ-UHFFFAOYSA-N 0.000 claims description 6
- 238000012935 Averaging Methods 0.000 claims description 4
- 241000196324 Embryophyta Species 0.000 claims description 4
- 238000002790 cross-validation Methods 0.000 claims description 4
- 238000004088 simulation Methods 0.000 claims description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 4
- 239000004695 Polyether sulfone Substances 0.000 claims description 3
- 206010036790 Productive cough Diseases 0.000 claims description 3
- VYPSYNLAJGMNEJ-UHFFFAOYSA-N Silicium dioxide Chemical compound O=[Si]=O VYPSYNLAJGMNEJ-UHFFFAOYSA-N 0.000 claims description 3
- 230000002159 abnormal effect Effects 0.000 claims description 3
- 239000011248 coating agent Substances 0.000 claims description 3
- 238000000576 coating method Methods 0.000 claims description 3
- 239000000428 dust Substances 0.000 claims description 3
- 230000002708 enhancing effect Effects 0.000 claims description 3
- 229910052736 halogen Inorganic materials 0.000 claims description 3
- 150000002367 halogens Chemical class 0.000 claims description 3
- 238000009616 inductively coupled plasma Methods 0.000 claims description 3
- 238000012886 linear function Methods 0.000 claims description 3
- 238000012417 linear regression Methods 0.000 claims description 3
- 239000000463 material Substances 0.000 claims description 3
- 239000004033 plastic Substances 0.000 claims description 3
- 229920006393 polyether sulfone Polymers 0.000 claims description 3
- 238000002203 pretreatment Methods 0.000 claims description 3
- 230000001681 protective effect Effects 0.000 claims description 3
- 230000009467 reduction Effects 0.000 claims description 3
- 210000003802 sputum Anatomy 0.000 claims description 3
- 208000024794 sputum Diseases 0.000 claims description 3
- 238000010025 steaming Methods 0.000 claims description 3
- 239000004575 stone Substances 0.000 claims description 3
- 238000012549 training Methods 0.000 claims description 3
- 241000361919 Metaphire sieboldi Species 0.000 claims description 2
- 230000009286 beneficial effect Effects 0.000 claims description 2
- 238000004422 calculation algorithm Methods 0.000 claims description 2
- 238000007865 diluting Methods 0.000 claims description 2
- -1 removing gravel Substances 0.000 claims description 2
- 238000007873 sieving Methods 0.000 claims description 2
- 238000012546 transfer Methods 0.000 claims description 2
- 238000005303 weighing Methods 0.000 claims description 2
- 238000011160 research Methods 0.000 description 26
- 230000000875 corresponding effect Effects 0.000 description 20
- 241000218378 Magnolia Species 0.000 description 10
- 238000010238 partial least squares regression Methods 0.000 description 7
- 238000006243 chemical reaction Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 238000012628 principal component regression Methods 0.000 description 6
- 239000000126 substance Substances 0.000 description 6
- 241000282414 Homo sapiens Species 0.000 description 5
- 230000008901 benefit Effects 0.000 description 5
- 238000012544 monitoring process Methods 0.000 description 5
- 238000000985 reflectance spectrum Methods 0.000 description 4
- 230000007547 defect Effects 0.000 description 3
- 229910052500 inorganic mineral Inorganic materials 0.000 description 3
- 239000011707 mineral Substances 0.000 description 3
- XEEYBQQBJWHFJM-UHFFFAOYSA-N Iron Chemical compound [Fe] XEEYBQQBJWHFJM-UHFFFAOYSA-N 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 239000003344 environmental pollutant Substances 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 231100000719 pollutant Toxicity 0.000 description 2
- 239000011148 porous material Substances 0.000 description 2
- 238000012847 principal component analysis method Methods 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 238000000844 transformation Methods 0.000 description 2
- BVKZGUZCCUSVTD-UHFFFAOYSA-L Carbonate Chemical compound [O-]C([O-])=O BVKZGUZCCUSVTD-UHFFFAOYSA-L 0.000 description 1
- 241001465754 Metazoa Species 0.000 description 1
- BPQQTUXANYXVAA-UHFFFAOYSA-N Orthosilicate Chemical compound [O-][Si]([O-])([O-])[O-] BPQQTUXANYXVAA-UHFFFAOYSA-N 0.000 description 1
- QAOWNCQODCNURD-UHFFFAOYSA-L Sulfate Chemical compound [O-]S([O-])(=O)=O QAOWNCQODCNURD-UHFFFAOYSA-L 0.000 description 1
- 230000002776 aggregation Effects 0.000 description 1
- 238000004220 aggregation Methods 0.000 description 1
- 238000012271 agricultural production Methods 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 239000004927 clay Substances 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 241001233061 earthworms Species 0.000 description 1
- 230000005274 electronic transitions Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000003912 environmental pollution Methods 0.000 description 1
- 239000003337 fertilizer Substances 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 229910052742 iron Inorganic materials 0.000 description 1
- UQSXHKLRYXJYBZ-UHFFFAOYSA-N iron oxide Inorganic materials [Fe]=O UQSXHKLRYXJYBZ-UHFFFAOYSA-N 0.000 description 1
- 239000011133 lead Substances 0.000 description 1
- QSHDDOUJBYECFT-UHFFFAOYSA-N mercury Chemical compound [Hg] QSHDDOUJBYECFT-UHFFFAOYSA-N 0.000 description 1
- 229910052753 mercury Inorganic materials 0.000 description 1
- 229910021645 metal ion Inorganic materials 0.000 description 1
- 239000004570 mortar (masonry) Substances 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- NDLPOXTZKUMGOV-UHFFFAOYSA-N oxo(oxoferriooxy)iron hydrate Chemical compound O.O=[Fe]O[Fe]=O NDLPOXTZKUMGOV-UHFFFAOYSA-N 0.000 description 1
- 239000000575 pesticide Substances 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 150000003839 salts Chemical class 0.000 description 1
- 238000010187 selection method Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000005527 soil sampling Methods 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
- 239000002699 waste material Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
- G01N21/31—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N21/55—Specular reflectivity
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- General Health & Medical Sciences (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Biochemistry (AREA)
- Analytical Chemistry (AREA)
- Chemical & Material Sciences (AREA)
- Health & Medical Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Mathematical Physics (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Pure & Applied Mathematics (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Evolutionary Computation (AREA)
- Probability & Statistics with Applications (AREA)
- Operations Research (AREA)
- Algebra (AREA)
- Artificial Intelligence (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
本发明公开了一种基于高光谱遥感技术估算土壤重金属含量的方法,实施步骤如下:1)获取待监测区域的土壤样品,在实验室自然风干与研磨处理后,分别测定土壤样品光谱反射率及土壤铬、钴、铜和砷四种重金属元素含量;2)对原始光谱数据进行Savitzky-Golay平滑滤波和光谱重采样后,利用不同的光谱变换方法,将反射率数据与重金属元素数据进行相关性分析,获得各重金属元素建模波段;3)利用不同的建模方法对反射率数据与重金属含量进行回归方程拟合,构建出高光谱估算模型,并依据建模决定系数、验证决定系数、均方根误差和相对分析误差四个评价指标对每种重金属元素下的模型进行评价,得到每种重金属元素最佳的估算模型。
Description
技术领域
本发明涉及高光谱遥感技术估算土壤重金属含量领域,具体涉及一种利用不同的建模方法对反射率数据与重金属含量进行回归方程拟合并构建出高光谱估算模型的研究方法。
背景技术
土壤是由矿物质、有机质、水分、空气等多种物质组成的复杂的自然综合体,不仅是大自然给予人类的宝贵又有限的自然资源,同时也是人类生存环境的重要载体,因此土壤的健康对于地球上生物的生存具有重要的意义。然而,随着工业 的发展和农业生产的现代化,大量废弃污染物通过各种媒介进入到土壤环境,使 得土壤受到不同程度的污染,其中重金属是最为重要的污染物质之一。重金属是指密度或比重大于5.0g/cm3的金属,比如汞、镉、铅、铬、砷、铜、钴等元素。重金属元素进入到土壤当中,难以被土壤中的生物降解,并且超过自然条件下含量的临界值即会在土壤中发生累积,其不仅会造成农作物减产,质量下降,并且会通过食物链对人体和动植物产生危害,土壤重金属污染已经成为中国最主要的土壤环境污染问题之一。
防治土壤重金属污染的核心问题在于如何快速准确地获取重金属含量的多少及其空间分布信息,传统的土壤重金属监测方法主要以野外采样和室内化学分析为主,其具有测量精度高、准确性强等优点,但是同时也费时费力,并且难以对大范围区域进行连续监测。高光谱遥感是指在可见光与红外波段范围内,获得很多窄的且光谱连续的数据的技术,利用该技术可以获取地物精细的光谱信息,且具有非破坏性、非接触元素、大范围快速监测等传统方法无可比拟的优点,因此高光谱遥感技术自从诞生以来便得到了十分广泛的应用,尤其是近几年运用高光谱技术监测土壤重金属含量,其原理是因为土壤成分的不同,反映在在光谱曲线上也会存在差异,而高光谱技术能够对土壤在不同光谱波段上进行连续的测量,且光谱分辨率高,能够探测到土壤成分间的细微信息。
在土壤重金属元素模型估算过程中,土壤光谱曲线不可避免地存在测量误差,从而对模型精度产生影响,且原始光谱数据所包含的信息量有限,因此为了充分挖掘出光谱信息,国内外学者们运用了多种光谱变换方法以期提高土壤重金属元素的反演结果,目前常用的光谱变换方法主要以下几类:①光谱重采样:由于原始光谱数据光谱分辨率高,建模时为了消除数据冗余,需要对光谱进行重采样处理。此外,因为土壤重金属元素其光谱敏感段比较狭长,因此不需要极高的光谱分辨率;②微分变换:微分变换是光谱变换中常用的方法之一,经过微分变换能够增强原始光谱中变化剧烈波段的光谱信息,而这些变化剧烈的波段极有可能是土壤重金属元素的敏感波段,因此在实际研究中被众多学者们采纳,但同时也需注意并非所有土壤重金属元素经过微分变换后能够取得较好的结果;③倒数与对数及均方根变换:前人研究表明土壤原始光谱经过对数、倒数或者均方根变换后 能够降低噪音引起的光谱误差,并且能够增强光谱之间的差异,从而能够在一定程度上提高建模的精度;④多元散射校正与连续统去除:多元散射校正和连续统去除能够降低土壤原始光谱中背景信息所造成的影响,增强土壤光谱中的吸收特征和反射特征;⑤光谱平滑滤波:由于原始光谱在测量时必定存在一定的误差,使得最终的原始光谱曲线存在许多“毛刺”噪声,因此对原始光谱数据进行滤波处理能够显著降低这些噪声的信息,目前最常使用的滤波方法是Savitzky-Golay滤波方法,虽然这些方法能够在一定程度上提高土壤重金属元素反演结果,但是在前人研究中可以发现对于不同的研究区不同的土壤重金属元素,光谱变换方法并没有普适性;因此在实际操作过程中需要分析不同光谱变换方法与土壤重金属元素之间的关系,以期能够选择出区域最合适的光谱变换方法。
由于土壤重金属元素含量一般很低,况且土壤是十分复杂的组成体,因此土壤重金属元素对于土壤光谱的响应过程也极其复杂,用物理模型很难反演得到土壤重金属含量。在实际操作过程中,常用统计方法来分析土壤光谱与重金属含量之间的相关性,并进一步通过相关建模方法反演得到区域不同土壤重金属元素的含量,常用的建模方法主要有以下几种:①统计回归分析法:统计回归方法主要分为一元回归分析和多元逐步回归分析;②主成分回归法:由于原始光谱数据信息冗余,而主成分回归分析是通过线性变换,将原始光谱信息进行压缩,使得在 降低数据量的同时最大化地将光谱信息压缩到几个互不相关的信息波段中,因 此在高光谱分析中被众多学者所采用;③偏最小二乘法:与传统的最小二乘法相 比,该方法综合了最小二乘法与主成分分析法,由于其偏爱与因变量有关的部分,也被称为偏最小二乘法;④地理加权回归法:该方法核心思路是采样点对其附近点特征的影响要大于离其较远的点,点位的回归系数不再是利用全局信息获得,而是利用邻近观测值的样本数据进行局部回归估计得到,且将数据的空间位置嵌入到回归参数中,使得对空间估值问题更合乎实际。
目前利用高光谱遥感技术估算土壤重金属元素含量的研究中,光谱处理方法及建模方法较多,如研究中发现一阶导数光谱变换方式和偏最小二次回归建模法效果比较好,但在实际过程中,由于土壤背景条件以及区域土壤重金属种类的不同,土壤最佳光谱变换方式以及最佳估算模型仍会存在差异,并且现有研究方法中缺乏彼此之间的对比分析,因此本发明在现有研究不足的基础上,分析了不同的光谱变换方式对研究区土壤重金属元素的影响以及探讨不同模型构建方法对研究区土壤重金属元素的适用性以及模型精度的影响,为进一步研究高光谱遥感技术估算土壤重金属含量等提供参考。
发明内容
本发明要解决的技术问题是提供一种基于高光谱遥感技术估算土壤重金属含量的方法,其具有估算结果准确、通用性好、建模结果可靠实用、适用范围广、直观性强的特点。
为解决上述技术问题,本发明采用的技术方案为:
一种基于高光谱遥感技术估算土壤重金属含量的方法,其实施步骤如下:
1)获取待监测区域的土壤样品,在实验室自然风干与研磨处理后,分别测定土壤样品光谱反射率及土壤铬、钴、铜和砷四种重金属元素含量;
2)对步骤1)获得的土壤样品光谱反射率及土壤铬、钴、铜和砷四种重金属元素含量的原始光谱数据进行Savitzky-Golay平滑滤波和光谱重采样后,利用不同的光谱变换方法,将反射率数据与重金属元素数据进行相关性分析,获得各重金属元素建模波段;
3)利用不同的建模方法对步骤1)获得的反射率数据与重金属含量进行回归方程拟合,构建出高光谱估算模型,并依据建模决定系数、验证决定系数、均方根误差和相对分析误差四个评价指标对每种重金属元素下的模型进行评价,得到每种重金属元素最佳的估算模型,在此基础上,对四种土壤重金属元素最佳估算模型和最佳光谱变换方式进行概况分析;
作为本发明上述技术方案的进一步改进:
所述步骤1)中获取待监测区域的土壤样品,在实验室自然风干与研磨处理,分别测定土壤样品光谱反射率及土壤铬、钴、铜和砷四种重金属元素含量的详细步骤包括:
A1)在ArcGIS10.3软件中采用网格随机抽样法确定采样点位置,并于室内结合待估算区域高分二号与北京一号高分辨率遥感影像、高程数据及土壤类型数据,采集土壤样本点,土壤采集过程中,取样质量为1kg,并且采集0-20cm的表层土壤,采集完成后立即用大号自封袋密封保存好,贴入样本编号,并依据Trimble Juno SB手持GPS进行定位并记录其坐标与土样周围地物类型的情况,之后放在实验室室内进行自然风干,并在表面盖上塑料防护罩,防止扬尘对土样进行污染,风干完成后,剔除砂砾、石块、植物根系及蚯蚓残渣,用玛瑙研钵研磨处理,将土壤样本分成两份,分别过100目筛和200目筛,进行后续的处理;
A2)土壤光谱数据采用美国Analytical Spectral Device公司生产的FieldSpec3光谱仪进行测量,首先,需要将仪器提前开启进行预热准备,时间为30分钟;其次,准备好三脚架,将1000W的卤光灯作为光源,并采用60度的光源照射角度与25度的视场角;最后,将过100目筛的土壤样品放置于光谱反射率极低的外表包有黑色背景材料的盛样皿中,且盛样皿距离光谱探头15cm,并利用40cm×40cm的白板进行定标,获取绝对反射率;准备工作结束后,依次测量每个土壤样本的光谱,各土壤样本获取十条光谱曲线,以消除测量过程中的不稳定性,此外,每测5个土壤样本,利用白板进行定标处理;测量结束后,利用ViewSpec Pro光谱处理软件查看各样本光谱并剔除异常光谱曲线,并对各土壤样本计算得到光谱反射率平均值,同时将其作为原始光谱反射率值;
A3)测定土壤重金属含量,将过200目筛的土壤样本装入自封袋中,并依次准确称取0.0400g样品放置于PTFE内胆中,加入1.5ml氢氟酸、0.5ml硝酸后,将内胆密封,并放置于涂有PTFE涂层的防腐高效溶样罐套内,于烘箱中150°C加热12小时;冷却后取出PTFE内胆,开盖,加入0.25m高氯酸,于150°C电热板上,蒸至近干之后,加入2ml高纯水和1ml硝酸,再次密封,并将内胆置于防腐高效溶样罐套内,于150°C烘箱中回溶12小时;冷却后,取出PTFE内胆,将溶液转移到40ml一次性痰杯中,用高纯水稀释至40ml后摇匀;之后,倒取10ml于一次性注射器中,并通过孔径为0.45m的聚醚砜水系滤头将其注射到10ml离心管中,再次摇匀并写上标号后,采用X-SerieII电感耦合等离子体质谱仪测定钒、铬、钴、镍、铜、锌、砷、镉和铅共9种重金属元素含量,其中每个样本在测量时共测量两次,选取平均值作为最终的测量结果;测量结束后,需要对其计算以得到最终的土壤重金属含量结果,根据得到土壤重金属含量,其中是土壤样本重金属含量值,是土壤样本测量值,是空白检验样本平均值,是定容重量,是土壤样品称取重量;
所述步骤2)中对原始光谱数据进行Savitzky-Golay平滑滤波和光谱重采样后,利用不同的光谱变换方法,将反射率数据与重金属元素数据进行相关性分析,获得各重金属元素建模波段的详细步骤包括:
B1)使用Savitzky-Golay平滑滤波方法处理光谱数据,Savitzky-Golay平滑滤波方法是由第一作者Savitzky于1964年提出的一种滤波方法,该方法是一种在时域内基于多项式,通过移动窗口,利用最小二乘法进行最佳拟合的方法;这是一种直接处理来自时间域内数据平滑问题的方法,与先在频域内定义特性再转到时间域的滤波器不同,计算机在滤波过程中仅是充当一个平滑噪音起伏的滤波器并尽量保证原始数据不失真;因此,相对于其他类似的平均方法而言,Savitzky-Golay滤波法能够保留相对极大值、极小值和宽度等分布特征;
B2)进行光谱重采样,由于FieldSpec3光谱仪在输出光谱结果时,是以1nm为间隔进行重采样,去除首尾噪音波段,共有2001个波段,这使得土壤光谱曲线相邻波段之间信息总存在重合,使得最终整个光谱数据存在冗余;因此,为了消除冗余,采用10nm为间隔进行光谱重采样,即每10个波段进行算术平均,并以第一个波段值作为计算后的波段值,如对400-409nm进行算术平均,计算结果为新的400nm波段的光谱值;通过光谱重采样处理,不仅能够降低冗余,所得曲线更加平滑,而且维持了SG滤波后初始光谱之间的特征;
B3)进行多元散射校正,多元散射校正是一种多变量散射校正技术,能有效地消除样品间散射影响所导致的基线平移和偏移现象从而提高光谱的信噪比;原理是通过数学方法将光谱中的组分吸收信息与光谱散射信息进行分离,通过校正光谱的散射从而得到较为理想的光谱曲线;具体计算过程如下,首先建立一个待测样品的“理想光谱”,以其为标准对其余样品的光谱曲线进行修正,其中包括基线平移和偏移校正,但是在实际应用中,难以得到“理想光谱”,因此对所有样本光谱曲线取平均值是完全可以实施的;其次将各土壤样本光谱曲线与所有样本光谱曲线平均值进行一元线性回归运算,得到一元线性拟合方程的截距和斜率,即线性平移量和倾斜偏移量;最后,将各土壤样本光谱曲线减去对应的线性平移量后再除以倾斜偏移量;由于各土壤光谱曲线均在标准光谱的参考下进行修正,而土壤样品中的成分含量所对应的光谱影响信息在整个校正过程中没有任何影响,因此校正结果提高了土壤光谱曲线的信噪比;
B4)进行光谱变换与导数分析,光谱变换是指通过对待变换光谱数据进行一系列的数学运算,如均方变换、倒数变换、倒数对数变换,变换结果能够在一定程度上更好地选择出敏感波段,从而能够提高建模精度;导数分析,是高光谱分析中常用的预处理方法,其对降低背景信息、分解混合光谱非常有效,并且经过导数分析的光谱结果能够提供比原始光谱更高的分辨率和更清晰的光谱轮廓变换,达到增强光谱信息的目的,其中,一阶导数和二阶导数是使用较多的变换方法,一阶导数的计算公式为,二阶导数的计算公式为其中是波段i的波长值,和分别是波长处的一阶和二阶导数,、、、分别是相应波长处的反射率值,是波长到的差值,由光谱采样间隔决定;
B5)进行波段深度及连续统去除,连续统去除法能够降低背景信息的影响,从而更加有效地突出光谱曲线的吸收和反射特征。连续统去除实质上就是用实际光谱反射率值去除包络线上相应波段反射率值,其特点就是将反射率光谱归一化到0-1之间的值,经过连续统去除后,原先峰值点上对应的反射率值均变为1,非峰值点上对应的反射率值则小于1;通过连续统去除可以间接得到波段深度,即每一样点的土壤反射率归一化到对应的光谱背景上,有利于光谱曲线之间特征波段的比较,根据和进行计算,其中为波长处连续统去除值,是经过SG滤波后波长处的初始光谱值,是对应波长的去包络线值,BD是波长处对应的波段深度;
所述步骤3)中利用不同的建模方法对反射率数据与重金属含量进行回归方程拟合,构建出高光谱估算模型的详细步骤包括:
C1)利用统计回归分析法估算土壤重金属含量,其核心是以回归方程为理论基础,通过多种数学方法建立土壤重金属元素与对应的高光谱数据之间的联系,一元回归分析方法计算较为简单,即在土壤高光谱数据中找出与土壤重金属含量相关系数绝对值最高的波段作为自变量,对应的土壤重金属含量作为因变量,之后通过一定的数学公式进行回归方程的拟合,具体有一元线性函数、幂函数、指数函数、对数函数、逆函数;多元逐步回归分析基本思路为先依据各自变量的重要性,之后每一步选择一个重要的变量进入回归方程;多元逐步回归不仅考虑到按贡献度大小逐一挑选重要变量,而且还考虑到较早进入回归方程的某些变量,随着其后一些变量的选入而失去原有的重要性,这样的变量应及时从回归方程中被剔除,从而使得回归方程中始终只保留重要的变量;
C2)利用主成分回归法估算土壤重金属含量,主成分回归法是通过线性变换,将原来多个指标组合成相互独立的少数几个能充分反映总体信息的主成分,从而在不丢掉重要信息的前提下,避开变量间的共线性,偏于进一步分析,其核心思想为降维;
C3)利用偏最小二乘法估算土壤重金属含量,偏最小二乘法的算法基础为最小二乘法,因为其尽可能提取包含自变量更多信息的基础上,保证了提取成分与因变量间最大相关性,即偏爱与因变量有关的部分,在使用偏最小二乘法建模时,需要正确选取主因子数量,因为主因子数量的多少与模型最终的实际预测能力有显著关系,在本发明中,所采用的确定最佳主因子数量的方法为常用的留一交互验证法,该方法的基本原理为若样本集中有N个样本,每次从中剔除一个样本,用剩下的N-1个样本作为训练集拟合模型,然后用这个拟合模型预测被剔除样本的组分含量,之后,用同样的方法对样本集N中每个样本都剔除一遍,并且得到它们的预测值;通过这一方法进行预测产生的误差和称为PRESS,根据 计算得到,其中是使用留一交叉验证法计算得到的第i个样本的预测值,是第i个样本的实测值,留一交叉验证法的均方根误差根据计算得到,通常情况下,当RMSECV取得最小值时,可以确定模型的最佳主因子数;
C4)利用地理加权回归法估算土壤重金属含量,地理加权回归法是由第一作者Fotheringham提出的一个空间变参数模型,该模型是一种对不同空间子区域自变量和因变量之间关系随空间变化进行建模的非参数局部空间回归分析方法;其核心思路是采样点对其附近点特征的影响要大于离其较远的点,点位的回归系数不再是利用全局信息获得,而是利用邻近观测值的样本数据进行局部回归估计得到,且将数据的空间位置嵌入到回归参数中,使得对空间估值问题更合乎实际;
所述步骤3)中利依据建模决定系数、验证决定系数、均方根误差和相对分析误差四个评价指标对每种重金属元素下的模型进行评价,得到每种重金属元素最佳的估算模型的详细步骤包括:
D1)采用决定系数评价模型的精度,决定系数主要用于评价模型建模与验证结果与1:1趋势线的拟合程度,并衡量模型精度的稳定性和准确度;决定系数的大小决定了相关的密切程度,当越接近于1,表示模型的效果越好,相反,越接近于0,则表示模型的效果越差,根据计算得到决定系数,其中是样本实测值,为样本预测值,为样本平均值,n为样本数量;
D2)采用均方根误差评价模型的反演能力,均方根误差主要评价模型的反演能力,其值越小表明实测值与预测值之间的偏差越小,两者的一致性也越好,也就是模型的模拟结果越精确可靠,根据计算得到均方根误差,其中是样本实测值,为样本预测值,n为样本数量;
D3)采用相对分析误差检验模型的预报能力,相对分析误差是验证样本实测值的标准差与验证样本预测值均方根误差的商,根据计算得到,第一作者Viscarra对RPD进行了标准的划分,当PRD<1.4时,模型无法对样本进行预测;当1.4≤RPD<2时,模型效果一般,可以用此模型对样本进行粗略评估;但RPD≥2时,模型具有极好的预测能力;
本发明具有下述优点:
1、本发明克服了传统的土壤重金属监测方法费时费力,且难以实现区域土壤重金属含量的快速监测的缺点,使用高光谱遥感技术快速获取地物精细的光谱信息,从而获取区域土壤重金属含量,实现区域土壤重金属含量的快速检测。
2、本发明通过对土壤样品反射率光谱曲线进行研究,且将其与土壤重金属元素进行相关性分析,获得各重金属元素的敏感波段,从而利用不同的建模方法进行高光谱模型的估算,为后续相关部门大范围提取土壤重金属污染信息提供一定的参考。
3、本发明运用了多种光谱变换方法提高土壤重金属元素的反演结果,克服了在土壤重金属元素模型估算过程中,土壤光谱曲线存在的测量误差,保证了模型的精度并充分挖掘出光谱信息。
4、本发明使用Savitzky-Golay平滑滤波方法处理光谱数据,相对于其他类似的平均方法而言,该方法能够保留相对极大值、极小值和宽度等分布特征,经过滤波后的光谱曲线明显消除了“毛刺”噪声,较好地保存原始光谱曲线的总体特征。
5、本发明以10nm为间隔进行光谱重采样,即每10个波段进行算术平均,并以第一个波段值作为计算后的波段值,不仅能够降低冗余,所得曲线更加平滑,而且维持了SG滤波后初始光谱之间的特征。
6、本发明使用倒数与对数及均方根变换方法降低了噪音引起的光谱误差,并且能够增强光谱之间的差异,从而能够在一定程度上提高建模的精度;使用多元散射校正与连续统去除法降低土壤原始光谱中背景信息所造成的影响,增强了土壤光谱中的吸收特征和反射特征。
7、本发明使用主成分回归法将原始光谱信息进行压缩,使得在降低数据量的同时最大化地将光谱信息压缩到几个互不相关的信息波段中,解决了原始光谱数据信息冗余的缺点。
8、本发明使用偏最小二乘法进行建模,与传统的最小二乘法相比,该方法综合了最小二乘法与主成分分析法,消除了多重共线性问题,能够最大程度地概括土壤光谱信息,提高模型精度,建模所得的结果可以广泛的应用于农业、环保、水利、国土等部门。
附图说明
图1为本发明实施例的基本流程示意图。
图2为本发明实施例得到的土壤采样点光谱反射率图。
图3为本发明实施例得到的Savitzky-Golay滤波后的土壤采样点光谱反射率图。
图4为本发明实施例得到的SG滤波后的土壤光谱反射率重采样结果图。
图5为本发明实施例得到的土壤样本光谱平均反射率结果图。
图6为本发明实施例得到的经过滤波和光谱重采样后的原始光谱及一阶导数和二阶导数与土壤钴含量的相关系数曲线图。
图7为本发明实施例得到的初始光谱均方根变换值及均方根变换后的一阶导数和二阶导数与土壤钴含量的相关系数曲线图。
图8为本发明实施例得到的初始光谱倒数变换值及倒数变换后的一阶导数和二阶导数与土壤钴含量的相关系数曲线图。
图9为本发明实施例得到的初始光谱波段深度及多元散射校正值与土壤钴含量的相关系数曲线图。
图10为本发明实施例得到的土壤钴元素不同回归模型实测值与预测值散点图,其中:a)为一元回归模型实测值与预测值散点图;b)为多元逐步回归模型实测值与预测值散点图;c)为偏最小二乘回归模型实测值与预测值散点图。
具体实施方式
如图1所示,本发明以木兰溪下游区为具体实施例来证明本发明的方法是可行的:
一、 土壤重金属含量测定与分析:
本发明实施例中使用的土壤样本数据是于2017年4月3日去木兰溪下游研究区进行实地观测,记录主要土地利用类型,并于室内结合研究区高分二号与北京一号高分辨率遥感影像、高程数据及土壤类型数据,在ArcGIS 10.3软件中采用网格随机抽样法确定采样点位置,并于2017年5月16 日与5月17日赴研究区采集43个土壤样本点得到的;土壤光谱数据则采用美国Analytical Spectral Device公司生产的FieldSpec3光谱仪进行测量,该光谱仪光谱波段测量范围为350-2500nm,光谱采样间隔在350-1000nm波段范围时为1.4nm;1000-2500nm波段范围时为2nm;光谱分辨率在350-1000nm范围内为3nm,1000-2500nm范围内为10nm。
1.1木兰溪下游土壤样品采集
土壤采集过程中,取样质量为1kg左右,并且主要采集0-20cm的表层土壤,采集完成后立即用大号自封袋密封保存好,贴入样本编号,并依据Trimble Juno SB手持GPS进行定位并记录其坐标与土样周围地物类型的情况,之后放在实验室室内进行自然风干,并在表面盖上塑料防护罩,防止扬尘对土样进行污染,风干完成后,剔除砂砾、石块、植物根系及蚯蚓等残渣,用玛瑙研钵研磨处理,将土壤样本分成两份,分别过100目筛和200目筛,进行后续的处理。需要注意的是,研究区土壤样本点土壤分布于木兰溪下游两岸1000m以内,其主要土地利用类型为水田和水浇地。由于采样点分布区域受人类扰动程度大,土壤利用程度高,尤其是农田区域化肥以及农药的使用量较多,因此其重金属元素的富集程度一般要高于人类扰动程度低的区域。
1.2木兰溪下游土壤光谱数据获取
土壤光谱数据采用美国Analytical Spectral Device公司生产的FieldSpec3光谱仪进行测量,由于外界环境对光谱测量的结果有着极大的影响,因此测量土壤光谱前,需要一系列的准备工作;首先,需要将仪器提前开启进行预热准备,时间为30分钟左右;其次,准备好三脚架,将1000W的卤光灯作为光源,并采用60度的光源照射角度与25度的视场角;最后,将过100目筛的土壤样品放置于光谱反射率极低的外表包有黑色背景材料的盛样皿中,且盛样皿距离光谱探头15cm,并利用40cm×40cm的白板进行定标,获取绝对反射率;准备工作结束后,依次测量每个土壤样本的光谱,各土壤样本获取十条光谱曲线,以消除测量过程中的不稳定性,此外,每测5个土壤样本,利用白板进行定标处理。测量结束后,利用ViewSpec Pro光谱处理软件查看各样本光谱并剔除异常光谱曲线,并对各土壤样本计算得到光谱反射率平均值,同时将其作为原始光谱反射率值。
1.3木兰溪下游土壤重金属含量测定
将过200目筛的土壤样本装入自封袋中,并依次准确称取0.0400g样品放置于PTFE内胆中,加入1.5ml氢氟酸、0.5ml硝酸后,将内胆密封,并放置于涂有PTFE涂层的防腐高效溶样罐套内,于烘箱中150°C加热12小时。冷却后取出PTFE内胆,开盖,加入0.25m高氯酸于150°C电热板上,蒸至近干之后,加入2ml高纯水和1ml硝酸,再次密封,并将内胆置于防腐高效溶样罐套内,于150°C烘箱中回溶12小时。冷却后,取出PTFE内胆,将溶液转移到40ml一次性痰杯中,用高纯水稀释至大约40ml后摇匀。之后,倒取约10ml于一次性注射器中,并通过孔径为0.45m的聚醚砜水系滤头将其注射到10ml离心管中,再次摇匀并写上标号后,采用X-SerieII电感耦合等离子体质谱仪测定钒、铬、钴、镍、铜、锌、砷、镉和铅共9种重金属元素含量,其中每个样本在测量时共测量两次,选取平均值作为最终的测量结果,测量结束后,需要对其计算以得到最终的土壤重金属含量结果,其表达式为:
式中是土壤样本重金属含量值,是土壤样本测量值,是空白检验样本平均值,是定容重量,是土壤样品称取重量。
1.4木兰溪下游土壤重金属含量分析
计算土壤重金属含量结束后,利用SPSS22.0统计软件对其进行描述性统计分析,结果如表1所示,其中背景值采用的是第一作者陈振金测量得到的福建省土壤背景值。
表1 木兰溪下游土壤重金属描述性统计
依据表1可以得到,研究区各重金属元素含量平均值分别为:钒为81.95mg/kg、铬为48.43mg/kg、钴为10.88mg/kg、镍为21.90mg/kg、铜为23.97mg/kg、锌为11.37mg/kg、砷为8.75mg/kg、镉为0.27mg/kg、铅为44.90mg/kg,将其与福建省土壤背景值进行比较,发现均高于土壤背景值。各重金属元素与土壤背景值之间的倍数关系分别为:钒为背景值1.04倍、铬为背景值的1.17倍、钴为背景值的1.46倍、镍是背景值的1.62倍、铜是背景值的1.10倍、锌是背景值的1.37倍、砷为背景值的1.51倍、镉是背景值的5.4倍、铅是背景值的1.28倍。变异系数是表征样品变异程度的指标,能在一定程度上反映样品受人为影响的程度。当变异系数小于10%时,表明样品呈现弱度变异;当变异系数介于10%-30%之间时,表明样品呈现中等变异;当变异系数大于30%时,表明样品呈现强度变异。根据划分标准,可以得到研究区钒、钴、锌、砷和铅为中等变异,其余重金属元素为强度变异,变异程度最大的为镍元素。整体表明研究区土壤重金属含量可能与区域人类活动及城市建设活动程度较为剧烈有一定的关系。根据峰度和偏度值,可以得到除镍、锌、铅和镉外,其余重金属元素峰度和偏度值均小于1,表明研究区土壤重金属含量近似服从正态分布。此外,大部分重金属元素标准差均较大,表明研究区土壤重金属空间分布呈现点源聚集性特征。
二、木兰溪下游区光谱数据处理与建模方法评价:
由于土壤是各种物理和化学性质各不相同的物质组成,因此在光谱曲线上也各不相同。研究表明影响土壤光谱的重要因素包括水分、有机质、氧化铁、土壤质地等,其中土壤中水分、铁、有机质等含量的增加会降低其反射率,盐分、粘粒、矿质粉粒含量的增加会提高其反射率。第一作者戴昌达通过对我国23 种土壤类型的光谱曲线进行研究,并依据光谱曲线的形状、深度、斜率等特征将我国土壤光谱曲线分为四种类型,分别是平直型光谱曲线、缓斜型光谱曲线、陡坎型光谱曲线与波浪型光谱曲线。通过对研究区43个土壤样木曲线进行平均值计算后,得到如图2所示的光谱曲线图,依据图2可以得到研究区土壤样本在350~2500 nm区间内,其反射光谱曲线整体趋势走向相似,即在可见光波段,反射光谱曲线斜率较大,在近红外波段,反射光谱曲线斜率较小,甚至接近于水平,之后随着波长的继续增加,土壤反射率逐渐下降,其中在1400 nm、1900 nm和2200 nm附近形成明显的光谱吸收峰。己有研究表明,土壤光谱在可见光与近红外波段的吸收特征主要是因为金属离子的电子跃迁,而在短波红外区域的吸收特征则主要是因为有机质、层状硅酸盐、碳酸盐、硫酸盐等矿物质的各类分子团中化学键的伸展、弯曲、变形等振动。
2.1木兰溪下游光谱数据处理
虽然土壤反射光谱数据是在室内严格的模拟环境下测量得到,但是仍会不可避免受到测试环境、仪器本身、样品背景以及杂散光等因素的影响,导致光谱曲线不够平滑,为了消除这些无关信息,充分挖掘出土壤光谱数据中隐含的有效信息,因此需要对土壤光谱数据采取一定的处理步骤。在本发明中,处理步骤主要有Savitzky-Golay平滑滤波、光谱重采样、多元散射校正、光谱变换与导数分析、波段深度与连续统去除。
平滑滤波
使用Savitzky-Golay平滑滤波方法处理光谱数据,Savitzky-Golay平滑滤波方法是由第一作者Savitzky于1964年提出的一种滤波方法,该方法是一种在时域内基于多项式,通过移动窗口,利用最小二乘法进行最佳拟合的方法。这是一种直接处理来自时间域内数据平滑问题的方法,与先在频域内定义特性再转到时间域的滤波器不同,计算机在滤波过程中仅是充当一个平滑噪音起伏的滤波器并尽量保证原始数据不失真。因此,相对于其他类似的平均方法而言,Savitzky-Golay滤波法能够保留相对极大值、极小值和宽度等分布特征;通过利用Origin 9.1软件对土壤光谱数据进行Savitzky-Golay滤波处理,滤波参数设置为Points of Window等于5,Polynomial Order等于3,即滤波器为“五点三次”滤波器,得到最终43个土壤样本光谱反射曲线滤波结果,如图3所示。由图3可以得到,经过Savitzky-Golay滤波后的土壤光谱反射率曲线表现更加平滑,通过对其中一个样本的土壤反射光谱曲线进行滤波前后对比,并以2000~2400 nm为对比波段,将其放大显示,可以发现,经过滤波后的光谱曲线明显消除了“毛刺”噪声,且 Savitzky-Golay滤波较好地保存原始光谱曲线的总体特征。
光谱重采样
由于FieldSpec3光谱仪在输出光谱结果时,是以1nm为间隔进行重采样,去除首尾噪音波段,共有2001个波段,这使得土壤光谱曲线相邻波段之间信息存在重合,使得最终整个光谱数据存在冗余。因此,为了消除冗余,一般采用10nm为间隔进行光谱重采样,即每10个波段进行算术平均,并以第一个波段值作为计算后的波段值,如对400-409nm进行算术平均,计算结果为新的400nm波段的光谱值。通过光谱重采样处理,不仅能够降低冗余,所得曲线更加平滑,而且维持了SG滤波后初始光谱之间的特征。本发明对以上预处理过程中生成的14种光谱预处理结果进行相同的光谱重采样操作,并以新生成的201个光谱波段进行后续的建模研究。以SG滤波后的初始光谱为例,对研究区43个样本进行光谱重采样,得到新的曲线如图4所示。
多元散射校正
多元散射校正,简称MSC法,是一种多变量散射校正技术,其可以有效地消除样品间散射影响所导致的基线平移和偏移现象从而提高光谱的信噪比。主要原理是通过数学方法将光谱中的组分吸收信息与光谱散射信息进行分离,通过校正光谱的散射从而得到较为理想的光谱曲线。具体计算过程如下,首先建立一个待测样品的“理想光谱”, 以其为标准对其余样品的光谱曲线进行修正,其中包括基线平移和偏移校正,但是在实际应用中,难以得到“理想光谱”,因此对所有样本光谱曲线取平均值是完全可以实施的;其次将各土壤样本光谱曲线与所有样本光谱曲线平均值进行一元线性回归运算,得到一元线性拟合方程的截距和斜率,即线性平移量和倾斜偏移量;最后,将各土壤样本光谱曲线减去对应的线性平移量后再除以倾斜偏移量。由于各土壤光谱曲线均在标准光谱的参考下进行修正,而土壤样品中的成分含量所对应的光谱影响信息在整个校正过程中没有任何影响,因此校正结果提高了土壤光谱曲线的信噪比。图5是43个土壤样本取平均值后的土壤反射光谱曲线,表2是43个土壤样本计算得到的线性平移量和倾斜偏移量。
表2 各土壤样本线性平移量和倾斜偏移量
2.1.4光谱变换与导数分析
光谱变换是指通过对待变换光谱数据进行一系列的数学运算,如均方变换、倒数变换、倒数对数变换等,变换结果能够在一定程度上可以更好地选择出敏感波段,从而能够提高建模精度。导数分析又称为微分分析,是高光谱分析中常用的预处理方法,其对降低背景信息、分解混合光谱非常有效,并且经过导数分析的光谱结果能够提供比原始光谱更高的分辨率和更清晰的光谱轮廓变换,达到增强光谱信息的目的,其中,一阶导数和二阶导数是使用较多的变换方法,一阶导数的表达式为:
二阶导数的表达式式为:
式中,是波段i的波长值,和分别是波长处的一阶和二阶导数,、、、分别是相应波长处的反射率值,是波长到的差值,由光谱采样间隔决定;本发明中所采用的光谱变换方法包括均方根变换、倒数变换、倒数对数变换、对数倒数变换、初始光谱一阶导数、均方根一阶导数、倒数一阶导数、倒数对数一阶导数、对数倒数一阶导数、初始光谱二阶倒数、均方根二阶导数、倒数二阶导数、倒数对数二阶导数、对数倒数二阶导数。利用木兰溪下游区土壤样本经过Savitzky-Golay滤波后的原始光谱数据及剩余16种光谱变换数据分别与土壤重金属钴含量进行相关性分析,并在Excel统计软件下作出相关系数曲线图,图6是经过滤波和光谱重采样后的原始光谱(以下记为初始光谱)及一阶导数和二阶导数与土壤钴含量的相关系数曲线图。由图6可以得到,通过0.01极显著性检验水平的波段较多,其中初始光谱自570nm波段之后,均通过了检验,且大部分波段的相关系数绝对值均大于0.6,由此可见初始光谱与钴元素之间具有非常好的拟合关系;初始光谱一阶导数通过0.01极显著性检验水平的波段主要位于470-780nm、860-990nm、1010nm、1030-1220nm、1240-1250nm、1350-1560nm、1590nm、1610-1620nm、1760nm、1860-1910nm、1940-2040nm、2130-2230nm和2300-2320nm,其中有些波段与钴元素呈现显著的正相关,有些波段则与钴元素呈现显著的负相关;初始光谱二阶导数与钴元素相关系数变化曲线震动特征较为明显,并没有像初始光谱和初始光谱一阶导数那样,存在着大范的波段通过0.01极显著性检验水平,而主要是分布在多个较窄的波段区间内。图7是初始光谱均方根变换值及均方根变换后的一阶导数和二阶导数与土壤钴含量的相关系数曲线图,由图7可以得到,经过均方根变换后,相关系数变化曲线变化趋势与变换前差别不大,大部分通过0.01极显著性检验水平的波段相同,但是相关系数具体数值却不一样。图8是初始光谱倒数变换值及倒数变换后的一阶导数和二阶导数与土壤钴含量的相关系数曲线图,由图8可以得到,倒数变换后,相关系数变化曲线基本上与变化前呈现相反的趋势,但是相关系数的绝对值通过0.01极显著性检验水平的波段基本相同。图9是初始光谱波段深度及多元散射校正值与土壤钴含量的相关系数曲线图。
波段深度及连续统去除
连续统去除法又称去包络线法,该方法能够降低背景信息的影响,从而更加有效地突出光谱曲线的吸收和反射特征。连续统去除实质上就是用实际光谱反射率值去除包络线上相应波段反射率值,其特点就是将反射率光谱归一化到0-1之间的值。经过连续统去除后,原先峰值点上对应的反射率值均变为1,非峰值点上对应的反射率值则小于1。通过连续统去除可以间接得到波段深度,即每一样点的土壤反射率归一化到对应的光谱背景上,有利于光谱曲线之间特征波段的比较,其表达式为:
式中,为波长处连续统去除值,是经过SG滤波后波长处的初始光谱值,是对应波长的去包络线值,BD是波长处对应的波段深度;以上连续统去除计算在ENVI5.3软件中完成。
2.2木兰溪下游土壤钴含量高光谱估算模型
在利用建模方法进行模型估算前,需要将43个土壤样本数据分为建模集和验证集两部分,这一步骤也被称为校正集的选取。选取合适的建模集对于模型的稳定性有着重要的影响,这是因为验证集,也被称为预测集,是依据建模集中样本的光谱值与土壤重金属含量之间的数学关系而建立起来的。目前建模集样本选取方法主要有随机法和含量梯度法。随机法是从样本集中随机选取一部分样本作为建模集,剩余样本作为验证集,随机法具有操作简单、不需要进行数据挑选等优点,因此在实际研究中具有广泛的应用,但是该方法也存在着一个重要的缺点,即无法保证建模集的代表性。含量梯度法是先将土壤样本的重金属含量按照大小顺序进行排序后,然后按照一定的间隔依次选择部分样本作为建模集。含量梯度法具有样本代表性高等优点,但是同时也有样本多容易造成巨大的工作量等缺点。因为只有43个土壤样本,因此本发明选择含量梯度法选择验证集,先将土壤重金属含量按照从小到大进行排序,后以3:1为采样间隔选取验证样本,即每隔3个样本选取1个样本作为验证样本。最后共获得33个建模样本,10个验证样本。
土壤钴含量统计回归分析模型
利用统计回归分析法估算土壤重金属含量,其核心是以回归方程为理论基础,通过多种数学方法建立土壤重金属元素与对应的高光谱数据之间的联系,比较常用的统计回归分析方法有一元回归分析和多元逐步回归分析,一元回归分析方法计算较为简单,即在土壤高光谱数据中找出与土壤重金属含量相关系数绝对值最高的波段作为自变量,对应的土壤重金属含量作为因变量,之后通过一定的数学公式进行回归方程的拟合,具体有一元线性函数、幂函数、指数函数、对数函数、逆函数等;多元逐步回归分析基本思路为先依据各自变量的重要性,之后每一步选择一个重要的变量进入回归方程。多元逐步回归不仅考虑到按贡献度大小逐一挑选重要变量,而且还考虑到较早进入回归方程的某些变量,有可能随着其后一些变量的选入而失去原有的重要性,这样的变量也应及时从回归方程中被剔除,从而使得回归方程中始终只保留重要的变量。通过挑选出的建模集的光谱变量值与钴元素在SPSS22.0统计软件中进行一元回归建模,模型结果如表3所示。
表3 土壤钴元素一元回归模型最佳建模结果
以挑选出的3个波段土壤光谱反射率分别作为自变量,土壤钴元素含量作为因变量进行多元逐步回归分析,并将51个波段全部进行多元逐步回归分析,各光谱变换方式最佳建模结果如表4所示。
表4 土壤钴元素多元逐步回归模型最佳建模结果
2.2.2土壤钴含量主成分回归模型
利用主成分回归法估算土壤重金属含量,主成分回归法是通过线性变换,将原来多个指标组合成相互独立的少数几个能充分反映总体信息的主成分,从而在不丢掉重要信息的前提下,避开变量间的共线性,偏于进一步分析,其核心思想为降维。将挑选出的建模集的光谱变量值与钴元素在SPS9.4统计软件中进行主成分回归建模,即将各种光谱变换下的3个波段进行主成分回归分析,最后运用多元逐步回归得到的对数倒数一阶导数920nm和波段深度730nm 两个波段进行建模。
土壤钴含量偏最小二乘回归模型
利用偏最小二乘法估算土壤重金属含量,偏最小二乘法的算法基础为最小二乘法,因为其尽可能提取包含自变量更多信息的基础上,保证了提取成分与因变量间最大相关性,即偏爱与因变量有关的部分,在使用偏最小二乘法建模时,需要正确选取主因子数量,因为主因子数量的多少与模型最终的实际预测能力有显著关系,在本发明中,所采用的确定最佳主因子数量的方法为常用的留一交互验证法,该方法的基本原理为若样本集中有N个样本,每次从中剔除一个样本,用剩下的N-1个样本作为训练集拟合模型,然后用这个拟合模型预测被剔除样本的组分含量,之后,用同样的方法对样本集N中每个样本都剔除一遍,并且得到它们的预测值。通过这一方法进行预测产生的误差和称为PRESS,表达式为:
式中是使用留一交叉验证法计算得到的第i个样本的预测值,是第i个样本的实测值,留一交叉验证法的均方根误差根据计算得到,通常情况下,当 RMSECV取得最小值时,可以确定模型的最佳主因子数。以挑选出的3个波段土壤光谱反射率分别作为自变量,土壤钴元素含量作为因变量进行偏最小二乘回归分析,并运用多元逐步回归得到的对数倒数一阶导数920 nm (LRFD920)和波段深度730 nm(BD730)两个波段进行建模,各光谱变换方式最佳建模结果如表5所示。
表5 土壤钴元素偏最小二乘回归最佳建模结果
2.2.4土壤钴含量地理加权回归模型
地理加权回归法是由第一作者Fotheringham提出的一个空间变参数模型,该模型是一种对不同空间子区域自变量和因变量之间关系随空间变化进行建模的非参数局部空间回归分析方法。其核心思路是采样点对其附近点特征的影响要大于离其较远的点,点位的回归系数不再是利用全局信息获得,而是利用邻近观测值的样本数据进行局部回归估计得到,且将数据的空间位置嵌入到回归参数中,使得对空间估值问题更合乎实际。将挑选出的建模集的光谱变量值与钴元素在GWR 4.0软件中进行地理加权回归建模,即将各种光谱变换下的3个波段进行偏最小二乘回归分析,最后运用多元逐步回归得到的对数倒数一阶导数920 nm和波段深度730 nm两个波段进行建模,由于地理加权回归模型考虑到建模样本的空间坐标值,本发明先在ArcGIS 10.3软件中导入采样点空间坐标,通过增加字段进行“计算几何”功能得到采样点在WGS84_UTM_50N下的投影坐标值,按照GWR 4.0软件格式准备数据进行后续的建模运算。
2.3 木兰溪下游土壤钴含量估算模型评价分析
为了对不同模型结果进行对比分析,需要釆用统一的模型评价标准对其进行评价,从而得到各重金属元素最佳的高光谱估算模型。本发明在综合前人研究的基础上,采用如下指标对模型进行精度评价,具体包括建模决定系数、验证决定系数、均方根误差和相对分析误差,其中建模决定系数和验证决定系数只是计算采用数据不一样,计算原理及公式相同。
决定系数
采用决定系数评价模型的精度,决定系数主要用于评价模型建模与验证结果与1:1趋势线的拟合程度,并衡量模型精度的稳定性和准确度。决定系数的大小决定了相关的密切程度,当越接近于1,表示模型的效果越好,相反,越接近于0,则表示模型的效果越差,表达式为:
式中是样本实测值,为样本预测值,为样本平均值,n为样本数量。
均方根误差
采用均方根误差评价模型的反演能力,均方根误差主要评价模型的反演能力,其值越小表明实测值与预测值之间的偏差越小,两者的一致性也越好,也就是模型的模拟结果越精确可靠,其表达式为:
式中是样本实测值,为样本预测值,n为样本数量。
相对分析误差
采用相对分析误差检验模型的预报能力,相对分析误差是验证样本实测值的标准差与验证样本预测值均方根误差的商,其表达式为:
第一作者Viscarra对RPD进行了标准的划分,当PRD<1.4时,模型无法对样本进行预测;当1.4≤RPD<2时,模型效果一般,可以用此模型对样本进行粗略评估;但RPD≥2时,模型具有极好的预测能力。本发明采用验证的10个土壤样本点对模型进行验证,按照建模决定系数与验证决定系数越大、建模均方根误差与验证均方根误差越小、相对分析误差越大和F检验值越大的原则综合挑选出各建模方法中的最佳模型,结果如表6所示,各种建模方法简称分别为一元回归、多元逐步回归和偏最小二乘回归。
表6 土壤钴元素不同建模方法最佳建模结果
根据表7,得到经过多元逐步回归后得到的两个波段对数倒数一阶导数920 nm和波段深度730 nm建模效果显著,分别在多元逐步回归与偏最小二乘回归中取得最优模型,且模型的建模决定系数和验证决定系数均高达0.6。以验证集实测值作为横坐标,验证集预测值作为纵坐标,作出各建模方法最佳估算模型验证集实测值与预测值散点图,如图10所示。由图10可以发现,对于土壤钴元素,不同建模方法中最佳回归模型效果相近,决定系数均高于0.6,并且趋势线均表现为与1:l参考线斜交,即预测值的低值部分高于实测值,高值部分低于实测值。结合估算方法对比表和验证集中实测值与预测值散点图拟合效果,本发明得出木兰溪下游沿岸土壤钴元素最终估算模型如下所示:
y = 12.899 + 2 575.429 × LRFD920 + 930.98× BD730
依据最终估算模型,得到模型验证样本集中实测值与验证值分布表7,可以得到相对误差介于1.51%~35.79%之间,经计算总平均相对误差为11.57%,表明该模型预测效果较好,能够用于估算研究区土壤钴元素含量。
表7 土壤钴元素最佳估算模型预测结果
Claims (1)
1.一种基于高光谱遥感技术估算土壤重金属含量的方法,其特征在于步骤如下:
1)获取待监测区域的土壤样品,在实验室自然风干与研磨处理后,分别测定土壤样品光谱反射率及土壤铬、钴、铜和砷四种重金属元素含量;
2)对步骤1)获得的土壤样品光谱反射率及土壤铬、钴、铜和砷四种重金属元素含量的原始光谱数据进行Savitzky-Golay平滑滤波和光谱重采样后,利用不同的光谱变换方法,将反射率数据与重金属元素数据进行相关性分析,获得各重金属元素建模波段;
3)利用不同的建模方法对步骤1)获得的反射率数据与重金属含量进行回归方程拟合,构建出高光谱估算模型,并依据建模决定系数、验证决定系数、均方根误差和相对分析误差四个评价指标对每种重金属元素下的模型进行评价,得到每种重金属元素最佳的估算模型,在此基础上,对四种土壤重金属元素最佳估算模型和最佳光谱变换方式进行概况分析;
所述步骤1)中获取待监测区域的土壤样品,在实验室自然风干与研磨处理,分别测定土壤样品光谱反射率及土壤铬、钴、铜和砷四种重金属元素含量的详细步骤包括:
A1)在ArcGIS10.3软件中采用网格随机抽样法确定采样点位置,并于室内结合待估算区域高分二号与北京一号高分辨率遥感影像、高程数据及土壤类型数据,采集土壤样本点,土壤采集过程中,取样质量为1kg,并且采集0-20cm的表层土壤,采集完成后立即用大号自封袋密封保存好,贴入样本编号,并依据Trimble Juno SB手持GPS进行定位并记录其坐标与土样周围地物类型的情况,之后放在实验室室内进行自然风干,并在表面盖上塑料防护罩,防止扬尘对土样进行污染,风干完成后,剔除砂砾、石块、植物根系及蚯蚓残渣,用玛瑙研钵研磨处理,将土壤样本分成两份,分别过100目筛和200目筛,进行后续的处理;
A2)土壤光谱数据采用美国Analytical Spectral Device公司生产的FieldSpec3光谱仪进行测量,首先,需要将仪器提前开启进行预热准备,时间为30分钟;其次,准备好三脚架,将1000W的卤光灯作为光源,并采用60度的光源照射角度与25度的视场角;最后,将过100目筛的土壤样品放置于光谱反射率极低的外表包有黑色背景材料的盛样皿中,且盛样皿距离光谱探头15cm,并利用40cm×40cm的白板进行定标,获取绝对反射率;准备工作结束后,依次测量每个土壤样本的光谱,各土壤样本获取十条光谱曲线,以消除测量过程中的不稳定性,此外,每测5个土壤样本,利用白板进行定标处理;测量结束后,利用ViewSpec Pro光谱处理软件查看各样本光谱并剔除异常光谱曲线,并对各土壤样本计算得到光谱反射率平均值,同时将其作为原始光谱反射率值;
A3)测定土壤重金属含量,将过200目筛的土壤样本装入自封袋中,并依次准确称取0.0400g样品放置于PTFE内胆中,加入1.5ml氢氟酸、0.5ml硝酸后,将内胆密封,并放置于涂有PTFE涂层的防腐高效溶样罐套内,于烘箱中150℃加热12小时;冷却后取出PTFE内胆,开盖,加入0.25m高氯酸,于150℃电热板上,蒸至近干之后,加入2ml高纯水和1ml硝酸,再次密封,并将内胆置于防腐高效溶样罐套内,于150℃烘箱中回溶12小时;冷却后,取出PTFE内胆,将溶液转移到40ml一次性痰杯中,用高纯水稀释至40ml后摇匀;之后,倒取10ml于一次性注射器中,并通过孔径为0.45μm的聚醚砜水系滤头将其注射到10ml离心管中,再次摇匀并写上标号后,采用X-SerieII电感耦合等离子体质谱仪测定钒、铬、钴、镍、铜、锌、砷、镉和铅共9种重金属元素含量,其中每个样本在测量时共测量两次,选取平均值作为最终的测量结果;测量结束后,需要对其计算以得到最终的土壤重金属含量结果,根据得到土壤重金属含量,其中MTv是土壤样本重金属含量值,MMv是土壤样本测量值,EAv是空白检验样本平均值,WCv是定容重量,WSw是土壤样品称取重量;
所述步骤2)中对原始光谱数据进行Savitzky-Golay平滑滤波和光谱重采样后,利用不同的光谱变换方法,将反射率数据与重金属元素数据进行相关性分析,获得各重金属元素建模波段的详细步骤包括:
B1)使用Savitzky-Golay平滑滤波方法处理光谱数据,Savitzky-Golay平滑滤波方法是由第一作者Savitzky于1964年提出的一种滤波方法,该方法是一种在时域内基于多项式,通过移动窗口,利用最小二乘法进行最佳拟合的方法;这是一种直接处理来自时间域内数据平滑问题的方法,与先在频域内定义特性再转到时间域的滤波器不同,计算机在滤波过程中仅是充当一个平滑噪音起伏的滤波器并尽量保证原始数据不失真;因此,相对于其他类似的平均方法而言,Savitzky-Golay滤波法能够保留相对极大值、极小值和宽度等分布特征;
B2)进行光谱重采样,由于FieldSpec3光谱仪在输出光谱结果时,是以1nm为间隔进行重采样,去除首尾噪音波段,共有2001个波段,这使得土壤光谱曲线相邻波段之间信息总存在重合,使得最终整个光谱数据存在冗余;因此,为了消除冗余,采用10nm为间隔进行光谱重采样,即每10个波段进行算术平均,并以第一个波段值作为计算后的波段值,如对400-409nm进行算术平均,计算结果为新的400nm波段的光谱值;通过光谱重采样处理,不仅能够降低冗余,所得曲线更加平滑,而且维持了SG滤波后初始光谱之间的特征;
B3)进行多元散射校正,多元散射校正是一种多变量散射校正技术,其能有效地消除样品间散射影响所导致的基线平移和偏移现象从而提高光谱的信噪比;原理是通过数学方法将光谱中的组分吸收信息与光谱散射信息进行分离,通过校正光谱的散射从而得到较为理想的光谱曲线;具体计算过程如下,首先建立一个待测样品的“理想光谱”,以其为标准对其余样品的光谱曲线进行修正,其中包括基线平移和偏移校正,但是在实际应用中,难以得到“理想光谱”,因此对所有样本光谱曲线取平均值是完全可以实施的;其次将各土壤样本光谱曲线与所有样本光谱曲线平均值进行一元线性回归运算,得到一元线性拟合方程的截距和斜率,即线性平移量和倾斜偏移量;最后,将各土壤样本光谱曲线减去对应的线性平移量后再除以倾斜偏移量;由于各土壤光谱曲线均在标准光谱的参考下进行修正,而土壤样品中的成分含量所对应的光谱影响信息在整个校正过程中没有任何影响,因此校正结果提高了土壤光谱曲线的信噪比;
B4)进行光谱变换与导数分析,光谱变换是指通过对待变换光谱数据进行一系列的数学运算,变换结果能够在一定程度上更好地选择出敏感波段,从而能够提高建模精度;导数分析,是高光谱分析中常用的预处理方法,其对降低背景信息、分解混合光谱非常有效,并且经过导数分析的光谱结果能够提供比原始光谱更高的分辨率和更清晰的光谱轮廓变换,达到增强光谱信息的目的,其中,一阶导数和二阶导数是使用较多的变换方法,一阶导数的计算公式为二阶导数的计算公式为其中λi是波段i的波长值,R′(λi)和R″(λi)分别是波长λi处的一阶和二阶导数,Ri-1、Ri+1、Ri-2、Ri+2分别是相应波长处的反射率值,Δλ是波长λi到λi+1的差值,由光谱采样间隔决定;
B5)进行波段深度及连续统去除,连续统去除法能够降低背景信息的影响,从而更加有效地突出光谱曲线的吸收和反射特征;连续统去除实质上就是用实际光谱反射率值去除包络线上相应波段反射率值,其特点就是将反射率光谱归一化到0-1之间的值;经过连续统去除后,原先峰值点上对应的反射率值均变为1,非峰值点上对应的反射率值则小于1;通过连续统去除可以间接得到波段深度,即每一样点的土壤反射率归一化到对应的光谱背景上,有利于光谱曲线之间特征波段的比较,根据和BD=1-RCR(λ)进行计算,其中RCR(λ)为波长λ处连续统去除值,R(λ)是经过SG滤波后波长λ处的初始光谱值,REL(λ)是对应波长λ的去包络线值,BD是波长处对应的波段深度;
所述步骤3)中利用不同的建模方法对反射率数据与重金属含量进行回归方程拟合,构建出高光谱估算模型的详细步骤包括:
C1)利用统计回归分析法估算土壤重金属含量,其核心是以回归方程为理论基础,通过多种数学方法建立土壤重金属元素与对应的高光谱数据之间的联系,一元回归分析方法计算较为简单,即在土壤高光谱数据中找出与土壤重金属含量相关系数绝对值最高的波段作为自变量,对应的土壤重金属含量作为因变量,之后通过一定的数学公式进行回归方程的拟合,具体有一元线性函数、幂函数、指数函数、对数函数、逆函数;多元逐步回归分析基本思路为先依据各自变量的重要性,之后每一步选择一个重要的变量进入回归方程;多元逐步回归不仅考虑到按贡献度大小逐一挑选重要变量,而且还考虑到较早进入回归方程的某些变量,随着其后一些变量的选入而失去原有的重要性,这样的变量应及时从回归方程中被剔除,从而使得回归方程中始终只保留重要的变量;
C2)利用主成分回归法估算土壤重金属含量,主成分回归法是通过线性变换,将原来多个指标组合成相互独立的少数几个能充分反映总体信息的主成分,从而在不丢掉重要信息的前提下,避开变量间的共线性,偏于进一步分析,其核心思想为降维;
C3)利用偏最小二乘法估算土壤重金属含量,偏最小二乘法的算法基础为最小二乘法,因为其尽可能提取包含自变量更多信息的基础上,保证了提取成分与因变量间最大相关性,在使用偏最小二乘法建模时,需要正确选取主因子数量,因为主因子数量的多少与模型最终的实际预测能力有显著关系,在本发明中,所采用的确定最佳主因子数量的方法为留一交互验证法,该方法的基本原理为若样本集中有N个样本,每次从中剔除一个样本,用剩下的N-1个样本作为训练集拟合模型,然后用这个拟合模型预测被剔除样本的组分含量,之后,用同样的方法对样本集N中每个样本都剔除一遍,并且得到它们的预测值;通过这一方法进行预测产生的误差和称为PRESS,根据计算得到,其中是使用留一交叉验证法计算得到的第i个样本的预测值,yi是第i个样本的实测值,留一交叉验证法的均方根误差根据计算得到,通常情况下,当RMSECV取得最小值时,可以确定模型的最佳主因子数;
C4)利用地理加权回归法估算土壤重金属含量,地理加权回归法是由第一作者Fotheringham提出的一个空间变参数模型,该模型是一种对不同空间子区域自变量和因变量之间关系随空间变化进行建模的非参数局部空间回归分析方法;其核心思路是采样点对其附近点特征的影响要大于离其较远的点,点位的回归系数不再是利用全局信息获得,而是利用邻近观测值的样本数据进行局部回归估计得到,且将数据的空间位置嵌入到回归参数中,使得对空间估值问题更合乎实际;
所述步骤3)中利依据建模决定系数、验证决定系数、均方根误差和相对分析误差四个评价指标对每种重金属元素下的模型进行评价,得到每种重金属元素最佳的估算模型的详细步骤包括:
D1)采用决定系数评价模型的精度,决定系数用于评价模型建模与验证结果与1:1趋势线的拟合程度,并衡量模型精度的稳定性和准确度;决定系数的大小决定了相关的密切程度,当R2越接近于1,表示模型的效果越好,相反,越接近于0,则表示模型的效果越差,根据计算得到决定系数,其中xi是样本实测值,Xi为样本预测值,为样本平均值,n为样本数量;
D2)采用均方根误差评价模型的反演能力,均方根误差主要评价模型的反演能力,其值越小表明实测值与预测值之间的偏差越小,两者的一致性也越好,也就是模型的模拟结果越精确可靠,根据计算得到均方根误差,其中xi是样本实测值,Xi为样本预测值,n为样本数量;
D3)采用相对分析误差检验模型的预报能力,相对分析误差是验证样本实测值的标准差与验证样本预测值均方根误差的商,根据计算得到,第一作者Viscarra对RPD进行了标准的划分,当PRD<1.4时,模型无法对样本进行预测;当1.4≤RPD<2时,模型效果一般,可以用此模型对样本进行粗略评估;但RPD≥2时,模型具有极好的预测能力。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111309905.3A CN114018833B (zh) | 2021-11-07 | 2021-11-07 | 基于高光谱遥感技术估算土壤重金属含量的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111309905.3A CN114018833B (zh) | 2021-11-07 | 2021-11-07 | 基于高光谱遥感技术估算土壤重金属含量的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114018833A CN114018833A (zh) | 2022-02-08 |
CN114018833B true CN114018833B (zh) | 2023-12-19 |
Family
ID=80061984
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202111309905.3A Active CN114018833B (zh) | 2021-11-07 | 2021-11-07 | 基于高光谱遥感技术估算土壤重金属含量的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114018833B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114486786A (zh) * | 2022-03-03 | 2022-05-13 | 上海园林绿化建设有限公司 | 土壤有机质测定方法及测定系统 |
CN114814167B (zh) * | 2022-04-08 | 2023-02-03 | 中科吉安生态环境研究院 | 融合多源环境变量与光谱信息的土壤重金属含量反演方法 |
CN114660105A (zh) * | 2022-04-08 | 2022-06-24 | 成都理工大学 | 土壤重金属Cd含量反演方法、系统、介质、计算机设备 |
CN115656074B (zh) * | 2022-12-28 | 2023-04-07 | 山东省科学院海洋仪器仪表研究所 | 一种海水cod光谱变量特征自适应选择估计方法 |
CN116304957A (zh) * | 2023-05-17 | 2023-06-23 | 成都交大光芒科技股份有限公司 | 一种供变电设备监测状态突变在线识别方法 |
CN117233114B (zh) * | 2023-11-07 | 2024-01-30 | 吉林大学 | 一种基于多源数据融合的土壤养分自动检测装置及方法 |
CN117647501B (zh) * | 2024-01-29 | 2024-04-26 | 中国计量科学研究院 | 基于icp光谱特征的食品中元素智能检测方法及系统 |
CN118090648B (zh) * | 2024-04-02 | 2024-06-21 | 常州建昊建筑鉴定检测有限公司 | 一种工程材料含水率检测方法及系统 |
CN118329807B (zh) * | 2024-06-14 | 2024-08-23 | 辽宁新华龙大有钼业有限公司 | 一种钼铁合金的成分含量精确测量方法及系统 |
Citations (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102854238A (zh) * | 2012-08-10 | 2013-01-02 | 同济大学 | 一种环境介质或生物样品中多种重金属的测定分析方法 |
CN103544550A (zh) * | 2013-11-08 | 2014-01-29 | 湖南科技大学 | 一种金属矿区土-水界面重金属污染负荷的预测方法 |
CN106680461A (zh) * | 2016-12-29 | 2017-05-17 | 贵州省土壤肥料研究所 | 一种测定贵州省红壤铅含量的方法及其校正因子 |
CN106918565A (zh) * | 2017-03-02 | 2017-07-04 | 中南大学 | 基于室内标样高光谱特征的土壤重金属Cd含量反演建模及其光谱响应特征波段识别方法 |
CN108960622A (zh) * | 2018-06-28 | 2018-12-07 | 河南理工大学 | 一种基于遥感影像的矿区复垦土地质量的自动化评价方法 |
CN110174359A (zh) * | 2019-05-27 | 2019-08-27 | 生态环境部南京环境科学研究所 | 一种基于高斯过程回归的航空高光谱影像土壤重金属浓度评估方法 |
CN110222475A (zh) * | 2019-07-03 | 2019-09-10 | 中国水利水电科学研究院 | 一种基于无人机多光谱遥感反演冬小麦植株含水率的方法 |
CN110596028A (zh) * | 2019-10-22 | 2019-12-20 | 中国地质科学院矿产综合利用研究所 | 一种沉积型稀土La元素含量的高光谱反演方法 |
CN110793923A (zh) * | 2019-10-31 | 2020-02-14 | 北京绿土科技有限公司 | 基于手机的高光谱土壤数据采集与分析方法 |
CN111062251A (zh) * | 2020-03-23 | 2020-04-24 | 乔红波 | 基于无人机成像的农田棉蚜为害等级模型的监测方法 |
CN112229817A (zh) * | 2020-09-29 | 2021-01-15 | 东北大学 | 一种苏打盐碱地重金属定量反演模型建立方法 |
CN112525829A (zh) * | 2020-04-30 | 2021-03-19 | 中国科学院地球化学研究所 | 一种重金属含量检测设备 |
-
2021
- 2021-11-07 CN CN202111309905.3A patent/CN114018833B/zh active Active
Patent Citations (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102854238A (zh) * | 2012-08-10 | 2013-01-02 | 同济大学 | 一种环境介质或生物样品中多种重金属的测定分析方法 |
CN103544550A (zh) * | 2013-11-08 | 2014-01-29 | 湖南科技大学 | 一种金属矿区土-水界面重金属污染负荷的预测方法 |
CN106680461A (zh) * | 2016-12-29 | 2017-05-17 | 贵州省土壤肥料研究所 | 一种测定贵州省红壤铅含量的方法及其校正因子 |
CN106918565A (zh) * | 2017-03-02 | 2017-07-04 | 中南大学 | 基于室内标样高光谱特征的土壤重金属Cd含量反演建模及其光谱响应特征波段识别方法 |
CN108960622A (zh) * | 2018-06-28 | 2018-12-07 | 河南理工大学 | 一种基于遥感影像的矿区复垦土地质量的自动化评价方法 |
CN110174359A (zh) * | 2019-05-27 | 2019-08-27 | 生态环境部南京环境科学研究所 | 一种基于高斯过程回归的航空高光谱影像土壤重金属浓度评估方法 |
CN110222475A (zh) * | 2019-07-03 | 2019-09-10 | 中国水利水电科学研究院 | 一种基于无人机多光谱遥感反演冬小麦植株含水率的方法 |
CN110596028A (zh) * | 2019-10-22 | 2019-12-20 | 中国地质科学院矿产综合利用研究所 | 一种沉积型稀土La元素含量的高光谱反演方法 |
CN110793923A (zh) * | 2019-10-31 | 2020-02-14 | 北京绿土科技有限公司 | 基于手机的高光谱土壤数据采集与分析方法 |
CN111062251A (zh) * | 2020-03-23 | 2020-04-24 | 乔红波 | 基于无人机成像的农田棉蚜为害等级模型的监测方法 |
CN112525829A (zh) * | 2020-04-30 | 2021-03-19 | 中国科学院地球化学研究所 | 一种重金属含量检测设备 |
CN112229817A (zh) * | 2020-09-29 | 2021-01-15 | 东北大学 | 一种苏打盐碱地重金属定量反演模型建立方法 |
Non-Patent Citations (3)
Title |
---|
张立福等.《高光谱遥感信息处理》.湖北科学技术出版社,2021,(第1版),第241-243页. * |
李艳坤等.《模式识别方法及其在复杂体系中的应用》.燕山大学出版社,2020,(第1版),第28-31页. * |
祝修高.平潭海岛土壤元素的高光谱遥感反演模型.《中国优秀硕博论文工程科技Ⅰ辑》.2018,(第6期),第4-18、45-48页. * |
Also Published As
Publication number | Publication date |
---|---|
CN114018833A (zh) | 2022-02-08 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN114018833B (zh) | 基于高光谱遥感技术估算土壤重金属含量的方法 | |
Hong et al. | Estimating lead and zinc concentrations in peri-urban agricultural soils through reflectance spectroscopy: Effects of fractional-order derivative and random forest | |
Zhou et al. | Hyperspectral inversion of soil heavy metals in Three-River Source Region based on random forest model | |
Hong et al. | Rapid identification of soil organic matter level via visible and near-infrared spectroscopy: Effects of two-dimensional correlation coefficient and extreme learning machine | |
CN104897592B (zh) | 基于高光谱技术的盐渍化土壤盐分离子含量监测方法 | |
Liu et al. | Estimation of soil organic matter content based on CARS algorithm coupled with random forest | |
Chen et al. | Development of a soil heavy metal estimation method based on a spectral index: Combining fractional-order derivative pretreatment and the absorption mechanism | |
Liu et al. | Study on the prediction of soil heavy metal elements content based on visible near-infrared spectroscopy | |
Lu et al. | Rapid inversion of heavy metal concentration in karst grain producing areas based on hyperspectral bands associated with soil components | |
Wan et al. | Estimation of soil pH using PXRF spectrometry and Vis-NIR spectroscopy for rapid environmental risk assessment of soil heavy metals | |
CN110596028B (zh) | 一种沉积型稀土La元素含量的高光谱反演方法 | |
CN113436153B (zh) | 一种基于高光谱成像和支持向量机技术的原状土壤剖面碳组分预测方法 | |
Sun et al. | Specific inherent optical quantities of complex turbid inland waters, from the perspective of water classification | |
Xia et al. | Moisture spectral characteristics and hyperspectral inversion of fly ash-filled reconstructed soil | |
CN113866102A (zh) | 一种基于光谱的土壤健康调查监测方法 | |
Xiong et al. | The total P estimation with hyper-spectrum–a novel insight into different P fractions | |
CN114814167B (zh) | 融合多源环境变量与光谱信息的土壤重金属含量反演方法 | |
Mao et al. | Research on the quantitative inversion model of heavy metals in soda saline land based on visible-near-infrared spectroscopy | |
Yang et al. | Prediction of soil heavy metal concentrations in copper tailings area using hyperspectral reflectance | |
CN116578851A (zh) | 一种高光谱土壤有效硼含量预测方法 | |
Li et al. | Optical imaging spectroscopy coupled with machine learning for detecting heavy metal of plants: A review | |
CN113971989A (zh) | 一种基于opls的森林土壤有机碳含量高光谱建模方法 | |
Zhang et al. | On the parsimony, interpretability and predictive capability of a physically− based model in the optical domain for estimating soil moisture content | |
Zheng et al. | Evolution of paddy soil fertility in a millennium chronosequence based on imaging spectroscopy | |
CN111595806A (zh) | 一种利用中红外漫反射光谱监测土壤碳组分的方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |