一种基于谐波幅度相位分析的粮油作物散射波谷识别方法

1.本发明属于星载雷达数据处理与遥感制图技术领域，具体涉及一种基于谐波幅度相位分析的粮油作物散射波谷识别方法。


背景技术：

2.粮油作物包括谷物、豆类和油料作物，这些作物分别含有丰富的淀粉、植物蛋白质和植物油脂，满足了人们基本的食物和能量需求。长时序遥感对地观测是大面积粮油作物种植动态调查的有效手段，光学卫星成像极易受到频繁的云雨天气影响，而搭载合成孔径雷达的卫星可以克服云雾干扰，提供连续稳定的长时序高质量雷达散射数据，因此是光学卫星的重要替代。然而，利用雷达散射时序监测粮油作物种植需要解决两个关键问题，分别是时序去噪和时序特征识别。
3.常用的时序去噪方法主要分为窗口统计法和局部拟合法，前者计算特定滑动窗口内的时序统计值(如均值、中值、最大值和最小值)以抑制异常值，后者对特定时间范围内的时序数值进行函数拟合(如高斯、逻辑斯蒂和多项式拟合)以抑制异常波动，两种方法均能实现对时序噪声的压制。窗口统计法依赖对滑动窗口尺寸的合理设置；局部拟合法依赖对局部时序趋势的合理假设。然而，散射时序存在复杂多样的随机噪声和局部趋势，特定滑动窗口和拟合函数对不同散射时序的适应性难以保证。
4.时序特征主要是指时序曲线的波峰(局部极大值)和波谷(局部极小值)，雷达散射波峰的物理意义并不明确，因为粮油作物抽穗、土壤粗糙度和含水量升高都可能导致散射波峰，而散射波谷一般代表地面裸露且平坦，对应粮油作物耕种前的平整土地时期，因此散射波谷是粮油作物种植监测中最可靠的时序特征。常用的时序特征识别方法主要是差分搜索法，通过计算时序数据的一阶差分，搜索一阶差分的零值点，结合零值点附近的曲线升高或降低趋势，判断零值点对应散射时序的波峰还是波谷。然而，差分搜索法需要对时序数据进行逐点运算，因此对运算内存和时间的消耗较大。此外，从更长的时间窗口来看，差分搜索法识别的波谷有时并不是散射低值，因此还需要对波谷进行复杂的二次判别，以确定波谷对应平坦地表。


技术实现要素：

5.针对上述存在的问题或不足，为解决现有散射波谷识别中的时序去噪适应性差和时序特征识别效率低的问题，本发明提供了一种基于谐波幅度相位分析的粮油作物散射波谷识别方法，通过谐波拟合实现散射时序重构和去噪，通过幅度相位分析实现散射波形的检测和散射波谷的快速识别。
6.一种基于谐波幅度相位分析的粮油作物散射波谷识别方法，具体步骤如下：
7.步骤1.时序谐波拟合；
8.对目标区域，获取雷达vh(vertically transmitting and horizontally receiving，垂直发射并水平接收)极化散射数据，构成散射时序s[t]，t代表时间，取值范围
归一化为0到1，对应全年的首日到尾日；
[0009][0010]
上式中a为常数项，代表散射年度均值。i依次取值1，2，3，代表余弦项的阶次。ai和分别依次为第i阶余弦项的幅度和相位，a、ai、的值未知，通过最小二乘拟合获得。
[0011]
步骤2.谐波指标计算；
[0012]
第i阶余弦项拥有i个波谷，波谷日期由计算得到。第1阶余弦项拥有1个波谷，设波谷日期为v1。第2阶余弦项拥有2个波谷，设波谷日期分别为v
2a
和v
2b
。第3阶余弦项拥有3个波谷，设波谷日期分别为v
3a
，v
3b
和v
3c
。
[0013][0014]
计算第i阶余弦项幅度占全部余弦项幅度之和的比例：
[0015]
pi＝ai/(a1+a2+a3)
[0016]
计算第i阶余弦项相对于第j阶余弦项的相位移动：
[0017][0018]
步骤3.散射波形检测；
[0019]
将散射波形分为三种类型：单波谷型、双波谷型和三波谷型。
[0020]
当pi和满足以下四种情况时，判定散射波形为单波谷型s型，具体包括：
[0021]
s0型：p1》2/3
[0022]
s1型：
[0023]
s2型：
[0024]
s3型：
[0025]
当pi和满足以下四种情况时，判定散射波形为双波谷型d型，具体包括：
[0026]
d0型：p1《1/3，p2≥1/3或p3≤1/2
[0027]
d1型：
[0028]
d2型：
[0029]
d3型：或
[0030]
当pi和满足以下情况时，判定散射波形为三波谷型t型，即：
[0031]
t0型：p1《1/3，p2《1/3，p3》1/2；
[0032]
步骤4.粮油作物散射波谷识别：包括基准分量确定、最近波谷查找、加权平均计算和波谷量级限制四步完成。
[0033]
基准分量确定：对于单波谷型，基准分量为第1阶余弦项；对于双波谷型，基准分量为第2阶余弦项；对于三波谷型，基准分量为第3阶余弦项。基准分量的波谷数量与散射波形的波谷数量相同，且散射波形中任一波谷的日期v
p
，都可由与之距离最近的基准分量波谷的日期v
x
来近似表示。
[0034]
最近波谷查找：根据步骤2中波谷日期的计算方法，求解得到三个余弦项中所有波谷的日期。对于v
p
相应的v
x
，分别在其余两种余弦项中查找与v
x
相应的基准分量波谷距离最近波谷的日期vy和vz。至此，成功查找到与散射波形中目标波谷日期v
p
距离最近的三种余弦项中的波谷相应日期v
x
，vy和vz。
[0035]
加权平均计算：设v
x
、vy和vz对应余弦项的阶次i分别为i
x
、iy和iz，幅度ai分别为a
x
、ay和az，则可以由v
x
，vy和vz加权平均计算得到v
p
1，v
p
1为v
p
的逼近值，公式如下：
[0036][0037]
波谷量级限制：将散射波谷日期v
p
1作为t值代入步骤1中的谐波拟合公式以获得散射波谷量级，将0.02设置为波谷量级的宽松阈值，只保留量级低于0.02的散射波谷，这些波谷即为粮油作物散射波谷。至此，完成粮油作物散射波谷识别。
[0038]
步骤1到4涉及的原理：
[0039]
时序谐波拟合原理：粮油作物的年度散射时序可以看作多个周期性余弦过程的叠加，农田中的粮油作物一年最多可以种植三轮，因此三种阶次余弦项可以在抑制更高阶随机噪声的同时，充分模拟粮油作物的多样化散射时序，它们的幅度表示周期性的耕地变化程度强弱，它们的相位表示周期性的耕地变化时间。
[0040]
谐波指标计算原理：初始状态下，三种阶次余弦项拥有相同的幅度和相位。进行最小二乘拟合后，三种阶次余弦项的幅度和相位发生改变，通过计算各余弦项的幅度占比可以确定各余弦项对散射时序的能量贡献，通过计算余弦项两两之间的相位移动可以确定余弦项波谷之间的位置关系。
[0041]
散射波形检测原理：散射波形是三种阶次余弦项共同作用的结果，余弦项的幅度占比越大，对散射波形的影响越大。此外，当余弦项之间的相位移动处于不同的区间时，余弦项波谷之间发生增强和削弱作用的时刻会发生变化。通过对余弦项的幅度占比和相位移动进行数值分段，可以实现散射波形的类型检测。
[0042]
粮油作物散射波谷识别：散射波形检测后可以获知散射时序中的波谷数量，每个散射波谷都由与之相邻的三种阶次余弦项的波谷共同贡献。三种阶次余弦项的波谷日期可以通过相位计算获得，散射波谷日期可以通过三种阶次余弦项的波谷日期的加权平均获得。波谷日期确定后可根据谐波拟合公式计算出散射波谷量级。粮油作物种植前平整土地导致散射强度非常低(0.01左右)，而其他波谷量级通常在0.02以上，因此只有当散射波谷量级小于0.02时，该波谷才与粮油作物种植前的平整土地相关。
[0043]
综上所述，本发明通过谐波拟合实现散射时序重构和去噪，通过幅度相位分析实现散射波形的检测和散射波谷的快速识别。现有的窗口统计法和局部拟合法难以确定具有适应性的滑动窗口和拟合函数，而差分搜索法需要复杂的时序逐点运算和二次判别。
附图说明
[0044]
图1为本发明的流程示意图；
[0045]
图2为实施例的时序谐波拟合图；
[0046]
图3为实施例的散射波形检测图；
[0047]
图4为实施例的粮油作物散射波谷识别图。
具体实施方式
[0048]
以下结合附图和实施例对本发明作进一步的详细说明。
[0049]
一种基于谐波幅度相位分析的粮油作物散射波谷识别方法(如图1所示)，本实施例开发环境为gee(google earth engine),编程语言为javascript。
[0050]
步骤1、检索sentinel-1在广西壮族自治区贵港市的2020年的升轨vh极化数据，使用imagecollection函数构建散射时序，使用image.linearregression函数实现谐波拟合，获得待定系数a,ai和
[0051]
步骤2、使用数学运算符add，subtract，multiply和divide，对以下参数进行逐像素计算：三种阶次余弦分量的波谷日期(v1，v
2a
，v
2b
，v
3a
，v
3b
和v
3c
)、各余弦项的幅度占比pi、余弦项两两之间的相位移动
[0052]
步骤3、使用逻辑运算符and，or，not，gt和lt，对pi和进行逐像素的数值分段，将符合不同数值区间的散射波形分为单波谷(s0，s1，s2和s3)型、双波谷(d0，d1，d2和d3)型和三波谷(t0)型。
[0053]
步骤4、根据散射波形的波谷类型逐像素地确定基准分量，使用逻辑运算符gt和lt搜索与基准分量波谷距离最近的其余两种余弦项中的波谷，使用数学运算符reduce(
‘
sum’)和reduce(
‘
mean’)计算波谷日期的加权平均值v
p
1，使用逻辑运算符gt和mask移除散射量级高于0.02的波谷。
[0054]
本实施例对2020年广西壮族自治区贵港市的sentinel-1数据进行处理，图2为实施例的时序谐波拟合图，图3为实施例的散射波形检测图，图4为实施例的粮油作物散射波谷识别图。
[0055]
通过以上实施例可见，本发明实现了2020年广西壮族自治区贵港市的粮油作物散射波谷识别，精度检验结果达到91.3％，本发明提供的方法在不预设滑动窗口尺寸和不假设局部拟合函数的情况下，灵活高效地实现了散射时序去噪，在不依赖复杂的时序逐点运算和二次判别的情况下，快速准确地识别了粮油作物散射波谷，为雷达数据的长时序大面积农业应用提供支撑。

技术特征：
1.一种基于谐波幅度相位分析的粮油作物散射波谷识别方法，其特征在于，具体步骤如下：步骤1.时序谐波拟合；对目标区域，获取雷达vh极化散射数据，构成散射时序s[t]，t代表时间，取值范围归一化为0到1，对应全年的首日到尾日；上式中a为常数项，代表散射年度均值；i依次取值1，2，3，代表余弦项的阶次；a
i
和分别依次为第i阶余弦项的幅度和相位，a、a
i
、的值未知，通过最小二乘拟合获得；步骤2.谐波指标计算；第i阶余弦项拥有i个波谷，波谷日期由计算得到；第1阶余弦项拥有1个波谷，设波谷日期为v1；第2阶余弦项拥有2个波谷，设波谷日期分别为v
2a
和v
2b
；第3阶余弦项拥有3个波谷，设波谷日期分别为v
3a
，v
3b
和v
3c
；计算第i阶余弦项幅度占全部余弦项幅度之和的比例：p
i
＝a
i
/(a1+a2+a3)计算第i阶余弦项相对于第j阶余弦项的相位移动：步骤3.散射波形检测；将散射波形分为三种类型：单波谷型、双波谷型和三波谷型；当p
i
和满足以下四种情况时，判定散射波形为单波谷型s型，具体包括：s0型：p1>2/3s1型：1/3≤p1≤2/3，p2>2p3，s2型：1/3≤p1≤2/3，p3>2p2，s3型：1/3≤p1≤2/3，p2≤2p3，p3≤2p2，当p
i
和满足以下四种情况时，判定散射波形为双波谷型d型，具体包括：d0型：p1<1/3，p2≥1/3或p3≤1/2d1型：1/3≤p1≤2/3，p2>2p3，d2型：1/3≤p1≤2/3，p3>2p2，
d3型：1/3≤p1≤2/3，p2≤2p3，p3≤2p2，或当p
i
和满足以下情况时，判定散射波形为三波谷型t型，即：t0型：p1<1/3，p2<1/3，p3>1/2；步骤4.粮油作物散射波谷识别：包括基准分量确定、最近波谷查找、加权平均计算和波谷量级限制四步完成；基准分量确定：对于单波谷型，基准分量为第1阶余弦项；对于双波谷型，基准分量为第2阶余弦项；对于三波谷型，基准分量为第3阶余弦项；基准分量的波谷数量与散射波形的波谷数量相同，且散射波形中任一波谷的日期v
p
，都可由与之距离最近的基准分量波谷的日期v
x
来近似表示；最近波谷查找：根据步骤2中波谷日期的计算方法，求解得到三个余弦项中所有波谷的日期；对于v
p
相应的v
x
，分别在其余两种余弦项中，查找与v
x
相应的基准分量波谷距离最近波谷的日期分别为v
y
和v
z
；至此，成功查找到与散射波形中目标波谷日期v
p
距离最近的三种余弦项中的波谷相应日期v
x
，v
y
和v
z
；加权平均计算：设v
x
、v
y
和v
z
对应余弦项的阶次i分别为i
x
、i
y
和i
z
，幅度a
i
分别为a
x
、a
y
和a
z
，则可以由v
x
，v
y
和v
z
加权平均计算得到v
p
1，v
p
1为v
p
的逼近值，公式如下：波谷量级限制：将散射波谷日期v
p
1作为t值代入步骤1中的谐波拟合公式以获得散射波谷量级，将0.02设置为波谷量级的宽松阈值，只保留量级低于0.02的散射波谷，这些波谷即为粮油作物散射波谷；至此，完成粮油作物散射波谷识别。2.如权利要求1所述基于谐波幅度相位分析的粮油作物散射波谷识别方法，其特征在于：所述步骤1中使用imagecollection函数构建散射时序，使用image.linearregression函数实现谐波拟合。

技术总结
本发明属于星载雷达数据处理与遥感制图技术领域，具体为一种基于谐波幅度相位分析的粮油作物散射波谷识别方法。本发明通过谐波拟合实现散射时序重构和去噪，通过幅度相位分析实现散射波形的检测和散射波谷的快速识别。本发明在不预设滑动窗口尺寸和不假设局部拟合函数的情况下，灵活高效地实现了散射时序去噪，在不依赖复杂的时序逐点运算和二次判别的情况下，快速准确地识别了粮油作物散射波谷，为雷达数据的长时序大面积农业应用提供支撑。有效克服现有的窗口统计法和局部拟合法难以确定具有适应性的滑动窗口和拟合函数，而差分搜索法需要复杂的时序逐点运算和二次判别，带来的相应行业问题。来的相应行业问题。来的相应行业问题。


技术研发人员：李世华 何泽
受保护的技术使用者：电子科技大学
技术研发日：2023.07.11
技术公布日：2023/10/19