如图 1所示，联合试验地点为中国科学院怀来遥感综合试验站，位于河北省怀来县(东经115°48′18"E, 北纬40°21′22"N)，该试验站所属区域具有华北平原和华北平原向蒙古高原过渡的双重生态地理特征，同时站点所在农牧区域的典型植被种类丰富，是定量遥感正、反演模型研发和验证基地。

图  1  试验区域示意图

试验时间为2014年7月3日至2014年9月30日，每天进行主被动同步观测试验。根据试验区域实际生长农作物选取玉米作为观测对象，为了配合主被动观测试验，在试验区域同步开展了地面测量，获取了大量实测数据，包括植被参数、地表参数以及环境参数等，其中利用烘干法获取农作物覆盖下土壤样本的土壤含水量作为真实值。

散射计系统包括4个波段：L波段(中心频率2.0 GHz)、S波段(中心频率3.1 GHz)、C波段(中心频率5.3 GHz)、X波段(中心频率10 GHz)，雷达体制为连续调频波(FM-CW)，由电子科技大学进行安装调试，双天线系统保持单基地雷达的性能测量后向散射系数的范围是-45~20 dB，测量目标的距离范围是10~100 m。散射测量系统配置4种不同的双抛物面天线，每个天线具有目标后向散射传输信道的线性极化(H或V)和接收信道检测的线性极化(H或V)。系统安装在塔吊上，可以进行方位角0°~360°，俯仰角0°~90°测量^[7-8]。主要参数如表 1所示。

辐射计系统由4部分组成：4频段8通道微波辐射计、液压升降平台、载运车辆和供电系统。辐射计4套天线并排置于云台顶端的支架上，云台底端固定于液压升降平台顶部，升降平台底座固定在卡车车厢底部。平台升降和天线旋转、俯仰以及数据采集由控制台和计算机控制^[9]。主要参数如表 2所示。

如图 2所示，散射计与辐射计所在的高度分别为15.6 m和6 m符合远场测量条件，并且在观测方向保持直线。作物在观测区域内生长条件相同可以认为是大面积均匀地物，塔吊轨道与农作物种植的行列方向保持平行。每次测量完毕，车载平台与塔吊平台沿着轨道方向(轨道长$L$ m)移动固定间距(2 m)到下一个测量点进行独立测量。为了保证测量准确度，每次测量选择5个点进行时只让一种传感器进行测量，其他设备保持关闭。入平台允许散射计和辐射计同时列距每隔一米进行一次重复测量。

图  2  辐射计与散射计空间位置

以散射计为例，本文中使用的样本数为：

式中，i是入射角度${{\mathit{\boldsymbol{\theta }}}_e}$的数量，不考虑接近垂直入射情况，$i = 10, 11, 12, \cdots , 70$；${n_p}$是极化数，本文是全极化测量，${n_p}$为常数4；轨道长度$L$实测为10 m。

融合有效性引入分级贝叶斯验证样本数为：

式中，$R_m^2$是期望的相关系数；$p$为验证的样本数，且验证的交叉相关系数不低于${\rm{pp}}$%。假设相关系数为0.5，验证交叉系数不低于80%，那么有效验证样本数为22个，因此测量的样本数量满足条件。

本文使用了两种基于贝叶斯理论的融合方法，融合过程如图 3所示。但是二者有着不同的优势，针对本文建立的土壤含水量融合，分级贝叶斯算法具备层次建模的优点，以条件概率建立各层之间的联系，把同化独立性放宽至条件独立性，更符合辐射计与散射计实际独立测量情况。在充分考虑了同化土壤含水量中的不确定性以及不确定性对同化结果的影响下，消除了线性和高斯假设的约束，但是相比于传统贝叶斯融合算法，其对于陆面数据融合发展程度较为落后，对于先验知识应用的合理性有待提高。

图  3  融合过程

1) 水云模型。水云模型由于输入参数较少，描述植被覆盖区的微波散射机制时较为实用，因而得到了非常广泛的应用。水云模型中描述的后向散射包含两部分：

① 农作物植被层直接后向散射所得的体散射部分；②经过植被层两次衰减所得的地表后向散射部分。水云模型为^[10]：

式中，${\mathit{\boldsymbol{\sigma}}} _{{\rm{can}}}^0({\mathit{\boldsymbol{\theta }}})$为植被覆盖区总的微波后向散射系数；${\mathit{\boldsymbol{\sigma}}} _{{\rm{veg}}}^0({\mathit{\boldsymbol{\theta }}})$为植被层的后向散射系数；${{\mathit{\boldsymbol{\gamma }}}^2}({\mathit{\boldsymbol{\theta }}})$为植被层的双层衰减因子；${\mathit{\boldsymbol{\sigma}}} _{{\rm{soil}}}^0({\mathit{\boldsymbol{\theta }}})$为土壤表面的直接后向散射系数；A和B参数跟植被类型相关；${{\mathit{\boldsymbol{m}}}_{{\rm{veg}}}}$为植被含水量。

2) ω-τ模型。当土壤表面覆盖有植被时，可以使用ω-τ模型来描述植被影响下的土壤含水量，土壤表面微波辐射模型具有周期性。光学厚度为：

式中，${\rm{VWC}}$是植被冠层的含水量；${{\mathit{\boldsymbol{\tau }}}_{{\rm{NAD}}}}$为星下点光学厚度；$b$和${C_{{\rm{pol}}}}$是经验系数，取值分别为0.2和0.4^[11]。单次散射反照率$\omega $取值为0.05。植被层不但衰减和散射下垫面土壤的辐射能量，而且自身也产生热辐射能量。

假设水云模型输出的土壤含水量序列为${{\mathit{\boldsymbol{Y}}}_t}$，${{\mathit{\boldsymbol{Y}}}_t} = {[{\mathit{\boldsymbol{Y}}}({V_1}, t){\mathit{\boldsymbol{Y}}}({V_1}, t) \cdots {\mathit{\boldsymbol{Y}}}({V_N}, t)]^{\rm{T}}}$，其中${\rm{T}}$为转置。ω-τ模型输出的土壤含水量序列为${{\mathit{\boldsymbol{X}}}_t}$，${{\mathit{\boldsymbol{X}}}_t} = {[{\mathit{\boldsymbol{X}}}({V_1}, t){\mathit{\boldsymbol{X}}}({V_2}, t) \cdots {\mathit{\boldsymbol{X}}}({V_N}, t)]^{\rm{T}}}$，其中${V_1}, $${V_2}, \cdots , {V_N}$为N个测量样本点。水云模型输出土壤含水量按极化分为4层：HH、VV、HV、VH，且每个测量时间点均同时有${{\mathit{\boldsymbol{X}}}_t}$输出用于土壤含水量融合。数据模型中${{\mathit{\boldsymbol{X}}}_t}$由真实概率分布和误差两部分组成，本文将真实过程定义为${\mathit{\boldsymbol{O}}}({{\mathit{\boldsymbol{S}}}_i}, t)$，${{\mathit{\boldsymbol{S}}}_i}$为测量点，数据以${{\mathit{\boldsymbol{X}}}_t}$作为真实过程模型主变量，以${{\mathit{\boldsymbol{Y}}}_t}$作为真实过程模型的协变量分析，同时误差模型包含测量误差${\varepsilon _t}$和全极化测量样本误差${\varepsilon _p}$，具体表达式为：

分级贝叶斯网络算法融合的过程模型中，散射计和辐射计独立测量且保持足够的样本数，因此考虑两种情况，将${{\mathit{\boldsymbol{X}}}_t}$或${{\mathit{\boldsymbol{Y}}}_t}$分别作为贝叶斯后验推理，选择其中一个作为主数据模型。

根据分级贝叶斯的已有过程模型，可以定义估计真实值的过程模型为^[12]：

式中，${\xi _{st}}$是一个常数项表征时空中土壤含水量的平均值；${\xi _p}$是表征极化分层中土壤含水量的平均值；$\rho {{\mathit{\boldsymbol{O}}}_{t - 1}}$是迭代过程中上一个时刻土壤含水量对当前土壤含水量的贡献大小，其中$0 < \rho < 1$，同时也保证了在时间序列上保持连续性；${{\mathit{\boldsymbol{\beta }}}_0}$和${{\mathit{\boldsymbol{\beta }}}_S}$作为随时间空间变化的融合因子将水云模型输出土壤含水量序列${{\mathit{\boldsymbol{Y}}}_t}$进行修正；${\mathit{\boldsymbol{\eta }}}$为过程模型随时空变化的误差项。令${{\mathit{\boldsymbol{V}}}_t}{\rm{ = }}{\xi _{st}} + {\xi _p} + \rho {{\mathit{\boldsymbol{O}}}_t}{\rm{ + }}({{\mathit{\boldsymbol{\beta }}}_0} + {{\mathit{\boldsymbol{\beta }}}_S}){{\mathit{\boldsymbol{Y}}}_t}$，表示所有的层状网络参数，已知观测数据以及待推理变量，对推理后验概率${\mathit{\boldsymbol{p}}}({\mathit{\boldsymbol{\theta}}} , w\left| z \right.)$取对数，有：

基于数据模型和过程模型产生的参数主要分为误差类参数、方差类参数、固有类参数等，它们均为随机变量，理论参数的后验概率分布推导都基于贝叶斯推理。但是在实际处理中，定义误差类参数为t分布，方差类参数为正态分布，固有类参数为均值分布。

利用经典贝叶斯融合方法融合散射计和辐射计估计土壤含水量的思路是首先建立理论值与测量值之间的关系，再引入乘性误差${{\mathit{\boldsymbol{R}}}_1}$和${{\mathit{\boldsymbol{R}}}_2}$，乘性误差${{\mathit{\boldsymbol{R}}}_1}$和${{\mathit{\boldsymbol{R}}}_2}$均和土壤含水量独立且不相关，以水平极化为例，给出数据模型为^[13]：

式中，${\mathit{\boldsymbol{\sigma}}} _{{\rm{HH}}}^0$为水平极化后向散射系数；${{\mathit{\boldsymbol{T}}}_b}_{{{\rm{H}}_m}}$为水平极化亮温。根据贝叶斯理论，条件概率密度函数可以写成：

式中，${{\mathit{\boldsymbol{M}}}_v}$是土壤含水量；${{\mathit{\boldsymbol{P}}}_{{\rm{prior}}}}({{\mathit{\boldsymbol{M}}}_v})$是土壤含水量先验信息，由归一化概率密度函数表示；分布函数${{\mathit{\boldsymbol{P}}}_{{\rm{post}}}}({{\mathit{\boldsymbol{T}}}_b}_{{\rm{H}}, }{\mathit{\boldsymbol{\sigma}}} _{{\rm{HH}}}^0|{{\mathit{\boldsymbol{M}}}_v})$是似然函数，描述了基于测量值和噪声计算获得的概率分布函数；${\mathit{\boldsymbol{P}}}({{\mathit{\boldsymbol{T}}}_b}_{{\rm{H}}, }{\mathit{\boldsymbol{\sigma}}} _{{\rm{HH}}}^0)$是一个标准化因子。

似然函数${{\mathit{\boldsymbol{P}}}_{{\rm{post}}}}({{\mathit{\boldsymbol{T}}}_b}_{{\rm{H}}, }{\mathit{\boldsymbol{\sigma}}} _{{\rm{HH}}}^0|{{\mathit{\boldsymbol{M}}}_v})$可以由概率密度函数${\mathit{\boldsymbol{P}}}({{\mathit{\boldsymbol{R}}}_1}, {{\mathit{\boldsymbol{R}}}_2})$表示，有：

随机变量${{\mathit{\boldsymbol{R}}}_1}$和${{\mathit{\boldsymbol{R}}}_2}$已经在前面给出，分别代表散射计和辐射计测量中引入的误差，因此该似然函数可以转化为联合概率密度函数。一般的测量场景中假设${{\mathit{\boldsymbol{R}}}_1}$和${{\mathit{\boldsymbol{R}}}_2}$服从高斯分布：

式中，${{\mathit{\boldsymbol{\mu }}}_1}$、${{\mathit{\boldsymbol{\mu }}}_2}$和标准偏差${{\mathit{\boldsymbol{\sigma}}} _1}$、${{\mathit{\boldsymbol{\sigma}}} _2}$均由${{\mathit{\boldsymbol{R}}}_1}$和${{\mathit{\boldsymbol{R}}}_2}$的联合概率密度决定，根据最大似然函数定义，${{\mathit{\boldsymbol{\mu }}}_1}$、${{\mathit{\boldsymbol{\mu }}}_2}$、${{\mathit{\boldsymbol{\sigma}}} _1}$和${{\mathit{\boldsymbol{\sigma}}} _2}$是联合分布函数的最大值。给出了特定的后向散射系数和亮温，也就确定了联合密度函数。

由于测量试验周期较长，随着测量时间的不同，农作物生长状态以及环境参数会发生变化，同时下垫面土壤含水量也会发生变化。如图 4所示，将融合前水云模型和ω-τ模型，以及分级贝叶斯和经典贝叶斯融合后的土壤含水量输出结果按照测量时间顺序比较，图中的土壤含水量真实值是利用烘干法将实测土样在温度85 ℃烘干10 h的结果。

图  4  观测周期实测土壤含水量与模型结果比较

从秧苗期至拔节期，水云模型反演结果相对于真实土壤含水量偏大，ω-τ模型虽然变化趋势与真实值接近，但是提取精度偏差较大。原因是生长前期植被矮小使地表土壤裸露部分较多，同时农作物前期生长需要翻土造成土壤中杂质进入地表。经过分级贝叶斯与经典贝叶斯融合后，提取明显更加接近真实土壤含水量，其中经典贝叶斯相比分级贝叶斯网络更加准确。从拔节期开始，土壤含水量变化后分级贝叶斯估计含水量的精度与经典贝叶斯一致。

从拔节期至抽雄期，土壤含水量变化频率开始增加，冠层厚度也在迅速增加，提取土壤含水量变得复杂。融合前ω-τ模型提取土壤含水量精度明显提高，变化趋势也同真实值保持较高相关性，水云模型提取结果存在一些波动，但是精度仍然好于ω-τ模型。通过对比两种贝叶斯融合估计土壤含水量表明，分级贝叶斯比经典贝叶斯算法更加准确，因为参数化层次模型能够更好地估计变化，而经典贝叶斯对含水量变化较大时估计有明显滞后。

从抽雄至成熟期，为了检验两种融合方法的估计精度，人为增加了灌溉次数，因此土壤含水量增长很快，而且由于该地区属于沙土，土壤持水性较差，导致含水量变化快。所以使用单一传感器在这种情况下土壤含水量提取精度明显下降，水云模型输出结果相对于真实土壤含水量偏低，ω-τ模型虽然变化趋势仍然保持较好，但是整个后期模型输出结果与真实土壤含水量相比都偏大。利用贝叶斯融合算法，两种贝叶斯融合针对后期复杂性土壤含水量的提取均表现出明显优势。两种融合相比，分级贝叶斯能够在精度以及变化响应方面更好。因为分级贝叶斯利用条件概率建立数据、过程、参数之间联系的特点，当土壤含水量发生变化时，能更好地在各层中描述状态量的变化并转化为已有数据条件下的推理过程和后验分布概率。

预测误差是模型值与实际值的差距，有几种常用的衡量预测误差的指标，但任何单一的一种指标都很难全面地评价一个预测模型。在实际中可将它们结合使用，根据实际的要求，如精度较高或偏离较低的模型，选择较为合适的模型。

针对模型输出结果和实测土壤含水量，为了更好地分析模型输出之间的误差大小，表 3给出了平均平方误差(MSE)、平均绝对差值(MAD)、平均相对误差(MRE)的统计结果。MSE是广泛使用的误差统计指标，能较好反应估计的精度，缺点是无法衡量无偏性。MAD容易求得，要求计算简单时可使用，且避免了误差相消的问题，可准确反映实际预测误差的大小。MRE能够较好的反应相对偏差大小及模型输出时间序列的土壤含水量偏离大小。

从拔节期至抽雄期，植被层稀疏且土壤含水量变化不大，由于辐射计视场较大受到噪声干扰较多，ω-τ模型反演精度最低，与真实值最大误差达到12.2%，其平均误差和相关系数分别为8.1%和0.78。水云模型反演土壤含水量相对准确最大误差为6.7%, 其平均误差和相关系数分别为3.3%和0.63，经过贝叶斯融合算法后，土壤含水量精度明显提高且相关系数均大于融合前。其中，分级贝叶斯融合后与真实值的最大误差为2.3%，与真实值比较接近；而经典贝叶斯算法融合结果比真实值偏小，最大误差为3.4%。从分蘖后期到扬花期，植被层开始增加且含水量开始逐渐增加，ω-τ模型反演土壤含水量精度有所提高但还是与真实值有差距，平均误差为9.8%，但是相关系数较高为0.73，保持了与真实土壤含水量相同的变化趋势。此时期分级贝叶斯与经典贝叶斯相比含水量平均误差减少了3.2%，相关系数也提高到81%。从抽雄至成熟期，植被生物量达到最大，此时土壤含水量随温度进行人工灌溉，因此土壤含水量有较大变化。可以看到，在土壤含水量发生明显变化的时候，水云模型反演土壤含水量精度下降，最大误差达到了11.2%，ω-τ模型反演土壤含水量平均误差和相关系数分别为5.3和0.76。这时期经典贝叶斯融合后土壤含水量明显低于真实值，平均误差较前一时期增加2.1%，与之相比分级贝叶斯融合结果较好，误差精度小于1.7%，相关系数达到91%。

在农作物完整生长期内，利用贝叶斯融合算法与真实值相比，平均平方误差(MSE)小于3.56、平均绝对差值(MAD)小于1.36、平均相对误差(MRE)小于13.92%，同时分级贝叶斯与经典贝叶斯同真实土壤含水量的决定系数为0.77和0.60。

本文提出了一种基于贝叶斯与主被动雷达传感器的真实土壤含水量估计方法，研究对象为主要生长季的农作物。总体来看，利用主被动传感器数据反演土壤含水量各有优势，被动辐射计测量对环境参数变化敏感，能够及时反应土壤含水量发生的变化。基于贝叶斯融合的土壤含水量平均平方误差(MSE)小于3.56、平均绝对差值(MAD)小于1.36、平均相对误差(MRE)小于13.92%，决定系数大于0.60，融合二者既能保证精度的情况下，也能增强对土壤含水量变化的响应。但是该方法仍然有一些不足，首先，模型本身的精度以及模拟中模型输入参数同实际情况的偏差可能是结果误差的主要来源；其次，虽然玉米生长后期生长冠层散射较强，但生长前期地面植被对反演结果造成的影响依然不可忽视。