CN105182362A - 一种自然水面偏振遥感成像仿真方法 - Google Patents

Abstract

本发明涉及一种自然水面偏振遥感成像仿真方法。其步骤如下：计算入射天空光的强度与偏振分布模式；基于Phillips波浪谱模型建立自然水面场景并得到波浪面元法向量；根据反射过程的矢量关系确定入射天空光的几何位置；计算天空光入射到微面元上的局部入射角；根据反射Mueller矩阵与入射光的偏振态计算反射光的Stokes参量，逐像素计算实现水面偏振遥感图像的仿真。本发明能够模拟不同时间、观测几何条件和不同风速下的自然水面偏振度和偏振角图像。本发明是基于物理模型和辐射传输模型推演的技术方法，具有模拟精度高、仿真速度快等优点，能够直观地了解海洋背景的偏振特性，可以为海洋目标偏振探测提供参考依据和偏振遥感器的设计提供验证数据。

Description

一种自然水面偏振遥感成像仿真方法

(一)所属技术领域

本发明涉及一种自然水面偏振遥感成像仿真方法，属于光学遥感探测领域，对海洋偏振探测的研究和在轨卫星偏振遥感器的设计具有重要意义和参考价值。

(二)背景技术

海洋探测与海上目标探测受到复杂波浪水面、太阳耀光等因素干扰，清晰成像尤为困难。波浪水面的反射光蕴含着大量的偏振信息，并且海上目标与海洋背景通常具有显著不同的偏振特性，这就为海洋目标探测提供了可靠依据。同时利用偏振手段进行探测可以有效地滤除太阳耀光。本发明中水面反射光的偏振特性研究对海洋偏振探测提供了理论依据。

虽然利用偏振特性进行遥感探测具有独特的优势，然而遥感探测特别是航天遥感成本高，传感器的研发制造难度大，阻碍目前偏振遥感成像技术的进展。借助计算机仿真技术与偏振遥感技术可以实现各种成像条件下地面场景的偏振遥感图像仿真。在偏振成像仿真方面，国内学者针对植被与土壤地表类型进行偏振特性的仿真建模研究。但迄今为止，尚未见到有效的波浪海面偏振遥感成像仿真方法的报道。本发明基于一种物理波浪谱模型，通过建立三维海面几何场景实现对成像背景的描述，建立场景中偏振反射物理模型，实现海面偏振遥感成像过程。因此本发明涉及一种自然水面偏振遥感成像仿真方法对偏振传感器的优化设计起到了指导作用并且大大节约了设计成本，可为偏振遥感应用提供有效的数据源。

(三)发明内容

本发明涉及一种自然水面偏振遥感成像仿真方法。技术解决方案是：首先，根据半经验瑞丽散射模型计算入射天空光的偏振分布模式；其次，基于Phillips波浪谱模型建立自然水面场景得到波浪面元法向量；根据反射过程的矢量辐射传输关系确定入射天空光的几何位置；计算入射到微面元的局部入射角；最后，根据菲涅耳反射定律，平静水面的Mueller矩阵乘以入射光的偏振态得到了反射光的Stokes参量，逐像素计算实现波浪水面偏振遥感图像的仿真。其具体步骤如下：

1一种自然水面偏振遥感成像仿真方法。其特征在于包含以下步骤：

(1)计算天空光的强度分布模式；

(2)计算天空光的偏振分布模式；

(3)基于Phillips波浪谱模型建立海洋高度场模型；

(4)计算网格点的波浪微面元法线方向；

(5)入射光与微面元形成的局部入射角；

(6)建立气-水界面的偏振反射率模型；

(7)计算波浪水面反射光的偏振态；

(8)模拟自然水面的偏振遥感图像。

2步骤(1)中所述的“计算天空光的强度分布模式”，建立晴朗天空光辐亮度分布模型，天空某观测方向的散射光强度可以表达为：

$N(\mathrm{\γ},{\mathrm{\θ}}_s,\mathrm{\θ})=(A+{\mathrm{Be}}^{-m\mathrm{\γ}}+{\mathrm{Ccos}}^2\mathrm{\γ}cos\mathrm{\θ})(1-e^{-{\mathrm{\ρsec\θ}}_s})(1-e^{-T\sec \mathrm{\θ}})$

其中，N为天空光的辐亮度，θ_s为太阳天顶角，θ为天空某位置的天顶角，ρ为大气单次散射反照率，τ为大气光学厚度，m为大气质量，γ为散射角。

3步骤(2)中所述的“计算天空光的偏振分布模式”，具体技术流程如下：

第一步：根据公式计算散射角γ，其中θ_s为太阳天顶角，θ为观测天顶角，φ为相对方位角；

第二步：计算天空光偏振度，如下式所示：

$\mathrm{\δ}={\mathrm{\δ}}_{max}\frac{{\sin }^2\mathrm{\γ}}{1+{\cos }^2\mathrm{\γ}}$

δ_max为特定太阳高度角时的最大偏振度，最大值取1；

第三步：计算天空光的偏振角，设散射面的法向量入射光矢量那么散射角aop＝|cos^-1m-90|，

其中 $m=(x^{\′}{z}_2-z^{\′}{x}_2)/\left(\sqrt{x_2^2+z_2^2}\sqrt{x^{\′2}+y^{\′2}+z^{\′2}}\right).$

4步骤(3)中所述的“基于Phillips波浪谱模型建立海洋高度场模型”，具体技术流程如下：

第一步：建立Phillips波浪谱模型，

$P_h\left(\stackrel{\→}{k}\right)=A\frac{\exp (-1/{\left(kl\right)}^2)}{k^4}{\left|\stackrel{\→}{k}\hat{w}\right|}^2\·\exp (-k^2{l}^2)$

其中，k为波数向量，w为风速向量，A，l为常量。

第二步：计算高度场的傅立叶波谱值，

$\tilde{h}\left(\right|\hat{k}|,0)=\frac{1}{\sqrt{2}}({\mathrm{\ϵ}}_r+{\mathrm{i\ϵ}}_i)\sqrt{P_h\left(\right|\hat{k}\left|\right)}$

其中ε_r、ε_i是两个相互独立的高斯随机数。

第三步：计算波浪频谱在时间t时的高度场，

$\tilde{h}(k,t)={\tilde{h}}_0\left(k\right)\exp \[iw\left(kt\right)\]+{{\tilde{h}}_0}^{*}(-k)\exp \[-iw\left(kt\right)\]$

第四步：计算t时刻任意水平位置X＝(x,z)上波浪水面的高度：

$\tilde{h}(x,t)=\mathrm{Re}\[\underset{k}{\Σ}\tilde{h}(k,t)\exp (ik\·x)\]$

上式是一个频域变换到空间域的过程。

5步骤(4)中所述的“计算网格点的波浪微面元法线方向”，利用临近网格点的有限差方法，具体技术流程如下：

第一步：面元斜坡在x，y，z方向的分量分别为：

$\begin{array}{c}\frac{dh}{dx}=\left(h\right(x,t\left)\right(m-1,n+1)+2\·h(x,t\left)\right(m,n+1)+h(x,t\left)\right(m+1,n+1\left)\right)\\ -\left(h\right(x,t\left)\right(m-1,n-1)+2\·h(x,t\left)\right(m,n-1)+h(x,t\left)\right(m+1,n-1\left)\right)\end{array}$

$\begin{array}{c}\frac{dh}{dy}=\left(h\right(x,t\left)\right(m+1,n-1)+2\·h(x,t\left)\right(m+1,n)+h(x,t\left)\right(m+1,n+1\left)\right)\\ -\left(h\right(x,t\left)\right(m-1,n-1)+2\·h(x,t\left)\right(m-1,n)+h(x,t\left)\right(m-1,n+1\left)\right)\end{array}$

$\frac{dh}{dz}=4\·\frac{Lx}{M}$

第二步：每一个网格点的法向量为：

$\mathrm{\η}=\sqrt{{\left(\frac{dh}{dx}\right)}^2+{\left(\frac{dh}{dy}\right)}^2+{\left(\frac{dh}{dz}\right)}^2}$

6步骤(5)中所述的“入射光与微面元形成的局部入射角”，具体技术流程如下：

第一步：入射光线ξ_i、面元法线n、反射光线ξ_r满足下式关系，且有反射角θ_r等于入射角θ_i；

${\mathrm{\ξ}}_r={\mathrm{\ξ}}_i-2({\mathrm{\ξ}}_i\·\stackrel{\→}{n})\·\stackrel{\→}{n}$

${\mathrm{\θ}}_r={\cos }^{-1}\left|{\mathrm{\ξ}}_i\·\stackrel{\→}{n}\right|={\mathrm{\θ}}_i$

第二步：确定入射天空光方向，并用空间三维坐标表示：

其中入射天顶角θ_i∈(0,π/2)，入射方位角

第三步：考虑到面元斜坡与水平面的倾角，入射到波浪微面元上的局部入射角的余弦为：

当天顶观测时，出射光矢量ξ_r＝(0,0,1)，则入射光天顶角为θ_i”＝2θ_i'。

7步骤(6)中所述的“建立气-水界面的偏振反射率模型”，具体技术流程如下：

第一步：平静水面反射Mueller矩阵为：

$M\left(\mathrm{\ω}\right)=\frac{1}{2}\left(\begin{array}{cccc}{r}_s^2+r_p^2& r_p^2-r_s^2& 0& 0\\ r_p^2-r_s^2& r_s^2+r_p^2& 0& 0\\ 0& 0& 2r_s{r}_p& 0\\ 0& 0& 0& 2r_s{r}_p\end{array}\right)$

其中r_s、r_p为反射E矢量的垂直分量和平行分量；

第二步：旋转参考面，旋转矩阵为：

$C\left(i\right)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos 2i& sin2i& 0\\ 0& -sin2i& \cos 2i& 0\\ 0& 0& 0& 1\end{array}\right]$

其中i为参考面旋转角度；

第三步：根据第一步和第二步求得面元的偏振反射率矩阵：

8步骤(7)中所述的“计算波浪水面反射光的偏振态”，具体技术流程如下：

第一步：用Stokes矢量表示反射光在某一点的偏振态：

S＝(I,Q,U,V)^T

其中I、Q、U和V为斯托克斯参量；

第二步：反射光偏振特性的Stokes矢量就可以表示为，

S_r＝M·S_i

第三步：根据Stokes参量计算面元反射光的偏振态，偏振度d和偏振角α由下式计算：

$\left\{\begin{array}{c}d=\frac{\sqrt{Q^2+U^2+V^2}}{I}\\ \mathrm{\α}=\frac{1}{2}ta{n}^{-1}\left(\frac{U}{Q}\right)\end{array}\right.$

9步骤(8)中所述的“模拟自然水面的偏振遥感图像”，具体技术流程：根据设定的海域面积和网格点数，逐像素计算每一微面元反射光的偏振态，利用C++/OpenGL实现图像的显示。

本发明与现有技术相比的优点在于：

(1)目前学术界缺少有效的波浪海面偏振遥感成像仿真的技术方法，本发明具有显著的创新性，拓展了偏振光学遥感领域的研究。

(2)本发明实现了场景—图像的仿真过程，建立了精确的波浪水面场景，有利于仿真结果更加真实准确。

(3)本发明逐像素定量计算了入射光经波浪海面反射时偏振特性的变化，可以模拟任意成像条件下的遥感图像，相比于其他仿真方法，本发明计算速度快，参数可调，适应性强。

(四)附图说明

图1为本发明的总体技术流程图。图2为基于本发明的方法仿真的不同成像条件下的自然水面的偏振图像，其中，

图2-(a)：ωs＝5m/s，θs＝30°垂直观测时的偏振度图像；

图2-(b)：ωs＝5m/s，θs＝60°垂直观测时的偏振度图像；

图2-(c)：ωs＝5m/s，θs＝40°，θv＝10°观测时的偏振角图像；

图2-(d)：ωs＝10m/s，θs＝40°，θv＝10°观测时的偏振角图像。

(五)具体实施方式

为了更好地说明本发明涉及的一种自然水面偏振遥感成像仿真方法，利用本发明的方法进行了仿真，取得了良好的效果，具体实施步骤如下：

(1)计算天空光的强度分布模式；

(2)计算天空光的偏振分布模式；

(3)基于Phillips波浪谱模型建立海洋高度场模型；

(4)计算网格点的波浪微面元法线方向；

(5)入射光与微面元形成的局部入射角；

(6)建立气-水界面的偏振反射率模型；

(7)计算波浪水面反射光的偏振态；

(8)模拟自然水面的偏振遥感图像。

Claims (9)

1.一种自然水面偏振遥感成像仿真方法。其特征在于包含以下步骤：

(1)计算天空光的强度分布模式；

(2)计算天空光的偏振分布模式；

(3)基于Phillips波浪谱模型建立海洋高度场模型；

(4)计算网格点的波浪微面元法线方向；

(5)入射光与微面元形成的局部入射角；

(6)建立气-水界面的偏振反射率模型；

(7)计算波浪水面反射光的偏振态；

(8)模拟自然水面的偏振遥感图像。

2.根据权利1要求所述的一种自然水面偏振遥感成像仿真方法，其特征在于：步骤(1)中所述的“计算天空光的强度分布模式”，建立晴朗天空光辐亮度分布模型，天空某观测方向的散射光强度可以表达为：

$N(\mathrm{\γ},{\mathrm{\θ}}_s,\mathrm{\θ})=(A+{\mathrm{Be}}^{-m\mathrm{\γ}}+{\mathrm{Ccos}}^2\mathrm{\γ}cos\mathrm{\θ})(1-e^{-{\mathrm{\ρsec\θ}}_s})(1-e^{-T\sec \mathrm{\θ}})$

其中，N为天空光的辐亮度，θ_s为太阳天顶角，θ为天空某位置的天顶角，ρ为大气单次散射反照率，τ为大气光学厚度，m为大气质量，γ为散射角。

3.根据权利要求1所述的一种自然水面偏振遥感成像仿真方法，其特征在于：步骤(2)中所述的“计算天空光的偏振分布模式”，具体技术流程如下：

第一步：根据公式计算散射角γ，其中θ_s为太阳天顶角，θ为观测天顶角，φ为相对方位角；

第二步：计算天空光偏振度，如下式所示：

$\mathrm{\δ}={\mathrm{\δ}}_{max}\frac{{\sin }^2\mathrm{\γ}}{1+{\cos }^2\mathrm{\γ}}$

δ_max为特定太阳高度角时的最大偏振度，最大值取1；

第三步：计算天空光的偏振角，设散射面的法向量入射光矢量那么散射角aop＝|cos^-1m-90|，

其中 $m=(x^{\′}{z}_2-z^{\′}{x}_2)/\left(\sqrt{x_2^2+z_2^2}\sqrt{x^{\′2}+y^{\′2}+z^{\′2}}\right).$

4.根据权利要求1所述的一种自然水面偏振遥感成像仿真方法，其特征在于：步骤(3)中所述的“基于Phillips波浪谱模型建立海洋高度场模型”，具体技术流程如下：

第一步：建立Phillips波浪谱模型，

$P_h\left(\stackrel{\→}{k}\right)=A\frac{\exp (-1/{\left(kl\right)}^2)}{k^4}{\left|\stackrel{\→}{k}\hat{w}\right|}^2\·\exp (-k^2{l}^2)$

其中，k为波数向量，w为风速向量，A，l为常量。

第二步：计算高度场的傅立叶波谱值，

$\tilde{h}\left(\right|\hat{k}|,0)=\frac{1}{\sqrt{2}}({\mathrm{\ϵ}}_r+{\mathrm{i\ϵ}}_i)\sqrt{P_h\left(\right|\hat{k}\left|\right)}$

其中ε_r、ε_i是两个相互独立的高斯随机数。

第三步：计算波浪频谱在时间t时的高度场，

$\tilde{h}(k,t)={\tilde{h}}_0\left(k\right)\exp \[iw\left(kt\right)\]+{{\tilde{h}}_0}^{*}(-k)\exp \[-iw\left(kt\right)\]$

第四步：计算t时刻任意水平位置X＝(x,z)上波浪水面的高度：

$\tilde{h}(x,t)=\mathrm{Re}\[\underset{k}{\Σ}\tilde{h}(k,t)\exp (ik\·x)\]$

上式是一个频域变换到空间域的过程。

5.根据权利要求1所述的一种自然水面偏振遥感成像仿真方法，其特征在于：步骤(4)中所述的“计算网格点的波浪微面元法线方向”，利用临近网格点的有限差方法，具体技术流程如下：

第一步：面元斜坡在x，y，z方向的分量分别为：

$\begin{array}{c}\frac{dh}{dx}=(h(x,t)(m-1,n+1)+2\·h(x,t)(m,n+1)+h(x,t)(m+1,n+1))\\ -\left(h\right(x,t)(m-1,n-1)+2\·h(x,t)(m,n-1)+h(x,t)(m+1,n-1))\end{array}$

$\begin{array}{c}\frac{dh}{dy}=(h(x,t)(m+1,n-1)+2\·h(x,t)(m+1,n)+h(x,t)(m+1,n+1))\\ -\left(h\right(x,t)(m-1,n-1)+2\·h(x,t)(m-1,n)+h(x,t)(m-1,n+1))\end{array}$

$\frac{dh}{dz}=4\·\frac{Lx}{M}$

第二步：每一个网格点的法向量为：

$\mathrm{\η}=\sqrt{{\left(\frac{dh}{dx}\right)}^2+{\left(\frac{dh}{dy}\right)}^2+{\left(\frac{dh}{dz}\right)}^2}$

6.根据权利要求1所述的一种自然水面偏振遥感成像仿真方法，其特征在于：步骤(5)中所述的“入射光与微面元形成的局部入射角”，具体技术流程如下：

第一步：入射光线ξ_i、面元法线n、反射光线ξ_r满足下式关系，且有反射角θ_r等于入射角θ_i；

${\mathrm{\ξ}}_r={\mathrm{\ξ}}_i-2({\mathrm{\ξ}}_i\·\stackrel{\→}{n})\·\stackrel{\→}{n}$

${\mathrm{\θ}}_r={\cos }^{-1}\left|{\mathrm{\ξ}}_i\·\stackrel{\→}{n}\right|={\mathrm{\θ}}_i$

第二步：确定入射天空光方向，并用空间三维坐标表示：

其中入射天顶角θ_i∈(0,π/2)，入射方位角

第三步：考虑到面元斜坡与水平面的倾角，入射到波浪微面元上的局部入射角的余弦为：

当天顶观测时，出射光矢量ξ_r＝(0,0,1)，则入射光天顶角为θ_i”＝2θ_i'。

7.根据权利要求1所述的一种自然水面偏振遥感成像仿真方法，其特征在于：步骤(6)中所述的“建立气-水界面的偏振反射率模型”，具体技术流程如下：

第一步：平静水面反射Mueller矩阵为：

$M\left(\mathrm{\ω}\right)=\frac{1}{2}\left(\begin{array}{cccc}{r}_s^2+r_p^2& r_p^2-r_s^2& 0& 0\\ r_p^2-r_s^2& r_s^2+r_p^2& 0& 0\\ 0& 0& 2r_s{r}_p& 0\\ 0& 0& 0& 2r_s{r}_p\end{array}\right)$

其中r_s、r_p为反射E矢量的垂直分量和平行分量；

第二步：旋转参考面，旋转矩阵为：

$C\left(i\right)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos 2i& sin2i& 0\\ 0& -sin2i& \cos 2i& 0\\ 0& 0& 0& 1\end{array}\right]$

其中i为参考面旋转角度；

第三步：根据第一步和第二步求得面元的偏振反射率矩阵：

8.根据权利要求1所述的一种自然水面偏振遥感成像仿真方法，其特征在于：步骤(7)中所述的“计算波浪水面反射光的偏振态”，具体技术流程如下：

第一步：用Stokes矢量表示反射光在某一点的偏振态：

S＝(I,Q,U,V)^T

其中I、Q、U和V为斯托克斯参量；

第二步：反射光偏振特性的Stokes矢量就可以表示为，

S_r＝M·S_i

第三步：根据Stokes参量计算面元反射光的偏振态，偏振度d和偏振角α由下式计算：

$\left\{\begin{array}{c}d=\frac{\sqrt{Q^2+U^2+V^2}}{I}\\ \mathrm{\α}=\frac{1}{2}ta{n}^{-1}\left(\frac{U}{Q}\right)\end{array}\right.$

9.根据权利要求1所述的一种自然水面偏振遥感成像仿真方法，其特征在于：步骤(8)中所述的“模拟自然水面的偏振遥感图像”，具体技术流程：根据设定的海域面积和网格点数，逐像素计算每一微面元反射光的偏振态，利用C++/OpenGL实现图像的显示。