一种基于GF-1WFV数据地表反照率反演方法

一种基于gf-1 wfv数据地表反照率反演方法
技术领域
1.本发明涉及空间信息遥感的技术领域，尤其涉及一种基于gf-1 wfv数据地表反照率反演方法。


背景技术：

2.地表反照率是地表辐射能量平衡和地球与大气相互作用的驱动因素之一。地表反照率是一个重要参数，广泛用于地表能量平衡、中长期天气预报和全球变化研究。目前，已有基于卫星遥感数据生成的反照率数据集，如：glass(global land surface satellite)、modis(moderate-resolution imaging spectroradiometer)等，但这些反照率产品空间分辨率较低，难以准确描述景观破碎度较大区域的地表反照率的时空变化，也难以观察到详细的地表细节特征信息。国内高分辨卫星地表反照率反演几乎都是针对landsat和hj-1影像建立的，但对国产高分系列卫星地表反照率反演并不多。
3.利用遥感技术反演地表反照率包含两种方法：一种是基于统计模型获得地表反照率，从反照率定义出发，运用modtran 4.0模型计算光通量来计算地表反照率。但这种方法没有考虑地表各向异性；另一种反演方法是考虑地表的各项异性特征，通过二向性反射分布函数(brdf)来反演地表反照率。brdf描述地表各向异性特征，是入射方向和反射方向的辐射亮度的微增量之间的比值。利用brdf模型对卫星数据进行反照率反演，这种反演方法需要多角度像元反射率信息，对于单角度的卫星来说，反演比较困难。
4.申请号为202210733214.4的发明专利公开了一种基于地表反射各向异性特征分类的反照率反演方法，在各向异性平整指数afx的基础上，提出一个与afx垂直的指数，命名为垂直各向异性平整指数pafx，联合两个指数对全球地表反射各向异性特性信息进行分类，形表反射各向异性特征先验知识进行反照率的反演与验证。该发明能够更精确的对地表反射各向异性特征进行分类，提高在观测数据不充足情况下的反照率反演精度，并且简单易操作，应用较便捷。但是，上述发明没有解决gf-1高分辨卫星反演地表反照率缺少多角度反射率数据，以此无法获得gf-1地表各向异性特征，难以进行地表反照率反演。


技术实现要素：

5.针对现有方法的地表反照率空间分辨率低和反演过程中很难在一段时间获得多角度反射率数据的技术问题，本发明提出一种基于gf-1 wfv数据地表反照率反演方法，基于gf-1的多光谱wfv数据和modis反照率产品(mcd43a1)两种影像的优势，gf-1号卫星拥有16m高空间分辨率，相对于目前大多数低分辨率地表反照率能够提供更为清晰的地表特征，而modis反照率产品可以获得地表各项异性特征(brdf)。
6.为了达到上述目的，本发明的技术方案是这样实现的：一种基于gf-1 wfv数据地表反照率反演方法，其步骤如下：
7.步骤1：从陆地观测卫星数据服务平台下载所需要的gf-1 wfv影像，然后下载研究区域的mcd43a1反照率产品；
8.步骤2：对gf-1 wfv影像进行预处理，并对预处理的数据进行裁剪，得到研究区域gf-1地表反射率数据集，提取mcd43a1反照率产品的地表异性先验知识；
9.步骤3：利用gf-1地表反射率数据集对mcd43a1反照率产品的地表异性先验知识进行前向计算，对地表异性先验知识进行调整得到gf-1 wfv数据的参数brdf，对参数brdf积分得到黑白空反照率，然后加权求和得到gf-1窄波段反照率；
10.步骤4：通过对宽波段反照率进行分析，筛选出地物光谱曲线，对步骤三得到的gf-1窄波段反照率通过多元线性回归分析得到gf-1 wfv数据的可见光宽波段转换模型，得到gf-1宽波段反照率；
11.步骤5：利用mcd43a1反照率产品计算得到的宽波段反照率对得到gf-1宽波段反照率的反演结果进行验证分析。
12.优选地，所述步骤2中对gf-1 wfv影像进行预处理包括辐射定标和大气校正，使用6s模型对gf-1 wfv影像进行辐射定标和大气校正得到gf-1 wfv反射率数据，对gf-1 wfv反射率数据进行wgs-84坐标系重投影获得标准的地理坐标系后裁剪出相应的研究区域，研究区域的gf-1 wfv反射率数据组合得到研究区域gf-1地表反射率数据集。
13.优选地，所述地表异性先验知识的提取方法为：用mrt工具将mcd43a1反照率产品由正弦投影转换为wgs-84坐标系重投影，得到标准的地理坐标系；mcd43a1反照率产品的地表异性先验知识为地表参数brdf'，地表参数brdf'选择7个波段中的各个像元的红、绿、蓝三个波段的核系数f
iso
、f
vol
、f
geo
。
14.优选地，所述步骤三的实现方法为：
15.设gf-1 wfv影像的反射率为rg，把mcd43a1反照率产品的地表异性先验知识带入核驱动模型中得到modis地表反射率rm，计算反射率rg与modis地表反射率rm的比值：
16.通过调整反射率比值c
λ
对mcd43a1反照率产品的地表异性先验知识的地表参数brdf'进行调整，得到待求gf-1 wfv数据的参数brdf＝c
λ
×
brdf
′
；即得到核系数f
iso
(λ)、f
vol
(λ)、f
geo
(λ)；
17.对得到gf-1 wfv数据的参数brdf进行积分得到相对应的黑空反照率a
bsa
和白空反照率a
wsa
分别为：
[0018][0019][0020]hi
(θi)＝g
0i
+g
1i
θ2+g
2i
θ3；
[0021]
其中，θ、g
0i
、g
1i
、g
2i
分别表示像元太阳天顶角和各个核的多项式拟合表达式的系数，hi(θ))是第i核在观测半球的积分值，hi是核在入射和观测半球的积分值，分别表示f
iso
(λ)、f
vol
(λ)、f
geo
(λ)；
[0022]
将黑空反照率a
bsa
(θ)和白空反照率a
wsa
进行加权组合，得到波段λ的窄波段反照率：a(θ,λ)＝[(1-s(τ))a
bsa
(θ)]+s(τ)a
wsa
；
[0023]
其中，s(τ)为天空散射光所占比分，其大小由大气校正后的参数估算得到。
[0024]
优选地，所述核驱动模型对mcd43a1反照率产品进行归一化处理，将mcd43a1反照率产品归一化相同的数量级，将各个像元的红、绿、蓝三个波段的核系数f
iso
、f
vol
、f
geo
归一化得到归一化处理的brdf'参数：
[0025][0026]
其中，r是地表二向反射率函数，θ为太阳天顶角，为观测天顶角，为相对方位角；λ为波长；k
vol
和k
geo
分别表示为体散射核和几何光学核，f
iso
(λ)、f
vol
(λ)、f
geo
(λ)都是和波长相关的常系数。
[0027]
优选地，所述体散射核k
vol
和几何光学核k
geo
分别选择rossthick核函数和litransit核函数；
[0028]
所述modis地表反射率rm的计算方法为：将mcd43a1反照率产品的brdf参数fi
iso
(λ)、f
vol
(λ)、f
geo
(λ)带入到下面核驱动模型得到：
[0029][0030]
优选地，所述各个核的多项式拟合表达式的系数取值为：
[0031][0032]
优选地，所述步骤4的实现方法为：通过sbdart模型在美国usgs数字光谱库里筛选出190个地物光谱曲线，其中包括草地、水体、农田、城市的地物类型，输入六种大气能见度、八个太阳天顶角和三种大气模式，建立多元线性回归分析得到gf-1 wfv数据的可见光宽波段转换模型：
[0033]
α
gf-1
＝0.443αb+0.317αg+0.240αr；
[0034]
其中，αb、αg、αr分别为蓝光波段、绿光波段、红光波段的窄波段反照率，通过波段λ的窄波段反照率获得；α
gf-1
表示gf-1宽波段反照率。
[0035]
优选地，所述宽波段反照率为在一定波长范围内的地表上行与下行辐射通量之比：
[0036][0037]
式中，α(θ,λ)是宽波段反照率，λ是波段λ1到λ2的波段范围，fu(θ,λ)和fd(θ,λ)分别是上行、下行辐射通量，通过sbdart模型得到；α(θ,λ)为波段λ的窄波段反照率，θ为太阳太顶角；
[0038]
所述太阳天顶角的范围为0-80
°
，每个间隔10
°
的值；三种大气模式包括热带、中纬度冬季、中纬度夏季；
[0039]
所述步骤5进行验证分析的方法为：
[0040]
随机抽取样本点进行均方根误差分析，均方根误差为：
[0041]
其中，n为样本点的数量，rmse小于0.05；a
gf-1
表示gf-1宽波段反照率，a
modis
为modis地表反照率。
[0042]
优选地，所述modis地表反照率a
modis
的计算方法为：由mcd43a1反照率产品中红、绿、蓝三个波段各个像元的核系数f
iso
、f
vol
、f
geo
带入到modis的brdf模型中进行积分，得到这三个波段的黑空反照率以及白空反照率，然后将这三个波段的黑空反照率和白空反照率分别代入到modis地表反照率产品自身的窄波段向宽波段转换的关系式a
modis
＝0.331ar+0.424ab+0.246ag中，得到a
modis
地表反照率；其中，ar、ab、ag为mcd43a1反照率产品红光波段、蓝光波段、绿光波段的窄波段反照率。
[0043]
本发明以核驱动模型为基础，首先对gf-1号数据进行辐射定标和大气校正得到地表反射率，然后选取和研究区相对应的modis的反照率产品(mcd43a1)，进行提取地表反射各向异性先验知识，通过gf-1号卫星数据对其进行拟合，最后进行窄波段向宽波段转换得到gf-1地表反照率。
[0044]
与现有技术相比，本发明的有益效果：基于gf-1的多光谱wfv数据和modis反照率产品(mcd43a1)两种影像的优势，以核驱动模型为基础，对gf-1 wfv数据进行预处理得到处理后的影像，对gf-1号数据进行辐射定标和大气校正得到地表反射率，从modis反照率产品mcd43a1中提取地表反射各向异性先验知识，作为下垫面数据集，通过gf-1 wfv影像数据对其进行拟合，进行窄波段向宽波段转换，从而得到gf-1 wfv数据下的brdf数据集，进而得到gf-1 wfv数据高分辨地表反照率产品。通过对比modis反照率产品，反演精度能够达到精度要求。
附图说明
[0045]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0046]
图1为本发明的流程示意图。
[0047]
图2为本发明的实验对比图之一，其中，(a)为gf-1地表反照率，(b)为modis地表反照率。
[0048]
图3为本发明的样本点的拟合曲线。
[0049]
图4为本发明gf-1与modis反照率的绝对误差图。
具体实施方式
[0050]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0051]
如图1所示，一种基于gf-1 wfv数据地表反照率反演方法，其步骤如下：
[0052]
步骤1：从陆地观测卫星数据服务平台下载所需要的gf-1 wfv影像。gf-1号四台wfv相机宽幅800公里，分辨率达到16m，实现了覆盖能力广、空间分辨率高的完美结合。然后从earthdata网站下载研究区域的mcd43a1反照率产品，mcd43a1反照率产品在全球范围内应用较为广泛，是国际遥感界较为认可的产品之一。
[0053]
步骤2：对gf-1 wfv影像进行预处理，并对预处理的数据进行裁剪，得到研究区域gf-1地表反射率数据集，提取mcd43a1反照率产品的地表异性先验知识。
[0054]
对下载的gf-1 wfv卫星影像进行预处理，包含：辐射定标和大气校正，使用6s模型对gf-1 wfv影像进行辐射定标和大气校正得到gf-1 wfv反射率数据，校正后的影像消除了大气影响变得更加清晰，并进行wgs-84坐标系重投影后裁剪出相应的研究区域，wgs坐标系重投影可以获得标准的地理坐标系，然后得到研究区域gf-1地表反射率数据集。
[0055]
用mrt工具将mcd43a1反照率产品由正弦投影转换为wgs-84坐标系进行重投影，得到标准地理坐标系；mcd43a1反照率产品提供了全球16天合成的500m地表参数brdf'，地表参数brdf'包括7个波段各个像元的核系数f
iso
、f
vol
、f
geo
，其中用mrt工具选择其中红、绿、蓝三个波段的核系数f
iso
、f
vol
、f
geo
。
[0056]
步骤3：利用gf-1地表反射率数据集对mcd43a1反照率产品的地表参数brdf'进行前向计算得到modis地表反射率数据rm，然后计算相同的观测几何方向的实测反射率和模拟方向反射率之比得到c
λ
，对地表参数brdf'进行调整，得到gf-1 wfv数据的参数brdf，对其积分得到黑白空反照率，然后以一定比例相加得到gf-1窄波段反照率。
[0057]
具体实现方法为：
[0058]
步骤3.1：基于核驱动模型，从mcd43a1反照率产品中提取先验知识，由于地表各向异性特征较为复杂，因此需要对mcd43a1反照率产品进行归一化处理，对mcd43a1反照率产品归一化相同的数量级，将其中f
iso
(λ)、f
vol
(λ)、f
geo
(λ)三个部分归一化得到归一化处理的brdf'参数：
[0059][0060]
其中，r是地表二向反射率，θ为太阳天顶角，为观测天顶角，为相对方位角；λ为波长；k
vol
和k
geo
分别表示为体散射核和几何光学核，f
iso
(λ)、f
vol
(λ)、f
geo
(λ)都是和波长相关的常系数，通过下面公式(4)获得gf-1的brdf(f
iso
(λ)、f
vol
(λ)、f
geo
(λ))。为了避免出现反射率为负的情况，在选择k
vol
和k
feo
两个核上面分别选择rossthick核函数和litransit核函数，litransit核函数可以有效避免了出现拟合反射率为负数的情况，其拟合能力强。
[0061]
步骤3.2：假设gf-1 wfv反射率为rg，通过把mcd43a1反照率参数f
iso
(λ)、f
vol
(λ)、f
geo
(λ)带入核驱动模型如公式(2)中，得到modis地表反射率rm，然后得到两者反射率比值：
[0062][0063][0064]
通过调整反射率比值c
λ
对mcd43a1反照率产品的地表异性先验知识的地表参数brdf'进行调整，由mcd43a1原型参数得到待求gf-1 wfv数据的brdf
[0065]
brdf＝c
λ
×
brdf
′ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)。
[0066]
步骤3.3：对得到gf-1 wfv数据的参数brdf进行积分得到相对应的黑空反照率a
bsa
和白空反照率a
wsa
分别为：
[0067][0068][0069]hi
(θi)＝g
0i
+g
1i
θ2+g
2i
θ3ꢀꢀꢀꢀꢀꢀꢀ
(7)
[0070]
其中，θ、g
0i
、g
1i
、g
2i
分别表示像元太阳天顶角和各个核的多项式拟合表达式的系数，取值为下表1所示，核积分hi(θ)为公式(5)进行计算的一项，hi为公式(6)计算的一项与观测角度无关，hi在rossthick和litransit对应的积分值分别是0.189184和－1.137762。fi(λ)包含f
iso
(λ)、f
vol
(λ)、f
geo
(λ)三项，由公式(4)计算得到的，将两者进行相乘都可以得到λ波段的黑空反照率和白空反照率。
[0071]
表1多项式拟合表达式的系数
[0072][0073]
步骤3.4：为了能够得到窄波段反照率，将两种黑白空反照率以一定的比例进行加权组合：
[0074]
a(θ,λ)＝[(1-s(τ))a
bsa
(θ)]+s(τ)a
wsa
ꢀꢀꢀꢀꢀꢀꢀ
(8)
[0075]
其中，s(τ)为天空散射光所占比分，其大小由步骤2大气校正后的参数估算得到。a(θ)为波段λ的窄波段反照率，后续进行宽波段计算。
[0076]
步骤4：宽波段反照率的定义为在一定波长范围内的地表上行与下行辐射通量之比：
[0077][0078]
式中，α(θ,λ)是宽波段反照率，λ是波段λ1到λ2的波段范围，fu(θ,λ)和fd(θ,λ)分别是上行、下行辐射通量，通过sbdart模型得到。α(θ,λ)为波段λ的窄波段反照率，θ为太阳太顶角。sbdart模型计算地球上行、下行辐射通量，sbdart有很好的接口，使我们可高效地建立模拟辐射通量。通过sbdart模型在美国usgs数字光谱库里筛选出190典型的地物光谱曲线，其中包括草地、水体、农田、城市等地物类型，输入六种大气能见度和八个太阳天顶角(0-80
°
)每个间隔10
°
的值和三种大气模式，其中包含热带、中纬度冬季、中纬度夏季，建立多元线性回归分析得到适用于gf-1 wfv数据的可见光宽波段转换模型：
[0079]
α
gf-1
＝0.443αb+0.317αg+0.240αrꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0080]
其中，αb、αg、αr分别为蓝光波段、绿光波段、红光波段的窄波段反照率，通过公式
(7)的a(θ)获得，α
gf-1
表示gf-1宽波段反照率。
[0081]
步骤5：利用mcd43a1计算得到的宽波段反照率对由公式(10)得到gf-1地表反照率反演结果进行验证分析。
[0082]
随机抽取样本点，进行均方根误差分析，均方根误差表示观测值与真值偏差的平方与观测次数n比值的平方根：
[0083][0084]
其中，rmse精度要求为小于0.05。a
gf-1
、a
modis
分别表示gf-1地表反照率和modis地表反照率，a
gf-1
由公式(10)得到。由modis mcd43a1产品中红、绿、蓝三个波段各个像元的核系数f
iso
、f
vol
、f
geo
带入到modis brdf模型中进行积分，得到这三个波段的黑空反照率以及白空反照率，然后将这三个波段的黑空反照率和白空反照率分别代入到modis地表反照率产品自身的窄波段向宽波段转换的关系式中如公式(12)，得到a
modis
地表反照率；其中，ar、ab、ag为mcd43a1反照率产品红光波段、蓝光波段、绿光波段的窄波段反照率。
[0085]amodis
＝0.331ar+0.424ab+0.246agꢀꢀ
(12)
[0086]
在本实施例中，所述gf-1为高时空卫星，其获得了大量数据集，携带了两台全色/多光谱(p/ms)和四台wfv相机，四台wfv相机宽幅800公里，分辨率达到16m。所述modis为aqua和terra卫星上的关键传感器，每隔一到两天可以观测整个地球表面，并获得36个波段，其中7个波段被用来生成brdf产品。
[0087]
在本实施例中，结合实际，选取我国河南省郑州航空港区作为研究区，对郑州航空港区进行地表反照率反演。首先对2020年3月18号的gf-1 wfv影像进行预处理，预处理后裁剪出来郑州航空港研究区，从mcd43a1产品中提取出可见光三个波段的brdf参数，然后把这些模型参数带入到它们自身转换系数关系式中，获得modis反照率产品。将两个反照率进行交叉验证。
[0088]
图2是通过带入gf-1波段转换函数获得，通过图2可以看出，gf-1反照率和modis反照率产品在影像分布上一致。然后在试验区随机选取400个样本点，对这些样本点进行拟合分析，如图3所示，得到的rmse能够小于0.05，能够达到对反照率的精度要求。图4为gf-1卫星反照率的绝对误差均值为0.0185。
[0089]
本发明首先对gf-1 wfv数据进行预处理得到处理后的影像，然后进行提取地表各向异性的先验知识，通过拟合得到gf-1的地表brdf参数，通过对核驱动模型进行积分得到黑白空地表反照率，最后进行窄波段向宽波段转化得到gf-1地表反照率产品。通过对比modis反照率产品，反演精度能够达到精度要求。
[0090]
本发明提出了使用modis反照率产品(mcd43a1)来获得地表各向异性特征，由此得到gf-1 wfv地表各向异性信息，然后进行计算gf-1窄波段反照率，最后进行窄波段反照率向宽波段反照率转换。以郑州航空港gf-1影像实验数据为例，通过提取mcd43a1地表各向异性先验知识，来克服gf-1高分辨率卫星无法获取多角度反射率数据，进而无法获得地表各向异性特征，本发明为gf-1 wfv数据反演地表反照率提供一种解决办法，尤其对gf-1卫星反演地表反照率具有重要的应用价值。
[0091]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。