CN111337434A - 一种矿区复垦植被生物量估算方法及系统 - Google Patents
一种矿区复垦植被生物量估算方法及系统 Download PDFInfo
- Publication number
- CN111337434A CN111337434A CN202010151513.8A CN202010151513A CN111337434A CN 111337434 A CN111337434 A CN 111337434A CN 202010151513 A CN202010151513 A CN 202010151513A CN 111337434 A CN111337434 A CN 111337434A
- Authority
- CN
- China
- Prior art keywords
- image
- remote sensing
- radar
- biomass
- fusion
- 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.)
- Pending
Links
- 239000002028 Biomass Substances 0.000 title claims abstract description 213
- 238000000034 method Methods 0.000 title claims abstract description 49
- 238000005065 mining Methods 0.000 title claims abstract description 31
- 230000003287 optical effect Effects 0.000 claims abstract description 178
- 230000004927 fusion Effects 0.000 claims abstract description 106
- 238000007781 pre-processing Methods 0.000 claims abstract description 40
- 238000007499 fusion processing Methods 0.000 claims abstract description 10
- 230000009466 transformation Effects 0.000 claims description 28
- 238000012937 correction Methods 0.000 claims description 19
- 238000004422 calculation algorithm Methods 0.000 claims description 13
- 230000001629 suppression Effects 0.000 claims description 10
- 230000001419 dependent effect Effects 0.000 claims description 9
- 238000000605 extraction Methods 0.000 claims description 8
- 238000010219 correlation analysis Methods 0.000 claims description 7
- 230000008901 benefit Effects 0.000 abstract description 5
- 238000012545 processing Methods 0.000 abstract description 2
- 238000005070 sampling Methods 0.000 description 16
- 230000003595 spectral effect Effects 0.000 description 15
- 230000010287 polarization Effects 0.000 description 14
- 238000002310 reflectometry Methods 0.000 description 7
- 239000011159 matrix material Substances 0.000 description 5
- 239000013598 vector Substances 0.000 description 5
- 230000005540 biological transmission Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 238000013507 mapping Methods 0.000 description 4
- 238000012544 monitoring process Methods 0.000 description 4
- 238000000513 principal component analysis Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000005855 radiation Effects 0.000 description 3
- 241000511730 Leymus chinensis Species 0.000 description 2
- DZSYJVXGONVNKA-UHFFFAOYSA-L NIR-1 dye Chemical compound [K+].[K+].C1=CC2=C(S([O-])(=O)=O)C=C(S([O-])(=O)=O)C=C2C(C2(C)C)=C1[N+](CC)=C2C=CC=CC=CC=C1C(C)(C)C2=CC(C(O)=O)=CC=C2N1CCCCS([O-])(=O)=O DZSYJVXGONVNKA-UHFFFAOYSA-L 0.000 description 2
- 238000012952 Resampling Methods 0.000 description 2
- 230000002146 bilateral effect Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 239000005416 organic matter Substances 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 235000003826 Artemisia Nutrition 0.000 description 1
- 235000003261 Artemisia vulgaris Nutrition 0.000 description 1
- 241000435900 Aster altaicus Species 0.000 description 1
- 241000045403 Astragalus propinquus Species 0.000 description 1
- 241000196324 Embryophyta Species 0.000 description 1
- 241000756137 Hemerocallis Species 0.000 description 1
- 101000583150 Homo sapiens Membrane-associated phosphatidylinositol transfer protein 3 Proteins 0.000 description 1
- 241000219823 Medicago Species 0.000 description 1
- 235000017587 Medicago sativa ssp. sativa Nutrition 0.000 description 1
- 102100030351 Membrane-associated phosphatidylinositol transfer protein 3 Human genes 0.000 description 1
- 241001668545 Pascopyrum Species 0.000 description 1
- 241001070186 Salsola collina Species 0.000 description 1
- 244000030166 artemisia Species 0.000 description 1
- 235000009052 artemisia Nutrition 0.000 description 1
- 235000006533 astragalus Nutrition 0.000 description 1
- 230000015556 catabolic process Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 229930002875 chlorophyll Natural products 0.000 description 1
- 235000019804 chlorophyll Nutrition 0.000 description 1
- ATNHDLDRLWWWCB-AENOIHSZSA-M chlorophyll a Chemical compound C1([C@@H](C(=O)OC)C(=O)C2=C3C)=C2N2C3=CC(C(CC)=C3C)=[N+]4C3=CC3=C(C=C)C(C)=C5N3[Mg-2]42[N+]2=C1[C@@H](CCC(=O)OC\C=C(/C)CCC[C@H](C)CCC[C@H](C)CCCC(C)C)[C@H](C)C2=C5 ATNHDLDRLWWWCB-AENOIHSZSA-M 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000006731 degradation reaction Methods 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 238000001035 drying Methods 0.000 description 1
- 230000005684 electric field Effects 0.000 description 1
- 230000035558 fertility Effects 0.000 description 1
- 238000003306 harvesting Methods 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000012535 impurity Substances 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 239000003595 mist Substances 0.000 description 1
- 238000007500 overflow downdraw method Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 241000894007 species Species 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
- 239000002699 waste material Substances 0.000 description 1
- 238000005303 weighing Methods 0.000 description 1
Images
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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/86—Combinations of radar systems with non-radar systems, e.g. sonar, direction finder
- G01S13/867—Combination of radar systems with cameras
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/25—Fusion techniques
- G06F18/253—Fusion techniques of extracted features
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
- G06V10/75—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
- G06V10/751—Comparing pixel values or logical combinations thereof, or feature values having positional relevance, e.g. template matching
-
- 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
- G01N2021/1793—Remote sensing
- G01N2021/1797—Remote sensing in landscape, e.g. crops
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Computer Networks & Wireless Communication (AREA)
- Evolutionary Computation (AREA)
- Artificial Intelligence (AREA)
- General Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Data Mining & Analysis (AREA)
- Health & Medical Sciences (AREA)
- Computing Systems (AREA)
- Analytical Chemistry (AREA)
- General Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Multimedia (AREA)
- Chemical & Material Sciences (AREA)
- Databases & Information Systems (AREA)
- Biochemistry (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Electromagnetism (AREA)
- Image Processing (AREA)
Abstract
本发明公开一种矿区复垦植被生物量估算方法及系统，涉及植被估算技术领域。该方法包括：获取地面的历史生物量数据与历史遥感数据；历史遥感数据包括：光学遥感数据和雷达遥感数据；对历史遥感数据依次进行预处理和融合处理得到融合图像；融合图像的遥感特征变量为融合图像的融合波段，建立生物量估算模型；对获取的待测地面的遥感数据进行处理得到待测融合图像；待测融合图像的待测遥感特征变量为待测融合图像的融合波段；将待测遥感特征变量输入生物量估算模型得到待测地面的生物量数据。本发明对历史遥感数据进行预处理后再进行融合处理，将光学遥感数据和雷达遥感数据的优势进行融合互补，以此提高植被生物量的估算精度。
Description
技术领域
本发明涉及植被估算技术领域，特别是涉及一种矿区复垦植被生物量估算方法及系统。
背景技术
矿区大量的开采活动会极大地破坏地表景观，导致土地肥力下降及土地沉陷，造成矿区生态系统的严重退化。植被的重建与复垦是矿区生态环境修复的主要技术之一，植被恢复情况是评价矿区生态环境修复的直观指标，因此，复垦植被生长状况的准确监测，有利于矿区生态环境的恢复与管理。常用于监测植被生长状况的参数包括植被覆盖度、叶面积指数、冠层叶绿素含量以及生物量等。植被生物量是指植物在一定数量的面积上所能积累的有机物的总量，可直接反映植被健康状况，体现了矿区生态系统获取能量的能力，对于衡量生态系统尤为重要。生物量的准确估算能直观地评价植被复垦效果。而在长期的植被复垦过程中，随着时间的推移植被物种会出现更替现象，土壤状况也会发生变化，说明植被复垦工作具有持续性和动态性，这就更加突出持续有效的监测复垦植被恢复信息对矿区复垦工程的重要性与现实意义。
但现有的生物量估算方法耗时长，需要大量人力物力财力以及部分采集操作对植被具有破坏性，所以对于大面积区域的生物量估算，具有一定的局限性。
发明内容
本发明的目的是提供一种矿区复垦植被生物量估算方法及系统，解决了现有生物量估算方法对于大面积区域生物量估算的局限问题。
为实现上述目的，本发明提供了如下方案：
一种矿区复垦植被生物量估算方法，包括：
获取地面的历史生物量数据以及与所述历史生物量数据对应的历史遥感数据；所述历史遥感数据包括：光学遥感数据和雷达遥感数据；
对所述历史遥感数据依次进行预处理和融合处理，得到融合图像；
获取所述融合图像的融合波段，并将所述融合波段作为所述融合图像的遥感特征变量；
利用回归模型、所述遥感特征变量和所述生物量数据建立生物量估算模型；
获取待测地面的遥感数据；
对所述遥感数据依次进行预处理和融合处理，得到待测融合图像；
获取所述待测融合图像的融合波段，并将所述待测融合图像的融合波段作为所述待测融合图像的待测遥感特征变量；
将所述待测遥感特征变量输入所述生物量估算模型，得到所述待测地面的生物量数据。
可选的，所述对所述历史遥感数据依次进行预处理和融合处理，得到融合图像，具体包括：
对所述光学遥感数据进行预处理，得到光学图像；
对所述雷达遥感数据进行预处理，得到雷达图像；
对所述光学图像和所述雷达图像进行图像配准，得到与所述光学图像对应的配准雷达图像；
利用小波主成分算法将所述光学图像和与所述光学图像对应的配准雷达图像进行图像融合，得到融合图像。
可选的，所述对所述光学遥感数据进行预处理，得到光学图像，具体包括：
对所述光学遥感数据依次进行辐射定标和大气校正，得到光学图像；
所述对所述雷达遥感数据进行预处理，得到雷达图像，具体包括：
对所述雷达遥感数据依次进行辐射定标、斑点噪声抑制、地形校正和超分辨率重建，得到雷达图像；
所述对所述光学图像和所述雷达图像进行图像配准，得到与所述光学图像对应的配准雷达图像，具体包括：
以所述光学图像为标准参考图像，以所述雷达图像为待配准图像，选取所述光学图像和所述雷达图像上的地物特征点；
根据所述地物特征点，确定所述光学图像和所述雷达图像的非线性关系；
利用所述非线性关系将所述雷达图像匹配到所述光学图像的坐标系下，得到与所述光学图像对应的配准雷达图像。
可选的，所述利用小波主成分算法将所述光学图像和与所述光学图像对应的配准雷达图像进行图像融合，得到融合图像，具体包括：
对所述光学图像进行主成分正变换，得到所述光学图像的主成分；所述主成分包括：第一主成分和剩余主成分；
将所述光学图像的第一主成分和所述配准雷达图像进行直方图匹配，得到匹配第一主成分和匹配雷达图像；
对所述匹配第一主成分和所述匹配雷达图像分别进行小波变换，得到多个主成分小波分量和多个雷达小波分量；
将多个所述主成分小波分量和多个所述雷达小波分量进行重新组合，得到多个小波分量，并对多个所述小波分量进行小波逆变换，得到第一主成分图像；
对所述第一主成分图像和所述剩余主成分进行逆变换，得到融合图像。
可选的，所述利用回归模型、所述遥感特征变量和所述生物量数据建立生物量估算模型，具体包括：
对所述遥感特征变量和所述生物量数据进行相关性分析，得到相关性值大于预设相关性值的遥感特征变量因子；
将所述遥感特征变量因子作为自变量，所述生物量数据作为因变量，利用回归模型建模，得到建立的生物量估算模型。
一种矿区复垦植被生物量估算系统，包括：
获取模块，用于获取地面的历史生物量数据以及与所述历史生物量数据对应的历史遥感数据；所述历史遥感数据包括：光学遥感数据和雷达遥感数据；
融合图像模块，用于对所述历史遥感数据依次进行预处理和融合处理，得到融合图像；
提取模块，用于获取所述融合图像的融合波段，并将所述融合波段作为所述融合图像的遥感特征变量；
建立模块，用于利用回归模型、所述遥感特征变量和所述生物量数据建立生物量估算模型；
遥感数据模块，用于获取待测地面的遥感数据；
待测融合图像模块，用于对所述遥感数据依次进行预处理和融合处理，得到待测融合图像；
待测遥感特征变量模块，用于获取所述待测融合图像的融合波段，并将所述待测融合图像的融合波段作为所述待测融合图像的待测遥感特征变量；
生物量数据模块，用于将所述待测遥感特征变量输入所述生物量估算模型，得到所述待测地面的生物量数据。
可选的，所述融合图像模块，具体包括：
光学图像单元，用于对所述光学遥感数据进行预处理，得到光学图像；
雷达图像单元，用于对所述雷达遥感数据进行预处理，得到雷达图像；
匹配单元，用于对所述光学图像和所述雷达图像进行图像配准，得到与所述光学图像对应的配准雷达图像；
融合图像单元，用于利用小波主成分算法将所述光学图像和与所述光学图像对应的配准雷达图像进行图像融合，得到融合图像。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明提供了一种矿区复垦植被生物量估算方法及系统。该方法包括：获取地面的历史生物量数据以及与历史生物量数据对应的历史遥感数据；历史遥感数据包括：光学遥感数据和雷达遥感数据；对历史遥感数据依次进行预处理和融合处理，得到融合图像；获取融合图像的融合波段，并将融合波段作为融合图像的遥感特征变量；利用回归模型、遥感特征变量和生物量数据建立生物量估算模型；获取待测地面的遥感数据；对遥感数据依次进行预处理和融合处理，得到待测融合图像；获取待测融合图像的融合波段，并将待测融合图像的融合波段作为待测融合图像的待测遥感特征变量；将待测遥感特征变量输入生物量估算模型，得到待测地面的生物量数据。本发明对历史遥感数据进行预处理后再进行融合处理，可以将光学遥感数据和雷达遥感数据的优势进行融合互补，以此提高植被生物量的估算精度，光学遥感数据尤其是高分辨率影像能为植被监测提供丰富的光谱信息和纹理特征，但只能获得植被表面信息，而雷达遥感数据能够穿透云雾，可在一定程度穿透植被冠层，获得植被冠层以下的结构信息，因此，将光学遥感数据和雷达遥感数据相融合可以同时获得植被表面和植被冠层以下的信息，从而可以提高植被生物量的估算精度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例所提供的矿区复垦植被生物量估算方法的流程图；
图2为本发明实施例所提供的矿区复垦植被生物量估算系统的系统图；
图3为本发明实施例所提供的采样地区及采样点位置图；
图4为本发明实施例所提供的数据预处理的流程图；
图5为本发明实施例所提供的融合过程的流程图；
图6为本发明实施例所提供的融合结果图；图6(a)为WV-3图像，图6(b)为配准后的Sentinel-1 SAR图像，图6(c)为融合图像；
图7为本发明实施例所提供的WV-3图像、配准后的Sentinel-1 SAR图像和融合图像的信息熵的柱形图；
图8为本发明实施例所提供的WV-3图像、配准后的Sentinel-1 SAR图像和融合图像的平均梯度的柱形图；
图9为本发明实施例所提供的WV-3图像、配准后的Sentinel-1 SAR图像和融合图像的空间相关系数和光谱扭曲度的折线图；
图10为本发明实施例所提供的生物量估算模型残差分布图；图10(a)为EVI模型的模型残差分布图；图10(b)为VH模型的模型残差分布图；图10(c)为NDVI+VHME模型的模型残差分布图；图10(d)为RHB8模型的模型残差分布图；
图11为本发明实施例所提供的四种模型生物量估算值等级分布图；图11(a)为EVI模型的生物量估算值等级分布图；图11(b)为VH模型的生物量估算值等级分布图；图11(c)为NDVI+VHME模型的生物量估算值等级分布图；图11(d)为RHB8模型的生物量估算值等级分布图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种矿区复垦植被生物量估算方法及系统，解决了现有生物量估算方法对于大面积区域生物量估算的局限问题。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
本实施例提供一种矿区复垦植被生物量估算方法，图1为本发明实施例所提供的矿区复垦植被生物量估算方法的流程图。参见图1，矿区复垦植被生物量估算方法包括：
步骤101，获取地面的历史生物量数据以及与历史生物量数据对应的历史遥感数据；历史遥感数据包括：光学遥感数据和雷达遥感数据。历史生物量数据为历史实际测量的地面单位面积内的有机物质总量。
步骤102，对历史遥感数据依次进行预处理和融合处理，得到融合图像。受大气、地形以及遥感卫星自身的影响，历史遥感影像中的历史遥感数据存在误差，通过预处理消除误差。
步骤102具体包括：对光学遥感数据进行预处理，得到光学图像。具体包括：
对光学遥感数据依次进行辐射定标和大气校正，得到光学图像。
利用公式(1)将光学遥感数据中的原始遥感影像像元亮度值(Digital Number，DN值)转换为辐射亮度值，得到辐射定标后的光学遥感数据：
Le(λe)＝Gain·DN+Offset (1)
式中Le(λe)表示辐射亮度值，单位为W·m-2·sr-1·μm-1；DN表示DN值。Gain表示定标斜率，Offset表示定标截距，单位均为W·m-2·sr-1·μm-1。
辐射定标的原理是通过一定的转换规则从遥感影像的原始遥感影像像元亮度值(Digital Number，DN值)中获得大气表观反射率或辐射亮度值，辐射定标通常可分为相对定标和绝对定标。本实施例采用的是绝对定标。
利用FLAASH模型对辐射定标后的光学遥感数据进行大气校正，得到光学图像。大气校正的方法主要包括基于图形特征模型、地面线形回归经验模型以及大气辐射传输模型，目前，国内外学者已经提出了许多基于大气辐射传输模型的大气校正方法，其中得到广泛应用的主要有LOWTRAN模型、中光谱分辨率(Moderate Resolution Transmission，MODTRAN)模型、6S(Second Simulation of the Satellite Signal in theSolarSpectrum)模型、FLAASH(Fast Line-of-sight-Atmospheric Analysis of SpectralHypercubes)模型以及基于MODTRAN4+辐射传输模型改进的FLAASH模型等。
对雷达遥感数据进行预处理，得到雷达图像，具体包括：对雷达遥感数据依次进行辐射定标、斑点噪声抑制、地形校正和超分辨率重建，得到雷达图像。
辐射定标为：利用公式(2)和后向散射截面计算雷达遥感数据辐射定标后像元的强度值，得到辐射定标后的雷达遥感数据。
上式中，pwr表示辐射定标后像元的强度值，α表示目标地物点的局部入射角，αref表示目标地物点的参考入射角，κ表示定标系数，σ0表示后向散射截面的后向散射系数，κ的取值和σ0的计算方法根据欧洲航天局的文件确定。
采用基于3×3的滑动窗口的Lee滤波方法对辐射定标后的雷达遥感数据进行降噪处理：将滑动窗口内辐射定标后的雷达遥感数据的所有像素值均输入进采用Lee滤波算法的滤波器，获取滑动窗口中心点，即3×3滑动窗口中第二行第二列的点的像素值，重复上述操作遍历辐射定标后的雷达遥感数据，得到斑点噪声抑制后的雷达遥感数据。
将雷达遥感数据对应区域的数字高程模型(Digital Elevation Model，DEM)数据从大地经纬度坐标系转换到地心直角坐标系下，得到第一DEM数据；然后利用合成孔径雷达(Synthetic Aperture Radar，SAR)成像时刻各个方位向的卫星状态矢量，计算雷达遥感数据对应区域与SAR系统之间的斜距及多普勒频率；再利用斜距、多普勒频率、基于距离-多普勒(Range-Doppler，RD)算法的定位模型和第一DEM数据，得到第一DEM数据与预设的模拟SAR图像之间的映射关系；最后利用预设的模拟SAR图像与斑点噪声抑制后的雷达遥感数据之间所对应的纹理特征对预设的模拟SAR图像与斑点噪声抑制后的雷达遥感数据进行配准，以此将斑点噪声抑制后的雷达遥感数据映射到预设的模拟SAR图像，再通过第一DEM数据与预设的模拟SAR图像的映射关系，将斑点噪声抑制后的雷达遥感数据映射到第一DEM数据上，得到第二DEM数据，并将第二DEM数据转换到大地经纬度坐标系进行表示，得到地形校正后的雷达遥感数据。
采用凸集投影(projection-over-convex-set，POCS)算法对地形校正后的雷达遥感数据进行超分辨率重建，得到雷达图像，以此来提高雷达图像的空间分辨率，为和高空间分辨率的光学图像进行融合做准备。
对光学图像和雷达图像进行图像配准，得到与光学图像对应的配准雷达图像。具体包括：以光学图像为标准参考图像，以雷达图像为待配准图像，选取光学图像和雷达图像上的地物特征点。地物特征点包括：路口和房角等图像上具有明显地物特征的点。
根据地物特征点，确定光学图像和雷达图像的非线性关系。
利用非线性关系将雷达图像匹配到光学图像的坐标系下，并采用三次卷积(CubicConvolution)插值方法对匹配到光学图像的坐标系下雷达图像进行重采样，得到与光学图像对应的配准雷达图像。
利用小波主成分(W-PCA)算法将光学图像和与光学图像对应的配准雷达图像进行图像融合，得到融合图像。主成分分析方法(Principal Component Analysis，PCA)能够较好的保留多波段光谱信息，而小波变换能够较好的保留纹理结构信息，故将两者结合实现优势互补，即为W-PCA算法。具体包括：
对光学图像进行主成分正变换，得到光学图像的主成分；主成分包括：第一主成分和剩余主成分，剩余主成分包括：第二主成分、第三主成分、第四主成分、第五主成分、第六主成分、第七主成分和第八主成分。具体包括：1)计算光学图像X的均值向量M和协方差矩阵Σ；2)计算协方差矩阵Σ的特征值λr和特征向量
r表示主成分序号，r＝1,2,....,8；3)将特征向量
按照对应的特征值λr由大到小的次序进行排列，构造变换矩阵φ；4)根据公式Y＝φX对光学图像进行变换，得到主成分正变换后的光学图像Y和光学图像的主成分，主成分正变换后的光学图像Y包含8个主成分。
将光学图像的第一主成分和配准雷达图像进行直方图匹配，得到匹配第一主成分和匹配雷达图像，使匹配第一主成分和匹配雷达图像具有相近的均值和方差。与光学图像的第一主成分进行直方图匹配的配准雷达图像为与光学图像对应的配准雷达图像。
对匹配第一主成分和匹配雷达图像分别进行小波变换，得到多个主成分小波分量和多个雷达小波分量。
将多个主成分小波分量和多个雷达小波分量进行重新组合，得到多个小波分量，并对多个小波分量进行小波逆变换，得到第一主成分图像。
对第一主成分图像和剩余主成分进行逆变换，得到融合图像，融合图像为包括8个融合波段的图像。
步骤103，获取融合图像的融合波段，并将融合波段作为融合图像的遥感特征变量。
步骤104，利用回归模型、遥感特征变量和生物量数据建立生物量估算模型。
步骤104具体包括：对遥感特征变量和生物量数据进行相关性分析，得到相关性值大于预设相关性值的遥感特征变量因子。
将遥感特征变量因子作为自变量，生物量数据作为因变量，利用回归模型建模，得到建立的生物量估算模型。具体包括：
回归模型包括：y＝a+bx、y＝axb、y＝ax2+bx+c、y＝ax3+bx2+cx+d、y＝a+blnx和y＝aebx，其中y表示因变量，x表示自变量，a、b、c、d表示回归模型的待求参数。
将遥感特征变量因子作为自变量，生物量数据作为因变量，分别代入回归模型，计算得到a、b、c、d，进而得到建立的6个初步生物量估算模型。
根据公式(2)和(3)计算6个初步生物量估算模型的均方根误差和估算精度：
上式中，RMSE表示均方根误差，AGRn表示实际测量的生物量数据，
表示初步生物量估算模型估算的生物量数据，N表示实际测量的生物量数据的样本数量，Ac表示估算精度，
表示实际测量的生物量数据的平均值，n表示实际测量的生物量数据的样本序号，n＝1,2,...,N。
比较6个初步生物量估算模型的均方根误差和估算精度，确定均方根误差低和估算精度高的1个初步生物量估算模型为生物量估算模型。均方根误差值越低，表示建立的模型越好；估算精度值越高，表示建立的模型越好。
步骤105，获取待测地面的遥感数据。遥感数据包括：待测光学遥感数据和待测雷达遥感数据。
步骤106，对遥感数据依次进行预处理和融合处理，得到待测融合图像。具体为根据步骤102的预处理和融合处理过程对遥感数据依次进行预处理和融合处理，得到待测融合图像。待测融合图像为包括8个融合波段的图像。
步骤107，获取待测融合图像的融合波段，并将待测融合图像的融合波段作为待测融合图像的待测遥感特征变量。
步骤108，将待测遥感特征变量输入生物量估算模型，得到待测地面的生物量数据。
本实施例提供一种矿区复垦植被生物量估算系统，图2为本发明实施例所提供的矿区复垦植被生物量估算系统的系统图。参见图2，矿区复垦植被生物量估算系统包括：获取模块201，用于获取地面的历史生物量数据以及与历史生物量数据对应的历史遥感数据；历史遥感数据包括：光学遥感数据和雷达遥感数据。
融合图像模块202，用于对历史遥感数据依次进行预处理和融合处理，得到融合图像。
融合图像模块202具体包括：光学图像单元，用于对光学遥感数据进行预处理，得到光学图像。
光学图像单元具体包括：光学图像子单元，用于对光学遥感数据依次进行辐射定标和大气校正，得到光学图像。
雷达图像单元，用于对雷达遥感数据进行预处理，得到雷达图像。
雷达图像单元具体包括：雷达图像子单元，用于对雷达遥感数据依次进行辐射定标、斑点噪声抑制、地形校正和超分辨率重建，得到雷达图像。
匹配单元，用于对光学图像和雷达图像进行图像配准，得到与光学图像对应的配准雷达图像。
匹配单元具体包括：地物特征点子单元，用于以光学图像为标准参考图像，以雷达图像为待配准图像，选取光学图像和雷达图像上的地物特征点。
非线性关系子单元，用于根据地物特征点，确定光学图像和雷达图像的非线性关系。
匹配子单元，用于利用非线性关系将雷达图像匹配到光学图像的坐标系下，并采用三次卷积(Cubic Convolution)插值方法对匹配到光学图像的坐标系下雷达图像进行重采样，得到与光学图像对应的配准雷达图像。
融合图像单元，用于利用小波主成分(W-PCA)算法将光学图像和与光学图像对应的配准雷达图像进行图像融合，得到融合图像。
融合图像单元具体包括：主成分子单元，用于对光学图像进行主成分正变换，得到光学图像的主成分；主成分包括：第一主成分和剩余主成分。
直方图匹配子单元，用于将光学图像的第一主成分和配准雷达图像进行直方图匹配，得到匹配第一主成分和匹配雷达图像，使匹配第一主成分和匹配雷达图像具有相近的均值和方差。与光学图像的第一主成分进行直方图匹配的配准雷达图像为与光学图像对应的配准雷达图像。
小波变换子单元，用于对匹配第一主成分和匹配雷达图像分别进行小波变换，得到多个主成分小波分量和多个雷达小波分量。
小波逆变换子单元，用于将多个主成分小波分量和多个雷达小波分量进行重新组合，得到多个小波分量，并对多个小波分量进行小波逆变换，得到第一主成分图像。
融合图像子单元，用于对第一主成分图像和剩余主成分进行逆变换，得到融合图像，融合图像为包括8个融合波段的图像。
提取模块203，用于获取融合图像的融合波段，并将融合波段作为融合图像的遥感特征变量。
建立模块204，用于利用回归模型、遥感特征变量和生物量数据建立生物量估算模型。
建立模块204具体包括：遥感特征变量因子单元，用于对遥感特征变量和生物量数据进行相关性分析，得到相关性值大于预设相关性值的遥感特征变量因子。
生物量估算模型单元，用于将遥感特征变量因子作为自变量，生物量数据作为因变量，利用回归模型建模，得到建立的生物量估算模型。
遥感数据模块205，用于获取待测地面的遥感数据。遥感数据包括：待测光学遥感数据和待测雷达遥感数据。
待测融合图像模块206，用于对遥感数据依次进行预处理和融合处理，得到待测融合图像。待测融合图像为包括8个融合波段的图像。
待测遥感特征变量模块207，用于获取待测融合图像的融合波段，并将待测融合图像的融合波段作为待测融合图像的待测遥感特征变量。
生物量数据模块208，用于将待测遥感特征变量输入生物量估算模型，得到待测地面的生物量数据。
本实施例以内蒙古呼伦贝尔宝日希勒露天矿排土场为采样地区，验证本发明矿区复垦植被生物量估算方法的估算精度。内蒙古呼伦贝尔宝日希勒露天矿排土场采样地区的植被类型与生长状况具有多样性，下垫面较为复杂，植被类型分布有披碱草、大籽蒿、黄花草苜蓿、阿尔泰狗娃花、羊草、冰草、柴胡、黄芪和猪毛菜等。
步骤1：地面实测生物量数据获取
获取该采样地区地面生物量的样方数据，在整个采样地区中，随机设置多个采样点，以采样点为中心确定1m×1m的小样方，需尽可能保证采样点分布覆盖全部采样地区，参见图3。利用全球定位系统(Global Positioning System，GPS)记录每个采样点的地理坐标，并记录小样方中的植被类型、植被总盖度和复垦年限；然后将1m×1m小样方内的植被全部齐地面收割，去除碎石等杂物，装入袋子做好标记，最后共获得32个小样方。将装有植被的袋子带回实验室，将袋子中的植被在80℃下烘干至恒重，然后用精度为0.01的电子天平称重，得到单位面积内的有机物质总量，即实测的生物量数据。图3中，N表示北方；0～300和0～2分别表示所在图像的地图比例，Km表示千米，2009～2016表示复垦年度。
步骤2：遥感数据预处理
获取遥感数据，其中光学遥感数据为2016年8月11日的重采样为2m空间分辨率的WV-38波段多光谱图像；雷达遥感数据为2016年7月29日、2016年8月5日、2016年8月10日和2016年8月17日的分辨率为5×20m2的Sentinel-1 SAR双极化C波段图像。遥感数据参数参见表1和表2。
表1光学遥感数据参数
表2雷达遥感数据参数
雷达信号可以传送水平(H)或者垂直(V)的电场矢量，接收水平(H)或者垂直(V)或者两者的返回信号。雷达遥感系统常用的四种极化方式为HH、VV、HV、VH，前两者为同向极化，后两者为异向(交叉)极化。
受大气、地形以及遥感卫星的传感器自身的影响，获取的遥感影像存在误差，需要对遥感影像进行遥感定量分析消除误差，即对遥感数据进行数据预处理，参见图4：1)光学遥感数据预处理。对WV-38波段多光谱图像依次进行辐射定标和大气校正，得到光学图像WV-3图像。本实施例利用绝对定标来完成辐射校正。本实施例利用FLAASH模型对辐射定标后的光学遥感数据进行大气校正，得到光学图像WV-3图像。
2)雷达数据预处理。对4幅Sentinel-1 SAR双极化C波段图像依次进行辐射定标、斑点噪声抑制和地形校正后得到4幅双极化C波段图像，对4幅双极化C波段图像进行超分辨率重建得到1幅雷达图像Sentinel-1 SAR图像。
3)地理配准。将雷达图像Sentinel-1 SAR图像与光学图像WV-3图像进行配准，得到与WV-3图像配准后的Sentinel-1 SAR图像。
图像配准的过程是基于一些特定的匹配规则，将多幅基于不同传感器、不同时间以及不同角度获取的同一场景图像，转换到相同的坐标系下，并在像素级上获取误差最小的最佳匹配结果。本实施例以WV-3图像为标准参考图像，以Sentinel-1 SAR图像为待配准图像，通过选取地物特征点建立WV-3图像和Sentinel-1 SAR图像的非线性关系，采用CubicConvolution插值方法将Sentinel-1 SAR图像重采样到WV-3图像的坐标系下，得到与WV-3图像配准后的Sentinel-1 SAR图像。本实施例共选取30个地物特征点，且所有地物特征点的误差均小于0.5个像素，地物特征点的均方根误差为0.054个像素。
步骤3：光学与雷达图像融合方法
本实施例利用W-PCA算法将WV-3图像和配准后的Sentinel-1 SAR图像进行融合，融合过程参见图5：(1)对WV-3图像进行主成分正变换，得到WV-3图像的8个主成分，其中第一主成分含有各波段图像最多的有用信息，因此提取WV-3图像的第一主成分PC1。图5中PC1、PC2、...、PC8分别表示第一主成分、第二主成分、...、第八主成分。
(2)将配准后的Sentinel-1 SAR图像和第一主成分PC1进行直方图匹配，得到匹配第一主成分和匹配雷达图像，使匹配PC1和匹配Sentinel-1 SAR图像具有相近的均值和方差。直方图匹配又称为直方图规定化，是指将一幅图像的直方图变成规定形状的直方图而进行的图像增强方法，即将某幅影像或某一区域的直方图匹配到另一幅影像上，使两幅影像的色调保持一致。
(3)对匹配PC1和匹配Sentinel-1 SAR图像分别进行小波变换，将匹配PC1和匹配Sentinel-1 SAR图像均分解成LL(低频部分)、HL(水平方向的小波系数)、LH(垂直方向的小波系数)、HH(对角方向的小波系数)4种小波分量，得到4个主成分小波分量和4个雷达小波分量，将4个主成分小波分量和4个雷达小波分量对应重新组合成新的小波分量LL、HL、LH、HH，即将主成分小波分量LL’和雷达小波分量LL”组合成新的LL，主成分小波分量HL’和雷达小波分量HL”组合成新的HL，主成分小波分量LH’和雷达小波分量LH”组合成新的LH，主成分小波分量HH’和雷达小波分量HH”组合成新的HH；采用最大绝对值原则重新组合LL，采用加权平均原则重新组合HL、LH和HH，并对新的小波分量进行小波逆变换，生成PC1图像。图5中LL’、HL’、LH’、HH’分别表示主成分小波分量；LL”、HL”、LH”、HH”分别表示雷达小波分量。
(4)将携带着Sentinel-1SAR信息的PC1图像和(1)中的另外7个主成分通过主成分逆变换X'＝φ-1Y'得到融合后的包括8个融合波段的图像，即融合图像。上式中X’表示融合后的融合图像，φ-1表示变换矩阵的逆矩阵，Y’表示携带着Sentinel-1SAR信息的PC1图像和(1)中的另外7个主成分。
本实施例基于影像信息丰富程度、影像清晰度提高程度及影像光谱特征保真程度几个方面，选取信息熵、平均梯度、空间相关系数及光谱扭曲程度4个评价指标对融合图像进行评价。
信息熵代表融合图像包含的信息量，无过多噪声影响时，一般信息熵越丰富，融合图像信息越丰富，细节信息越好，根据下式计算信息熵H：
其中，H表示信息熵，M表示融合图像灰度级，m表示像元的灰度值，Pm代表灰度值为m的像元出现的频率，hm表示灰度值为m的像素数，h表示融合图像中像素总数。
平均梯度反映融合图像的清晰度，表现出融合图像中微小的细节变化及纹理反差特征，一般平均梯度值越大，融合图像视觉上越清晰，根据下式计算平均梯度
其中，
分别表示融合图像的像素值在横坐标方向、纵坐标方向的差分，P、Q分别表示融合图像的总行数、总列数，p、q分别表示融合图像的行数、列数；平均梯度越大，表示融合图像细节纹理信息越清晰，故平均梯度为地质复杂融合图像的重要指标。
空间相关系数反映融合图像与原始图像的空间相关性，其空间相关系数越大，表示原始图像的纹理特征保持越好。本实施例中原始图像包括WV-3图像和配准后的Sentinel-1 SAR图像。根据下式计算空间相关系数R：
其中，P、Q分别表示图像的总行数、总列数，p、q分别表示图像的行数、列数；fA(p,q)、fB(p,q)分别表示融合图像与原始图像在(p,q)处的灰度值，
分别表示融合图像与原始图像的平均灰度值。
光谱扭曲度(Twisting)反映融合图像相对于原始图像光谱信息的失真程度，其值越小，失真程度越低，说明光谱保真度越高，一般光谱扭曲度的值越大，融合图像的光谱变化越严重，根据下式计算光谱扭曲度DFA：
其中，P、Q分别表示图像的总行数、总列数，p、q分别表示图像的行数、列数；fA(p,q)、fB(p,q)分别表示融合图像与原始图像在(p,q)处的灰度值。
评价指标参见图7-9，较高的信息熵值反映了融合图像与WV-3图像相比具有更多的细节。平均梯度越高，说明融合图像的清晰度和纹理信息都优于配准后的Sentinel-1SAR图像。与WV-3图像相比，融合图像的7波段和8波段具有较高的空间相关系数和较低的光谱扭曲度表示融合图像的纹理特征更多，光谱变化更少。图7-9中B1、B2、B3、B4、B5、B6、B7、B8分别表示融合图像的8个波段图像。图6为本发明实施例所提供的融合结果图，图6(a)为WV-3图像，图6(b)为配准后的Sentinel-1 SAR图像，图6(c)为融合图像。参见图6，可见融合图像的清晰度以及光谱和结构纹理特征在视觉上优于单一的WV-3图像或配准后的Sentinel-1 SAR图像。
步骤4：遥感特征变量提取
将融合图像的8个融合波段作为遥感特征变量。
步骤5：生物量估算模型
将遥感特征变量因子作为自变量，生物量数据作为因变量，分别代入回归模型，计算得到a、b、c、d，进而得到建立的6个初步生物量估算模型。
计算6个初步生物量估算模型的均方根误差和估算精度，并比较6个初步生物量估算模型的均方根误差和估算精度，确定均方根误差低和估算精度高的1个初步生物量估算模型作为生物量估算模型。
步骤6：验证生物量估算模型的精度
为验证本实施例生物量估算模型的精度，提取光学图像和雷达图像的遥感特征变量，并根据光学图像和雷达图像的遥感特征变量建立单一光学生物量估算模型、单一雷达生物量估算模型以及光学雷达联合生物量估算模型，并对单一光学生物量估算模型、单一雷达生物量估算模型、光学雷达联合生物量估算模型以及本实施例建立的生物量估算模型进行精度分析。
首先，获取光学图像的波段反射率和雷达图像的后向散射系数。
基于光学图像，利用光学图像的波段反射率计算得到光学图像的植被指数。本实施例选取了7个植被指数，植被指数包括：归一化植被指数(NDVI)、差值植被指数(DVI)、比值指数(RVI)、归一化绿度指数(NDGI)、抗大气植被指数(ARVI)、增强型植被指数(EVI)和红边归一化指数(NDVI705)。利用光学图像的各波段反射率根据表3中的公式计算得到光学图像的植被指数。
表3 7种植被指数计算公式及其优点
表3中，R、G、B、NIR1、RE分别表示光学图像在红色波段、绿色波段、蓝色波段、近红外1波段、红色边缘波段的反射率。
基于光学图像和雷达图像，选取纹理值中最常用的基于灰度共生矩阵的8个纹理特征统计量：均值(Mean，ME)、方差(Variance，VA)、同质性(Homogeneity，HO)、对比度(Contrast，CO)、差异性(Dissimilarity，DI)、信息熵(Entropy，EN)、相关性(Correlation，CR)、二阶矩(Second Moment，SM)，分别计算光学图像的纹理特征和雷达图像的纹理特征。灰度共生矩阵是表示像素距离与像素角度的关系函数，通过计算图像中任意两个像素点灰度在某一既定方向和距离的相关性，来表现图像的方向、间隔以及变化幅度和速度的综合信息。本实施例利用3×3滑窗口提取纹理特征，8个纹理特征的计算公式参见表4。
表4纹理特征计算公式
表4中，L表示图像的像元数；i、j分别表示像元所在的行数、列数；Pi,j表示第i行第j列的像元值；ui、uj分别表示像元所在的行数的均值、列数的均值；
表示分别表示像元所在的行数的方差、列数的方差。
光学图像的遥感特征变量包括：光学图像的波段反射率、纹理特征和植被指数；雷达图像的遥感特征变量包括：雷达图像的后向散射系数(VV极化、VH极化)和纹理特征；光学雷达联合的遥感特征变量包括：光学图像的波段反射率、光学图像的纹理特征、植被指数、雷达图像的后向散射系数(VV极化和VH极化)和雷达图像的纹理特征。
分别对光学图像的遥感特征变量、雷达图像的遥感特征变量、光学雷达联合的遥感特征变量、融合图像的遥感特征变量与生物量进行相关性分析。表5表示显著性p＝0.01和p＝0.05时的遥感特征变量及相关性，可见植被指数和光学图像的纹理特征与生物量的相关性数值高于WV-3图像的波段反射率，配准后的Sentinel-1SAR图像的VV极化、VH极化和配准后的Sentinel-1SAR图像的纹理特征均与生物量有显著的相关性。融合波段与生物量的相关性均大大提高，且均呈现显著相关性。选取每个遥感特征变量的数据集中显著性高或相关性值大于预设相关性值的遥感特征变量因子建立相应的生物量估算模型。
表5生物量与各遥感特征变量相关性分析
**表示在显著性p＝0.01的水平(双侧)上显著相关；*表示在显著性p＝0.05的水平(双侧)上显著相关；WV-3图像的纹理特征的变量分别表示WV-3图像波段4-8对应的纹理特征；BO表示波段，下角标表示纹理特征，表5包括第4-8波段的部分纹理特征，如BO7VA表示第7波段的方差；配准后的Sentinel-1SAR图像的极化及纹理特征的变量分别表示两种极化方式(VV、VH)以及两种极化方式对应的纹理特征。
然后将上述选取的遥感特征变量因子作为自变量，实测的生物量数据作为因变量，采用6种回归模型，分别得到6个初步单一光学生物量估算模型、6个初步单一雷达生物量估算模型以及6个初步光学雷达联合生物量估算模型。计算6个初步单一光学生物量估算模型的均方根误差与估算精度，确定均方根误差低和估算精度高的1个初步单一光学生物量估算模型作为单一光学生物量估算模型。计算6个初步单一雷达生物量估算模型的均方根误差与估算精度，确定均方根误差低和估算精度高的1个初步单一雷达生物量估算模型作为单一雷达生物量估算模型。计算6个初步光学雷达联合生物量估算模型的均方根误差与估算精度，确定均方根误差低和估算精度高的1个初步光学雷达联合生物量估算模型作为光学雷达联合生物量估算模型。
对单一光学生物量估算模型、单一雷达生物量估算模型、光学雷达联合生物量估算模型以及本实施例建立的生物量估算模型进行精度分析。生物量估算模型的建模精度和验证精度如表6所示，单一变量建立的模型中单一光学生物量估算模型即EVI模型精度最高，建模精度：R2＝0.8163，RMSE＝17.3187g·m-2，Ac＝80.56％；验证精度：R2＝0.7098，RMSE＝24.2018g·m-2，Ac＝73.12％。其中R2表示决定系数。相对于单一变量建立的模型，光学雷达联合生物量估算模型即NDVI+VHME模型提高了建模精度，建模精度：R2＝0.8340，RMSE＝16.4646g·m-2，Ac＝81.52％。相比于单一光学生物量估算模型、单一雷达生物量估算模型和光学雷达联合生物量估算模型，融合图像的生物量估算模型即RHB8模型提高了验证精度：R2＝0.7983，RMSE＝22.8283g·m-2，Ac＝74.64％。建模精度还包括模型残差，模型残差为生物量估算模型估算的生物量与实测的生物量的差值。图10为本发明实施例所提供的生物量估算模型残差分布图；图10(a)为EVI模型的模型残差分布图；图10(b)为VH模型的模型残差分布图；图10(c)为NDVI+VHME模型的模型残差分布图；图10(d)为RHB8模型的模型残差分布图。图10表明，EVI模型和VH模型均在实测生物量低于50g·m-2的采样点处生物量明显被高估；在实测生物量高于120g·m-2的采样点处生物量明显被低估。NDVI+VHME模型明显降低了实测生物量低于50g·m-2和高于150g·m-2采样点处的模型残差，而融合图像的RHB8模型降低了实测生物量在50-100g·m-2采样点处的模型残差。估算的生物量大于实测的生物量时，生物量被高估，模型残差值体现为正值；估算的生物量小于实测的生物量时，生物量被低估，模型残差值体现为负值；四种模型全部存在生物量被高估或低估的现象，但NDVI+VHME模型和RHB8模型可以减弱被高估或低估的程度，提高了植被生物量估算的精度。
表6各生物量估算模型及其精度
表6中，y’表示模型估算的生物量；RHB8表示融合图像的第8波段。
图11(a)为EVI模型的生物量估算值等级分布图；图11(b)为VH模型的生物量估算值等级分布图；图11(c)为NDVI+VHME模型的生物量估算值等级分布图；图11(d)为本实施例融合图像的生物量估算模型即RHB8模型估算的生物量估算值；图例表示生物量估算值；柱形图的横坐标是生物量估算值，纵坐标是生物量估算值等级在采样地区出现的像元数。图11将利用生物量估算模型得到的估算结果以空间分布图的形式展现出来，图11中生物量估算值分布仅体现植被信息提取结果中有植被的区域，植被信息提取结果中无植被的区域以白色体现。根据实测生物量的统计分布情况将表6中的四种类型的生物量估算模型估算出的生物量估算值分为7个等级：0-40g·m-2，40-60g·m-2，60-80g·m-2，80-100g·m-2，100-150g·m-2，150-200g·m-2，>200g·m-2，结合植被信息提取结果，整个采样地区的生物量估算值等级分布图如图11所示，生物量估算值在150-200g·m-2和>200g·m-2区间主要分布在2013年和2015年复垦区域，生物量的均值在61g·m-2左右；2009年、2010年和2011年复垦植被生物量均值在30g·m-2左右；2012年及2016年复垦植被生物量均值最低，在22g·m-2左右。
遥感技术手段的时间周期短，地域覆盖全面，且具有快速、客观、经济的特点，能够实现对整个矿区的全方位、大范围快速监测，为矿区复垦植被生物量监测与评估提供准确可靠的数据基础。光学影像尤其是高分辨率影像可以为植被监测提供丰富的光谱信息和纹理特征，但只能获得植被表面信息而无法得到植被的垂直信息；微波遥感能够穿透云雾，具有多极化特征，反映地面地物的后向散射强度，可在一定程度穿透植被冠层，获得植被冠层以下的结构信息，却容易被植被下垫面所影响。因此，利用单一数据源，例如光学或者雷达，进行植被信息提取及参数反演，都具有一定的局限性。而将光学数据和雷达数据的协同使用，可以弥补单独使用某种数据源的不足，光学数据和微波遥感的结合可以提高植被生物量估算的精度，提高遥感影像的地物解译能力。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。
Claims (10)
1.一种矿区复垦植被生物量估算方法，其特征在于，包括：
获取地面的历史生物量数据以及与所述历史生物量数据对应的历史遥感数据；所述历史遥感数据包括：光学遥感数据和雷达遥感数据；
对所述历史遥感数据依次进行预处理和融合处理，得到融合图像；
获取所述融合图像的融合波段，并将所述融合波段作为所述融合图像的遥感特征变量；
利用回归模型、所述遥感特征变量和所述生物量数据建立生物量估算模型；
获取待测地面的遥感数据；
对所述遥感数据依次进行预处理和融合处理，得到待测融合图像；
获取所述待测融合图像的融合波段，并将所述待测融合图像的融合波段作为所述待测融合图像的待测遥感特征变量；
将所述待测遥感特征变量输入所述生物量估算模型，得到所述待测地面的生物量数据。
2.根据权利要求1所述的矿区复垦植被生物量估算方法，其特征在于，所述对所述历史遥感数据依次进行预处理和融合处理，得到融合图像，具体包括：
对所述光学遥感数据进行预处理，得到光学图像；
对所述雷达遥感数据进行预处理，得到雷达图像；
对所述光学图像和所述雷达图像进行图像配准，得到与所述光学图像对应的配准雷达图像；
利用小波主成分算法将所述光学图像和与所述光学图像对应的配准雷达图像进行图像融合，得到融合图像。
3.根据权利要求2所述的矿区复垦植被生物量估算方法，其特征在于，所述对所述光学遥感数据进行预处理，得到光学图像，具体包括：
对所述光学遥感数据依次进行辐射定标和大气校正，得到光学图像；
所述对所述雷达遥感数据进行预处理，得到雷达图像，具体包括：
对所述雷达遥感数据依次进行辐射定标、斑点噪声抑制、地形校正和超分辨率重建，得到雷达图像；
所述对所述光学图像和所述雷达图像进行图像配准，得到与所述光学图像对应的配准雷达图像，具体包括：
以所述光学图像为标准参考图像，以所述雷达图像为待配准图像，选取所述光学图像和所述雷达图像上的地物特征点；
根据所述地物特征点，确定所述光学图像和所述雷达图像的非线性关系；
利用所述非线性关系将所述雷达图像匹配到所述光学图像的坐标系下，得到与所述光学图像对应的配准雷达图像。
4.根据权利要求3所述的矿区复垦植被生物量估算方法，其特征在于，所述利用小波主成分算法将所述光学图像和与所述光学图像对应的配准雷达图像进行图像融合，得到融合图像，具体包括：
对所述光学图像进行主成分正变换，得到所述光学图像的主成分；所述主成分包括：第一主成分和剩余主成分；
将所述光学图像的第一主成分和所述配准雷达图像进行直方图匹配，得到匹配第一主成分和匹配雷达图像；
对所述匹配第一主成分和所述匹配雷达图像分别进行小波变换，得到多个主成分小波分量和多个雷达小波分量；
将多个所述主成分小波分量和多个所述雷达小波分量进行重新组合，得到多个小波分量，并对多个所述小波分量进行小波逆变换，得到第一主成分图像；
对所述第一主成分图像和所述剩余主成分进行逆变换，得到融合图像。
5.根据权利要求4所述的矿区复垦植被生物量估算方法，其特征在于，所述利用回归模型、所述遥感特征变量和所述生物量数据建立生物量估算模型，具体包括：
对所述遥感特征变量和所述生物量数据进行相关性分析，得到相关性值大于预设相关性值的遥感特征变量因子；
将所述遥感特征变量因子作为自变量，所述生物量数据作为因变量，利用回归模型建模，得到建立的生物量估算模型。
6.一种矿区复垦植被生物量估算系统，其特征在于，包括：
获取模块，用于获取地面的历史生物量数据以及与所述历史生物量数据对应的历史遥感数据；所述历史遥感数据包括：光学遥感数据和雷达遥感数据；
融合图像模块，用于对所述历史遥感数据依次进行预处理和融合处理，得到融合图像；
提取模块，用于获取所述融合图像的融合波段，并将所述融合波段作为所述融合图像的遥感特征变量；
建立模块，用于利用回归模型、所述遥感特征变量和所述生物量数据建立生物量估算模型；
遥感数据模块，用于获取待测地面的遥感数据；
待测融合图像模块，用于对所述遥感数据依次进行预处理和融合处理，得到待测融合图像；
待测遥感特征变量模块，用于获取所述待测融合图像的融合波段，并将所述待测融合图像的融合波段作为所述待测融合图像的待测遥感特征变量；
生物量数据模块，用于将所述待测遥感特征变量输入所述生物量估算模型，得到所述待测地面的生物量数据。
7.根据权利要求6所述的矿区复垦植被生物量估算系统，其特征在于，所述融合图像模块，具体包括：
光学图像单元，用于对所述光学遥感数据进行预处理，得到光学图像；
雷达图像单元，用于对所述雷达遥感数据进行预处理，得到雷达图像；
匹配单元，用于对所述光学图像和所述雷达图像进行图像配准，得到与所述光学图像对应的配准雷达图像；
融合图像单元，用于利用小波主成分算法将所述光学图像和与所述光学图像对应的配准雷达图像进行图像融合，得到融合图像。
8.根据权利要求7所述的矿区复垦植被生物量估算系统，其特征在于，所述光学图像单元，具体包括：
光学图像子单元，用于对所述光学遥感数据依次进行辐射定标和大气校正，得到光学图像；
所述雷达图像单元，具体包括：
雷达图像子单元，用于对所述雷达遥感数据依次进行辐射定标、斑点噪声抑制、地形校正和超分辨率重建，得到雷达图像；
所述匹配单元，具体包括：
地物特征点子单元，用于以所述光学图像为标准参考图像，以所述雷达图像为待配准图像，选取所述光学图像和所述雷达图像上的地物特征点；
非线性关系子单元，用于根据所述地物特征点，确定所述光学图像和所述雷达图像的非线性关系；
匹配子单元，用于利用所述非线性关系将所述雷达图像匹配到所述光学图像的坐标系下，得到与所述光学图像对应的配准雷达图像。
9.根据权利要求8所述的矿区复垦植被生物量估算系统，其特征在于，所述融合图像单元，具体包括：
主成分子单元，用于对所述光学图像进行主成分正变换，得到所述光学图像的主成分；所述主成分包括：第一主成分和剩余主成分；
直方图匹配子单元，用于将所述光学图像的第一主成分和所述配准雷达图像进行直方图匹配，得到匹配第一主成分和匹配雷达图像；
小波变换子单元，用于对所述匹配第一主成分和所述匹配雷达图像分别进行小波变换，得到多个主成分小波分量和多个雷达小波分量；
小波逆变换子单元，用于将多个所述主成分小波分量和多个所述雷达小波分量进行重新组合，得到多个小波分量，并对多个所述小波分量进行小波逆变换，得到第一主成分图像；
融合图像子单元，用于对所述第一主成分图像和所述剩余主成分进行逆变换，得到融合图像。
10.根据权利要求9所述的矿区复垦植被生物量估算系统，其特征在于，所述建立模块，具体包括：
遥感特征变量因子单元，用于对所述遥感特征变量和所述生物量数据进行相关性分析，得到相关性值大于预设相关性值的遥感特征变量因子；
生物量估算模型单元，用于将所述遥感特征变量因子作为自变量，所述生物量数据作为因变量，利用回归模型建模，得到建立的生物量估算模型。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010151513.8A CN111337434A (zh) | 2020-03-06 | 2020-03-06 | 一种矿区复垦植被生物量估算方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010151513.8A CN111337434A (zh) | 2020-03-06 | 2020-03-06 | 一种矿区复垦植被生物量估算方法及系统 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN111337434A true CN111337434A (zh) | 2020-06-26 |
Family
ID=71182166
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010151513.8A Pending CN111337434A (zh) | 2020-03-06 | 2020-03-06 | 一种矿区复垦植被生物量估算方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111337434A (zh) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111966952A (zh) * | 2020-07-21 | 2020-11-20 | 电子科技大学 | 一种叶面积指数数据过滤方法 |
| CN112580484A (zh) * | 2020-12-14 | 2021-03-30 | 中国农业大学 | 基于深度学习的遥感影像玉米秸秆覆盖识别方法及装置 |
| CN112711989A (zh) * | 2020-12-15 | 2021-04-27 | 中国农业大学 | 基于雷达遥感与光学遥感的玉米秸秆覆盖度估算方法 |
| CN113269429A (zh) * | 2021-05-19 | 2021-08-17 | 青岛星科瑞升信息科技有限公司 | 一种基于水生态效益的生态环境质量评价方法 |
| CN113468361A (zh) * | 2021-06-30 | 2021-10-01 | 中国科学院地理科学与资源研究所 | Ndvi时序数据补偿重建方法、装置和电子设备 |
| CN113505635A (zh) * | 2021-05-24 | 2021-10-15 | 中国农业大学 | 基于光学和雷达的冬小麦与大蒜混种区识别方法及装置 |
| CN113962248A (zh) * | 2021-12-01 | 2022-01-21 | 中国农业大学 | 一种确定草原地上生物量的方法、装置以及存储介质 |
| CN114581348A (zh) * | 2022-02-16 | 2022-06-03 | 三峡大学 | 一种基于植物群落行为的图像融合方法 |
| CN115993336A (zh) * | 2023-03-23 | 2023-04-21 | 山东省水利科学研究院 | 一种输水渠渠道两侧植被损坏监测方法及预警方法 |
| CN117368118A (zh) * | 2023-10-09 | 2024-01-09 | 太原理工大学 | 一种基于多光谱和点云数据处理的矿区生物量监测方法 |
| RU2822145C2 (ru) * | 2022-07-22 | 2024-07-02 | Публичное Акционерное Общество "Сбербанк России" (Пао Сбербанк) | Устройство и способ определения урожайности сельскохозяйственных культур |
-
2020
- 2020-03-06 CN CN202010151513.8A patent/CN111337434A/zh active Pending
Non-Patent Citations (7)
| Title |
|---|
| NISHA BAO ET.AL: "Biomass Estimation for Semiarid Vegetation and Mine Rehabilitation Using Worldview-3 and Sentinel-1 SAR Imagery", 《REMOTE SENSING》 * |
| 匡纲要等: "《极化合成孔径雷达基础理论及其应用》", 30 June 2011, 国防科技大学出版社 * |
| 向海燕等: "基于Sentinel-1A极化SAR数据与面向对象方法的山区地表覆被分类", 《自然资源学报》 * |
| 焦明连等著: "《合成孔径雷达干涉测量理论与应用》", 31 March 2009, 测绘出版社 * |
| 王一丁等: "《数字图像处理》", 31 August 2015, 西安电子科技大学 * |
| 王野乔等: "《长白山地理系统与生态安全 第4辑》", 30 April 2015, 东北师范大学出版社 * |
| 龙四春著: "《DlnSAR改进技术及其在沉降监测中的应用》", 31 July 2012, 测绘出版社 * |
Cited By (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111966952A (zh) * | 2020-07-21 | 2020-11-20 | 电子科技大学 | 一种叶面积指数数据过滤方法 |
| CN111966952B (zh) * | 2020-07-21 | 2023-04-18 | 电子科技大学 | 一种叶面积指数数据过滤方法 |
| CN112580484A (zh) * | 2020-12-14 | 2021-03-30 | 中国农业大学 | 基于深度学习的遥感影像玉米秸秆覆盖识别方法及装置 |
| CN112580484B (zh) * | 2020-12-14 | 2024-03-29 | 中国农业大学 | 基于深度学习的遥感影像玉米秸秆覆盖识别方法及装置 |
| CN112711989B (zh) * | 2020-12-15 | 2024-03-05 | 中国农业大学 | 基于雷达遥感与光学遥感的玉米秸秆覆盖度估算方法 |
| CN112711989A (zh) * | 2020-12-15 | 2021-04-27 | 中国农业大学 | 基于雷达遥感与光学遥感的玉米秸秆覆盖度估算方法 |
| CN113269429A (zh) * | 2021-05-19 | 2021-08-17 | 青岛星科瑞升信息科技有限公司 | 一种基于水生态效益的生态环境质量评价方法 |
| CN113269429B (zh) * | 2021-05-19 | 2022-03-18 | 青岛星科瑞升信息科技有限公司 | 一种基于水生态效益的生态环境质量评价方法 |
| CN113505635A (zh) * | 2021-05-24 | 2021-10-15 | 中国农业大学 | 基于光学和雷达的冬小麦与大蒜混种区识别方法及装置 |
| CN113505635B (zh) * | 2021-05-24 | 2024-05-31 | 中国农业大学 | 基于光学和雷达的冬小麦与大蒜混种区识别方法及装置 |
| CN113468361A (zh) * | 2021-06-30 | 2021-10-01 | 中国科学院地理科学与资源研究所 | Ndvi时序数据补偿重建方法、装置和电子设备 |
| CN113962248B (zh) * | 2021-12-01 | 2022-03-18 | 中国农业大学 | 一种确定草原地上生物量的方法、装置以及存储介质 |
| CN113962248A (zh) * | 2021-12-01 | 2022-01-21 | 中国农业大学 | 一种确定草原地上生物量的方法、装置以及存储介质 |
| CN114581348A (zh) * | 2022-02-16 | 2022-06-03 | 三峡大学 | 一种基于植物群落行为的图像融合方法 |
| CN114581348B (zh) * | 2022-02-16 | 2024-04-30 | 三峡大学 | 一种基于植物群落行为的图像融合方法 |
| RU2822145C2 (ru) * | 2022-07-22 | 2024-07-02 | Публичное Акционерное Общество "Сбербанк России" (Пао Сбербанк) | Устройство и способ определения урожайности сельскохозяйственных культур |
| CN115993336A (zh) * | 2023-03-23 | 2023-04-21 | 山东省水利科学研究院 | 一种输水渠渠道两侧植被损坏监测方法及预警方法 |
| CN117368118A (zh) * | 2023-10-09 | 2024-01-09 | 太原理工大学 | 一种基于多光谱和点云数据处理的矿区生物量监测方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111337434A (zh) | 一种矿区复垦植被生物量估算方法及系统 | |
| Li et al. | An evaluation of the use of atmospheric and BRDF correction to standardize Landsat data | |
| WO2023087630A1 (zh) | 一种利用遥感构建指数估算秸秆残留农田土壤盐分的方法 | |
| Pacheco et al. | Evaluating multispectral remote sensing and spectral unmixing analysis for crop residue mapping | |
| Rainey et al. | Mapping intertidal estuarine sediment grain size distributions through airborne remote sensing | |
| Roy et al. | Conterminous United States demonstration and characterization of MODIS-based Landsat ETM+ atmospheric correction | |
| Jin et al. | Evaluation of topographic effects on multiscale leaf area index estimation using remotely sensed observations from multiple sensors | |
| CN110136194A (zh) | 基于星载多光谱遥感数据的积雪覆盖度测算方法 | |
| CN107229913A (zh) | 基于高分卫星遥感数据结合建筑高度的人口密度分析系统 | |
| Persson et al. | Combining TanDEM-X and Sentinel-2 for large-area species-wise prediction of forest biomass and volume | |
| CN114022783A (zh) | 基于卫星图像的水土保持生态功能遥感监测方法和装置 | |
| CN114819737B (zh) | 公路路域植被的碳储量估算方法、系统及存储介质 | |
| Zhang et al. | Reconstruction of GF-1 soil moisture observation based on satellite and in situ sensor collaboration under full cloud contamination | |
| Ullmann et al. | Data Processing, feature extraction, and time-series analysis of Sentinel-1 Synthetic Aperture Radar (SAR) imagery: examples from Damghan and Bajestan Playa (Iran) | |
| Zhou et al. | Assessment of bidirectional reflectance effects on desert and forest for radiometric cross-calibration of satellite sensors | |
| CN117826112B (zh) | 一种基于sar的土壤含水量的反演方法 | |
| Oyoshi | Hourly LST monitoring with Japanese geostationary satellite MTSAT-1R over the Asia-Pacific region | |
| Neuhauser et al. | Multi-scale statistical properties of disaggregated SMOS soil moisture products in Australia | |
| Gill et al. | Estimates of bare ground and vegetation cover from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) short-wave-infrared reflectance imagery | |
| Brivio et al. | Urban pattern characterization through geostatistical analysis of satellite images | |
| CN115830476A (zh) | 地形因子空间降尺度方法 | |
| Song | Cross-sensor calibration between Ikonos and Landsat ETM+ for spectral mixture analysis | |
| Hong et al. | Estimating within-field variations in soil properties from airborne hyperspectral images | |
| Vancutsem et al. | An assessment of three candidate compositing methods for global MERIS time series | |
| Al-Ali et al. | Potionential of Spectral Indices for Halophyte Vegetation Cover Detection in Arid and Salt-Affected Landscape |
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 | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20200626 |