CN103512663A - 一种起伏月面微波辐射亮温的计算方法 - Google Patents

Abstract

本发明公开了一种起伏月面微波辐射亮温的计算方法，包括以下步骤：获取月球表面的微波数据参数，并利用平面拟合、坐标转换以及遮蔽函数对月球表面的微波数据参数进行计算，以获得月球表面有效太阳辐照度，根据月壤的物理参数和步骤（1）中获得的有效太阳辐照度并使用热传导理论和月壤分层模型获得月壤不同深度的温度T，根据步骤（2）获得的温度剖面并利用burke多层平面分层亮温模型以及电磁波极化理论计算月球表面亮温。本发明能够解决现有模型中忽略地形起伏对月壤微波辐射亮温影响这一技术问题。

Description

一种起伏月面微波辐射亮温的计算方法

技术领域

本发明属于微波遥感领域，更具体地，涉及一种起伏月面微波辐射亮温的计算方法。

背景技术

“嫦娥一号”（CE-1）和“嫦娥二号”(CE-2)对月探测的一个重要科学任务就是记录月壤微波辐射亮度温度数据，从而反演月壤层厚度，并在此基础上评估月球上氦3资源含量。月壤厚度的研究是月壤研究的一项重要的内容，它对以后的月球探测、载人登月、月球基地选址以及开发和利用月球资源等都具有非常重要的意义。要想提高反演月壤厚度的精度，首先要建立一个准确的月壤微波辐射亮温正向模型。

目前对粗糙月面微波辐射亮温的研究，主要是基于平面分层媒质，并未考虑地形起伏对微波辐射亮温的影响。事实上，对于星载遥感应用，微波辐射计的探测范围在几十公里以上，其中的月表可能有地形起伏。这种大尺度的地形起伏将直接导致太阳入射角的变化，引起区域内月球表面物理温度的升高和降低。对于月球表面存在的小尺度地形起伏（如辐射计探测范围内相对于基准面的高度起伏），因为随机粗糙面的遮蔽效应会导致太阳光的照射情况的变化，从而引起月球表面物理温度的降低。因此有必要对粗糙月面的地形起伏加以妥善处理，建立更符合实际情况的粗糙月面分层微波辐射亮温模型。

在地形起伏较大的区域，如月陆地区，地形起伏对于月球表面温度的影响是不能忽略的，因为这些区域的倾斜、遮蔽等特性对太阳光照、有效太阳辐照度、天线观测角等有重要的影响，进而对表面温度有较大影响，因此现有的平面模型会导致月球表面亮温的计算值与实测值相差较大（可达十几K），从而造成反演月壤厚度的精确度较低。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种起伏月面微波辐射亮温的计算方法，其目的在于解决现有模型中忽略地形起伏对月壤微波辐射亮温影响这一技术问题。

为实现上述目的，按照本发明的一个方面，提供了一种起伏月面微波辐射亮温的计算方法，包括以下步骤：

（1）获取月球表面的微波数据参数，并利用平面拟合、坐标转换以及遮蔽函数对月球表面的微波数据参数进行计算，以获得月球表面有效太阳辐照度；本步骤包括以下子步骤：

（1-1）通过微波辐射计获取月球表面被测点的微波数据参数，包括经度lgt、纬度lat、高程H、太阳入射角θ、太阳方位角φ、对应的月球时间t、月表接收的太阳辐照度I₀、月壤中氧化铁和二氧化钛的含量S以及微波辐射计的工作频率f，这里的高程H和太阳入射角θ、太阳方位角φ都是在全局坐标系下得到的；

（1-2）根据步骤（1-1）中获得的微波数据参数并使用最小二乘法将微波辐射计天线的月面足印内小面元中所获取的所有点拟合成一个平面，并获取该平面的倾斜角和方位角；

（1-3）依据步骤（1-1）和（1-2）中获得的参数，以步骤（1-2）中的拟合平面作为基准面建立本地坐标系，根据坐标转换获得太阳本地入射角θ′和太阳本地方位角φ′，其计算公式为：

cosθ′＝cosαcosθ+sinαsinθcos(φ-β)

sinθ′＝(1-cos²θ′)^1/2

sinφ′＝sinθsin(φ-β)/sinθ′

cosφ′＝(cosαsinθcos(φ-β)-sinαcosθ)/sinθ′

（1-4）依据步骤（1-1）至（1-3）中获得的参数获得微波辐射计探测范围内的月球表面斜率p、q以及月球表面的局部法向矢量

表面均方根高度σ、表面相关长度l_x和l_y、表面均方根斜率w；

(1-5)依据步骤（1-3）和（1-4）中获得的参数并根据遮蔽函数理论获得微波辐射计探测范围内的月球表面的二维双向遮蔽函数S(p,q,θ′,φ′,θ₁)，其中观测角θ₁＝α；

（1-6）依据步骤（1-1）至（1-5）获得的参数计算太阳辐照度的衰减因子IN和太阳有效辐照度I_eff，其计算公式为：

$\mathrm{IN}=\∫\∫\mathrm{dpdq}\frac{1}{{2\mathrm{\πw}}_x{w}_y}\exp (-\frac{p^2}{{2w}_x^2}-\frac{q^2}{{2w}_y^2})\·S(p,q,{\mathrm{\θ}}^{\′},{\mathrm{\φ}}^{\′},{\mathrm{\θ}}_1)\·(p{\sin \mathrm{\θ}}^{\′}\cos {\mathrm{\φ}}^{\′}+q\sin {\mathrm{\θ}}^{\′}\sin {\mathrm{\φ}}^{\′}+\cos {\mathrm{\θ}}^{\′})$

I_eff＝I₀·IN

（2）根据月壤的物理参数和步骤（1）中获得的有效太阳辐照度并使用热传导理论和月壤分层模型获得月壤不同深度的温度T；

（2-1）将月球表面的月壤层划分成N层，第二层到第N层的厚度满足以下等式：

d_m＝A·e^B·m

其中m表示层数，e为自然对数函数的底数，A和B均为常数，其取值与月壤的厚度d成正比；

（2-2）获取月球表面物理参数的模型，月壤的物理参数包括月壤密度ρ、月壤的复介电常数ε^*、月壤的比热C和月壤的导热率K；

（2-3）根据步骤（2-2）获得的参数并根据热传导理论求解热传导方程，得到月壤不同深度的温度；

（3）根据步骤（2）获得的温度剖面并利用burke多层平面分层亮温模型以及电磁波极化理论计算月球表面亮温。

优选地，步骤（1-2）中，曲面拟合平面的具体方法如下：

（1-2-1）将小面元中n个点拟合成的平面z的表达式为

z＝b₁+b₂x+b₃y

其中，b₁、b₂和b₃为待求系数，则有

$\left(\begin{array}{c}{H}_1\\ H_2\\ \·\\ \·\\ \·\\ H_n\end{array}\right)=b_1\·\left(\begin{array}{c}1\\ 1\\ \·\\ \·\\ \·\\ 1\end{array}\right)+b_2\·\left(\begin{array}{c}{\mathrm{lgt}}_1\\ {\mathrm{lgt}}_2\\ \·\\ \·\\ \·\\ {\mathrm{lgt}}_n\end{array}\right)+b_3\·\left(\begin{array}{c}{\mathrm{lat}}_1\\ {\mathrm{lat}}_2\\ \·\\ \·\\ \·\\ {\mathrm{lat}}_n\end{array}\right)$

（1-2-2）拟合平面z＝b₁+b₂x+b₃y与x轴、y轴、z轴的交点分别为

a₃(0,0,b₁)，

为拟合平面的外法向矢量，则有

$\left\{\begin{array}{c}{\stackrel{\→}{n}}_1\·\stackrel{\→}{a_1{a}_3}=0\\ {\stackrel{\→}{n}}_1\·\stackrel{\→}{a_2{a}_3}=0\end{array}\right.$

（1-2-3）解得 ${\stackrel{\→}{n}}_1=({-b}_2\·z,-b_3\·z,z),$ 即 ${\stackrel{\→}{n}}_1=(-b_2,{-b}_3,1);$

（1-2-4）获得拟合平面的倾斜角α和方位角β为：

$\cos \mathrm{\α}=\frac{1}{\sqrt{b_2^2+b_3^2+1}}$

$\tan \mathrm{\β}=\frac{b_3}{b_2}$

优选地，步骤（1-4）的计算公式为：

$p=\frac{\∂H}{\∂x}$

$q=\frac{\∂H}{\∂y}$

${\stackrel{\→}{n}}_2=\frac{-p\hat{x}-q\hat{y}+\hat{z}}{\sqrt{1+p^2+q^2}}$

其中

分别为x轴、y轴和z轴的单位方向矢量；

表面均方根高度σ为：

$\mathrm{\σ}={\left[\frac{1}{M_1\·M_2}(\underset{m_1=1}{\overset{M_1}{\mathrm{\Σ}}}\underset{m_2=1}{\overset{M_2}{\mathrm{\Σ}}}{H}_{m_1{m}_2}^2-M_1\·M_2\·{\stackrel{\‾}{H}}^2)\right]}^{1/2}$

其中

$\stackrel{\‾}{H}=\frac{1}{M_1\·M_2}\underset{m_1=1}{\overset{M_1}{\mathrm{\Σ}}}\underset{m_2=1}{\overset{M_2}{\mathrm{\Σ}}}{H}_{m_1{m}_2}$

M₁，M₂分布为微波辐射计探测范围内的x和y方向的高程数目；

$l_x=\frac{1}{M_1}\underset{m_1=1}{\overset{M_1}{\mathrm{\Σ}}}{l}_{{\mathrm{xm}}_1}$

$l_y=\frac{1}{M_2}\underset{m_2=1}{\overset{M_2}{\mathrm{\Σ}}}{l}_{{\mathrm{ym}}_2}$

对于微波辐射计探测范围内的二维表面，方向和y方向的均方根斜率分别为：

$w_x=\sqrt{2}\frac{\mathrm{\σ}}{l_x}$

$w_y=\sqrt{2}\frac{\mathrm{\σ}}{l_y}$

二维均方根斜率为：

$w=\sqrt{{(w_x\·{\cos \mathrm{\φ}}^{\′})}^2+{(w_y\·\sin {\mathrm{\φ}}^{\′})}^2}$

优选地，步骤（1-5）的计算公式为：

当θ₁＞θ′时，

$S(p,q,{\mathrm{\θ}}^{\′},{\mathrm{\φ}}^{\′},{\mathrm{\θ}}_1)=\frac{\mathrm{\γ}(\stackrel{\‾}{\mathrm{\μ}}-p\sin {\mathrm{\φ}}^{\′}-q\cos {\mathrm{\φ}}^{\′})}{[\mathrm{\Λ}\left(\stackrel{\‾}{\mathrm{\μ}}\right)+1]}$

其中

$\mathrm{\&gamma;}(\stackrel{\&OverBar;}{\mathrm{\&mu;}}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})=\left\{\begin{array}{cc}1,& \stackrel{\&OverBar;}{\mathrm{\&mu;}}\&GreaterEqual;p\sin {\mathrm{\&phi;}}^{\&prime;}+{q\cos \mathrm{\&phi;}}^{\&prime;}\\ 0,& \stackrel{\&OverBar;}{\mathrm{\&mu;}}<p\sin {\mathrm{\&phi;}}^{\&prime;}+q\cos {\mathrm{\&phi;}}^{\&prime;}\end{array}\right.,$ $\mathrm{\&Lambda;}\left(\stackrel{\&OverBar;}{\mathrm{\&mu;}}\right)=\frac{1}{2}\&CenterDot;[{\left(\frac{2}{\mathrm{\&pi;}}\right)}^{1/2}\&CenterDot;\frac{w}{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{e}^{-{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}^2/{2w}^2}-\mathrm{erfc}\left(\frac{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{\sqrt{2}w}\right)],$ erfc是误差函数的余集函数，定义为：

$\mathrm{erfc}\left(\frac{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{\sqrt{2}w}\right)=\frac{2}{\sqrt{\mathrm{\&pi;}}}{\&Integral;}_{\frac{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{\sqrt{2}w}}^{\&infin;}{e}^{-x^2}\mathrm{dx}$

当0＜θ₁＜θ′时：

$S(p,q,{\mathrm{\&theta;}}^{\&prime;},{\mathrm{\&phi;}}^{\&prime;},{\mathrm{\&theta;}}_1)=\frac{\mathrm{\&gamma;}(\mathrm{\&mu;}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})}{[\mathrm{\&Lambda;}\left(\mathrm{\&mu;}\right)+1]}$

其中

$\mathrm{\&Lambda;}\left(\mathrm{\&mu;}\right)=\frac{1}{2}\&CenterDot;[{\left(\frac{2}{\mathrm{\&pi;}}\right)}^{1/2}\&CenterDot;\frac{w}{\mathrm{\&mu;}}{e}^{-{\mathrm{\&mu;}}^2/{2w}^2}-\mathrm{erfc}\left(\frac{\mathrm{\&mu;}}{\sqrt{2}w}\right)],$ $\mathrm{\&gamma;}(\mathrm{\&mu;}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})=\left\{\begin{array}{cc}1,& \mathrm{\&mu;}\&GreaterEqual;p\sin {\mathrm{\&phi;}}^{\&prime;}+{q\cos \mathrm{\&phi;}}^{\&prime;}\\ 0,& \mathrm{\&mu;}<p\sin {\mathrm{\&phi;}}^{\&prime;}+q\cos {\mathrm{\&phi;}}^{\&prime;}\end{array}\right.,$

当θ₁＜0时：

$S(p,q,{\mathrm{\&theta;}}^{\&prime;},{\mathrm{\&phi;}}^{\&prime;},{\mathrm{\&theta;}}_1)=\frac{\mathrm{\&gamma;}(\stackrel{\&OverBar;}{\mathrm{\&mu;}}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})\&CenterDot;\mathrm{\&gamma;}(\mathrm{\&mu;}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})}{[\mathrm{\&Lambda;}\left(\stackrel{\&OverBar;}{\mathrm{\&mu;}}\right)+1]\&CenterDot;[\mathrm{\&Lambda;}\left(\mathrm{\&mu;}\right)+1]}$

其中μ＝cotθ′， $\mathrm{\&Lambda;}\left(\mathrm{\&mu;}\right)=\frac{1}{2}\&CenterDot;[{\left(\frac{2}{\mathrm{\&pi;}}\right)}^{1/2}\&CenterDot;\frac{w}{\mathrm{\&mu;}}{e}^{-{\mathrm{\&mu;}}^2/{2w}^2}-\mathrm{erfc}\left(\frac{\mathrm{\&mu;}}{\sqrt{2}w}\right)]$ $\mathrm{\&Lambda;}\left(\stackrel{\&OverBar;}{\mathrm{\&mu;}}\right)=\frac{1}{2}\&CenterDot;[{\left(\frac{2}{\mathrm{\&pi;}}\right)}^{1/2}\&CenterDot;\frac{w}{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{e}^{-{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}^2/{2w}^2}-\mathrm{erfc}\left(\frac{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{\sqrt{2}w}\right)].$

优选地，步骤（2-2）包括以下子步骤：

（2-2-1）根据登月点实测和经验公式获得月壤密度ρ，具体采用以下公式：

$\mathrm{\&rho;}=1.92\frac{z+12.2}{z+18}$

其中ρ是月壤密度，z是月壤深度；

（2-2-2）根据月球样品的回归分析获得复介电常数ε^*，其余与月壤密度ρ和氧化铁与二氧化钛的含量S的关系如下：

ε^*＝ε₀（ε′-jε″）

ε′＝1.919^ρ

ε″＝ε′·10^{0.038S+0.312ρ-3.260}

其中，ε₀＝0.8854F/m， ε′表示复介电常数的实部，ε″表示复介电常数的虚部，S表示氧化铁与二氧化钛的含量；

（2-2-3）根据实验获得月壤比热与月壤温度关系为：

C＝c₁T³+c₂T²+c₃T+c₄

其中T为月壤的温度。

（2-2-4）根据Vasavada模型获得月壤热导率与月壤温度的关系为：

$K=K_c\&CenterDot;(1+\mathrm{\&chi;}\&CenterDot;{\left(\frac{T}{T_{350}}\right)}^3)$

其中K_c表示固体导热率，χ是辐射热导率与固体热导率的比值。

优选地，步骤（2-3）具体为：根据能量守恒定律导出的一维热传导方程为：

$\mathrm{\&rho;}(z,T)C(z,T)\frac{\&PartialD;T}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;z}\left[K(z,T)\frac{\&PartialD;T}{\&PartialD;z}\right]+Q(z,t)$ 其中，ρ(z,T)表示密度（kg/m²），C(z,T)表示比热(J/(kg·K))，K(z,T)表示热导率(W/(m·K))，Q(z,t)表示部分透明介质由于吸收太阳的辐射而产生的源项，且上述方程满足以下边界条件：

（2-3-1）在月壤表层：

$K_s\frac{\&PartialD;T}{\&PartialD;z}{|}_s=\mathrm{\&epsiv;}{\mathrm{\&sigma;}}_B{T}_s^4-(1-A_b)[I_{\mathrm{eff}}+E]+J_0$

其中,

表示表面处的温度梯度，K_s为表面热导率，表示表示传入次表面的能量；ε为红外表面发射率(一般设定为0.90-1.0之间)，σ_B为Stefan-Boltzman常数5.6703×10^-8Wm^-2·K^-4，T_s为月壤表层温度，

表示月球表面辐射的红外能量；A_b为月表热辐射反照度0.12，E表示地球反射的太阳辐照度；J₀表示月球内部发射的热通亮；

（2-3-2）在热平衡深度Z₀

$\frac{\&PartialD;T}{\&PartialD;z}{|}_{\mathrm{depth}}=\frac{-J_0}{K_{\mathrm{depth}}}<<1$

其中，K_depth为在Z₀时的温度梯度，J₀是远小于1的常数。

优选地，步骤（3）包括以下子步骤：

（3-1）根据步骤（2）中获得月壤第i层的物理温度为T_i，复介电常数为

获得月壤第i层的复波阻抗η_i和复波数k_i，计算公式如下：

${\mathrm{\&eta;}}_i=\sqrt{\frac{{\mathrm{\&mu;}}_0}{{\mathrm{\&epsiv;}}_i^{*}}},$ $k_i=2\mathrm{\&pi;f}\sqrt{{\mathrm{\&mu;}}_0{\mathrm{\&epsiv;}}_i^{*}}$

其中μ₀＝4π×10^-7H/m；

（3-2）根据Burke推导出的具有不同温度、介电常数的多层平面分层结构辐射亮温的一阶近似解，获得辐射计在观测角θ₁接收到的亮温如下：

${\mathrm{TB}}_r=(1-{\mathrm{\&Gamma;}}_{1,r}\left({\mathrm{\&theta;}}_1\right))\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{T}_{i+1}(1-\frac{1}{L_{i+1}})(1+\frac{{\mathrm{\&Gamma;}}_{i+1,r}}{L_{i+1}})\underset{m_3=2}{\overset{i}{\mathrm{\&Pi;}}}\left(\frac{1-{\mathrm{\&Gamma;}}_{m_3,r}}{L_{m_3}}\right)$

其中，r表示极化状态，可取h或v极化，即由上式可以得到h极化的亮温TB_h和v极化的亮温TB_v，T_i+1是第i+1层的物理温度，

是辐射在第i+1层的功率损耗因子，

为功率吸收系数，

为第i+1层的复介电常数，d_i+1为第i+1层的厚度，θ_i+1是传播光线在第i+1层的方向角，满足

Γ_i+1是第i+1层边界上的菲涅尔反射率，Γ_i+1,r＝|R_i+1,r|²，在水平极化和铅垂极化下，分别为：

$R_{i+1,h}=\frac{{\mathrm{\&eta;}}_{i+2}\cos {\mathrm{\&theta;}}_i-{\mathrm{\&eta;}}_{i+1}\cos {\mathrm{\&theta;}}_{i+1}}{{\mathrm{\&eta;}}_{i+2}\cos {\mathrm{\&theta;}}_i+{\mathrm{\&eta;}}_{i+1}\cos {\mathrm{\&theta;}}_{i+1}},$ $R_{i+1,v}=\frac{{\mathrm{\&eta;}}_{i+1}\cos {\mathrm{\&theta;}}_i-{\mathrm{\&eta;}}_{i+2}\cos {\mathrm{\&theta;}}_{i+1}}{{\mathrm{\&eta;}}_{i+1}\cos {\mathrm{\&theta;}}_i+{\mathrm{\&eta;}}_{i+2}\cos {\mathrm{\&theta;}}_{i+1}}$

其中η_i+1为第i+1层的复波阻抗；

（3-3）根据电磁波极化理论获得最终的中的亮温；具体为：根据步骤

（3-2）获得的h极化的亮温TB_h和v极化的亮温TB_v，得到天线接收到的微波亮温TB，计算公式如下：

TB＝TB_v·cos²α+TB_h·sin²α

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

本发明对地形起伏粗糙面进行建模，通过定义粗糙表面“平均有效太阳辐照度”的概念，建立了考虑地形起伏（包括大尺度和小尺度）的粗糙月面微波辐射亮温正向模型，以满足微波辐射遥感定量信息反演的精度要求。具体而言，分为以下两部分：

1、通过定量计算基准面的二维倾斜情况，可以考虑大尺度起伏对亮温的影响。具体而言，不考虑大尺度地形起伏时，基准面是垂直于CE-2卫星观测方向的水平面，由于地形起伏，使得辐射计探测范围内的基准面产生倾斜，这种倾斜基准面将直接改变太阳入射角，从而影响有效太阳辐照度，引起月球表面温度的升高或者降低，最终对月球表面的亮温产生影响。小尺度起伏是通过遮蔽效应来改变太阳光的照射情况，从而影响有效太阳辐照度。

2、对于月球表面存在的小尺度地形起伏（如辐射计探测范围内月球表面相对于倾斜基准面的高度起伏），本发明利用随机粗糙面的遮蔽理论，将滩羊辐照度和遮蔽函数结合起来，定义了“有效太阳辐照度”，这样就能定量的刻画小尺度起伏对亮温的影响。

附图说明

图1是本发明起伏月面微波辐射亮温的计算方法的流程图；

图2是全局坐标系下二维随机粗糙月面示意图；

图3是微波辐射计的探测范围示意图；

图4是拟合的倾斜平面示意图；

图5是拟合平面倾斜角和方位角示意图；

图6是本地坐标系与全局坐标系中太阳入射角的对比图；

图7是本地坐标系与全局坐标系的相对关系；

图8是微波辐射计的探测范围的二维表面离散高程示意图

图9是月表平面分层结构模型示意图；

图10是A15附近地区，月球时间13.64点，CE2实测亮温和模拟亮温随纬度变化的示意图，其中图10（a）是19.35GHz对比图，10（b）是37GHz对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明的整体思路在于：本发明对地形起伏粗糙面进行建模，通过定义粗糙表面“平均有效太阳辐照度”的概念，将有地形起伏（包括大尺度和小尺度）的月壤的温度剖面计算转化为平面分层模型的温度剖面计算。具体而言，大尺度起伏是通过定量计算基准面的二维倾斜情况，即通过基准面的倾斜来改变太阳入射角和方位角，从而影响有效太阳辐照度；小尺度起伏是通过遮蔽效应来改变太阳光的照射情况，从而影响有效太阳辐照度。综合考虑以上两种尺度的起伏，对有效太阳辐照度进行修正，并将修正后的有效太阳辐照度带入热传导方程中，求解得到月壤的温度剖面，最后利用分层亮温模型以及电磁波的极化理论得到考虑地形起伏的月面微波辐射亮温。

如图1所示，本发明起伏月面微波辐射亮温的计算方法包括以下步骤：

（1）获取月球表面的微波数据参数，并利用平面拟合、坐标转换以及遮蔽函数对月球表面的微波数据参数进行计算，以获得月球表面有效太阳辐照度；本步骤包括以下子步骤：

（1-1）通过微波辐射计获取月球表面被测点的微波数据参数，包括经度lgt、纬度lat、高程H、太阳入射角θ、太阳方位角φ、对应的月球时间t、月表接收的太阳辐照度I₀、月壤中氧化铁和二氧化钛的含量S以及微波辐射计的工作频率f，这里的高程H和太阳入射角θ、太阳方位角φ都是在全局坐标系下得到的，x轴方向为经度方向，y轴方向为纬度方向，如图2所示。

（1-2）根据步骤（1-1）中获得的微波数据参数并使用最小二乘法将微波辐射计天线的月面足印内小面元中所获取的所有点拟合成一个平面，并获取该平面的倾斜角和方位角，如图3所示，其为微波辐射计未考虑地形起伏与考虑地形起伏的探测范围示意图；具体而言，是将小面元视作一个倾斜基准面和其上随机起伏的叠加。

本步骤中，曲面拟合平面的具体方法如下：

已知小面元中n个点的坐标（lgt₁,lat₁,H₁），（lgt₂,lat₂,H₂）……（lgt_n,lat_n,H_n），通过这n个点最终拟合成的平面z的表达式为

z＝b₁+b₂x+b₃y

其中，b₁、b₂和b₃为待求系数，则有

$\left(\begin{array}{c}{H}_1\\ H_2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ H_n\end{array}\right)=b_1\&CenterDot;\left(\begin{array}{c}1\\ 1\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ 1\end{array}\right)+b_2\&CenterDot;\left(\begin{array}{c}{\mathrm{lgt}}_1\\ {\mathrm{lgt}}_2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{lgt}}_n\end{array}\right)+b_3\&CenterDot;\left(\begin{array}{c}{\mathrm{lat}}_1\\ {\mathrm{lat}}_2\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{lat}}_n\end{array}\right)$

拟合平面z＝b₁+b₂x+b₃y与x轴、y轴、z轴的交点分别为

a₃(0,0,b₁)，如图4所示，

为拟合平面的外法向矢量，则有

$\left\{\begin{array}{c}{\stackrel{\&RightArrow;}{n}}_1\&CenterDot;\stackrel{\&RightArrow;}{a_1{a}_3}=0\\ {\stackrel{\&RightArrow;}{n}}_1\&CenterDot;\stackrel{\&RightArrow;}{a_2{a}_3}=0\end{array}\right.$

解得

因为取

为外法向量，其z分量应该为正，所以取

拟合平面的倾斜角α和方位角β，如图5所示。由几何关系得

$\cos \mathrm{\&alpha;}=\frac{1}{\sqrt{b_2^2+b_3^2+1}}$

$\tan \mathrm{\&beta;}=\frac{b_3}{b_2}$

至此，可求出月球表面任意区域拟合平面的倾斜角α和方位角β。

（1-3）依据步骤（1-1）和（1-2）中获得的参数，以步骤（1-2）中的拟合平面作为基准面建立本地坐标系，本地坐标系与全局坐标系中入射角的区别如图6所示，本地坐标系(X,Y,Z)与全局坐标系(x,y,z)的相对关系如图7所示，即(X,Y,Z)可由(x,y,z)经过两步得到：第一步是z轴不动，将x轴和y轴同时逆时针旋转β；第二步是y轴不动，将z轴和x轴同时向下旋转α，根据坐标转换获得太阳本地入射角θ′和太阳本地方位角φ′，其计算公式为：

cosθ′＝cosαcosθ+sinαsinθcos(φ-β)

sinθ′＝(1-cos²θ′)^1/2

sinφ′＝sinθsin(φ-β)/sinθ′

cosφ′＝(cosαsinθcos(φ-β)-sinαcosθ)/sinθ′

（1-4）依据步骤（1-1）至（1-3）中获得的参数获得微波辐射计探测范围内的月球表面斜率p、q以及月球表面的局部法向矢量表面均方根高度σ、表面相关长度l_x和l_y、表面均方根斜率w，其计算公式为：

$p=\frac{\&PartialD;H}{\&PartialD;x}$

$q=\frac{\&PartialD;H}{\&PartialD;y}$

${\stackrel{\&RightArrow;}{n}}_2=\frac{-p\hat{x}-q\hat{y}+\hat{z}}{\sqrt{1+p^2+q^2}}$

其中

分别为x轴、y轴和z轴的单位方向矢量。

对于微波辐射计探测范围内的二维表面离散高程数据，如图8所示，表面均方根高度σ为：

$\mathrm{\&sigma;}={\left[\frac{1}{M_1\&CenterDot;M_2}(\underset{m_1=1}{\overset{M_1}{\mathrm{\&Sigma;}}}\underset{m_2=1}{\overset{M_2}{\mathrm{\&Sigma;}}}{H}_{m_1{m}_2}^2-M_1\&CenterDot;M_2\&CenterDot;{\stackrel{\&OverBar;}{H}}^2)\right]}^{1/2}$

其中

$\stackrel{\&OverBar;}{H}=\frac{1}{M_1\&CenterDot;M_2}\underset{m_1=1}{\overset{M_1}{\mathrm{\&Sigma;}}}\underset{m_2=1}{\overset{M_2}{\mathrm{\&Sigma;}}}{H}_{m_1{m}_2}$

M₁，M₂分布为微波辐射计探测范围内的x和y方向的高程数目。

对于微波辐射计探测范围内的离散数据，如图8所示，在x方向，将第m₁(1≤m₁≤M₁)行的所有高程当成一维情况，则其在相距x′＝(m₃-1)Δx，m₃≥1的情况下的归一化自相关函数为

$F_1\left(x^{\&prime;}\right)=\frac{\underset{m_2=1}{\overset{M_2+1-m_3}{\mathrm{\&Sigma;}}}{H}_{m_1{m}_2}{H}_{m_1(m_2+m_3-1)}}{\underset{m_2=1}{\overset{M_2}{\mathrm{\&Sigma;}}}{H}_{m_1{m}_2}^2}$

当自相关函数F₁(x′)等于1/e时的间隔x′值被定义为表面相关长度

即

$F_1\left(l_{{\mathrm{xm}}_1}\right)=\frac{1}{e}$

最后取平均值即可得到l_x，即：

$l_x=\frac{1}{M_1}\underset{m_1=1}{\overset{M_1}{\mathrm{\&Sigma;}}}{l}_{{\mathrm{xm}}_1}$

在y方向，与x方向同理，将第m₂(1≤m₂≤M₂)列的所有高程当成一维情况，在相距y′＝(m₄-1)Δy，m₄≥1的情况下的归一化自相关函数为

$F_2\left(y^{\&prime;}\right)=\frac{\underset{m_1=1}{\overset{M_1+1-m_4}{\mathrm{\&Sigma;}}}{H}_{m_1{m}_2}{H}_{(m_1+m_4-1)m_2}}{\underset{m_1=1}{\overset{M_1}{\mathrm{\&Sigma;}}}{H}_{m_1{m}_2}^2}$

$F_2\left(l_{{\mathrm{ym}}_2}\right)=\frac{1}{e}$

$l_y=\frac{1}{M_2}\underset{m_2=1}{\overset{M_2}{\mathrm{\&Sigma;}}}{l}_{{\mathrm{ym}}_2}$

对于微波辐射计探测范围内的二维表面，可以得到x方向和y方向的均方根斜率分别为：

$w_x=\sqrt{2}\frac{\mathrm{\&sigma;}}{l_x}$

$w_y=\sqrt{2}\frac{\mathrm{\&sigma;}}{l_y}$

二维均方根斜率为：

$w=\sqrt{{(w_x\&CenterDot;{\cos \mathrm{\&phi;}}^{\&prime;})}^2+{(w_y\&CenterDot;\sin {\mathrm{\&phi;}}^{\&prime;})}^2}$

(1-5)依据步骤（1-3）和（1-4）中获得的参数并根据遮蔽函数理论获得微波辐射计探测范围内的月球表面的二维双向遮蔽函数S(p,q,θ′,φ′,θ₁)，其中观测角θ₁＝α，其计算公式根据θ₁的取值分为三种情况（这里规定太阳入射光线和观测方向位于局部法向同侧时，观测角θ₁取值为正，否则为负）：

（a）θ₁＞θ′（同侧）

$S(p,q,{\mathrm{\&theta;}}^{\&prime;},{\mathrm{\&phi;}}^{\&prime;},{\mathrm{\&theta;}}_1)=\frac{\mathrm{\&gamma;}(\stackrel{\&OverBar;}{\mathrm{\&mu;}}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})}{[\mathrm{\&Lambda;}\left(\stackrel{\&OverBar;}{\mathrm{\&mu;}}\right)+1]}$

其中

$\mathrm{\&gamma;}(\stackrel{\&OverBar;}{\mathrm{\&mu;}}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})=\left\{\begin{array}{cc}1,& \stackrel{\&OverBar;}{\mathrm{\&mu;}}\&GreaterEqual;p\sin {\mathrm{\&phi;}}^{\&prime;}+{q\cos \mathrm{\&phi;}}^{\&prime;}\\ 0,& \stackrel{\&OverBar;}{\mathrm{\&mu;}}<p\sin {\mathrm{\&phi;}}^{\&prime;}+q\cos {\mathrm{\&phi;}}^{\&prime;}\end{array}\right.,$ $\mathrm{\&Lambda;}\left(\stackrel{\&OverBar;}{\mathrm{\&mu;}}\right)=\frac{1}{2}\&CenterDot;[{\left(\frac{2}{\mathrm{\&pi;}}\right)}^{1/2}\&CenterDot;\frac{w}{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{e}^{-{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}^2/{2w}^2}-\mathrm{erfc}\left(\frac{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{\sqrt{2}w}\right)],$ erfc是误差函数的余集函数，定义为：

$\mathrm{erfc}\left(\frac{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{\sqrt{2}w}\right)=\frac{2}{\sqrt{\mathrm{\&pi;}}}{\&Integral;}_{\frac{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{\sqrt{2}w}}^{\&infin;}{e}^{-x^2}\mathrm{dx}$

（b）0＜θ₁＜θ′（同侧）

$S(p,q,{\mathrm{\&theta;}}^{\&prime;},{\mathrm{\&phi;}}^{\&prime;},{\mathrm{\&theta;}}_1)=\frac{\mathrm{\&gamma;}(\mathrm{\&mu;}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})}{[\mathrm{\&Lambda;}\left(\mathrm{\&mu;}\right)+1]}$

其中

$\mathrm{\&Lambda;}\left(\mathrm{\&mu;}\right)=\frac{1}{2}\&CenterDot;[{\left(\frac{2}{\mathrm{\&pi;}}\right)}^{1/2}\&CenterDot;\frac{w}{\mathrm{\&mu;}}{e}^{-{\mathrm{\&mu;}}^2/{2w}^2}-\mathrm{erfc}\left(\frac{\mathrm{\&mu;}}{\sqrt{2}w}\right)],$

$\mathrm{\&gamma;}(\mathrm{\&mu;}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})=\left\{\begin{array}{cc}1,& \mathrm{\&mu;}\&GreaterEqual;p\sin {\mathrm{\&phi;}}^{\&prime;}+{q\cos \mathrm{\&phi;}}^{\&prime;}\\ 0,& \mathrm{\&mu;}<p\sin {\mathrm{\&phi;}}^{\&prime;}+q\cos {\mathrm{\&phi;}}^{\&prime;}\end{array}\right.$

（c）θ₁＜0（异侧）

$S(p,q,{\mathrm{\&theta;}}^{\&prime;},{\mathrm{\&phi;}}^{\&prime;},{\mathrm{\&theta;}}_1)=\frac{\mathrm{\&gamma;}(\stackrel{\&OverBar;}{\mathrm{\&mu;}}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})\&CenterDot;\mathrm{\&gamma;}(\mathrm{\&mu;}-p\sin {\mathrm{\&phi;}}^{\&prime;}-q\cos {\mathrm{\&phi;}}^{\&prime;})}{[\mathrm{\&Lambda;}\left(\stackrel{\&OverBar;}{\mathrm{\&mu;}}\right)+1]\&CenterDot;[\mathrm{\&Lambda;}\left(\mathrm{\&mu;}\right)+1]}$

其中

μ＝cotθ′， $\mathrm{\&Lambda;}\left(\mathrm{\&mu;}\right)=\frac{1}{2}\&CenterDot;[{\left(\frac{2}{\mathrm{\&pi;}}\right)}^{1/2}\&CenterDot;\frac{w}{\mathrm{\&mu;}}{e}^{-{\mathrm{\&mu;}}^2/{2w}^2}-\mathrm{erfc}\left(\frac{\mathrm{\&mu;}}{\sqrt{2}w}\right)]$ $\mathrm{\&Lambda;}\left(\stackrel{\&OverBar;}{\mathrm{\&mu;}}\right)=\frac{1}{2}\&CenterDot;[{\left(\frac{2}{\mathrm{\&pi;}}\right)}^{1/2}\&CenterDot;\frac{w}{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{e}^{-{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}^2/{2w}^2}-\mathrm{erfc}\left(\frac{\stackrel{\&OverBar;}{\mathrm{\&mu;}}}{\sqrt{2}w}\right)].$

（1-6）依据步骤（1-1）至（1-5）获得的参数计算太阳辐照度的衰减因子IN和太阳有效辐照度I_eff，其计算公式为：

$\mathrm{IN}=\&Integral;\&Integral;\mathrm{dpdq}\frac{1}{{2\mathrm{\&pi;w}}_x{w}_y}\exp (-\frac{p^2}{{2w}_x^2}-\frac{q^2}{{2w}_y^2})\&CenterDot;S(p,q,{\mathrm{\&theta;}}^{\&prime;},{\mathrm{\&phi;}}^{\&prime;},{\mathrm{\&theta;}}_1)\&CenterDot;(p{\sin \mathrm{\&theta;}}^{\&prime;}\cos {\mathrm{\&phi;}}^{\&prime;}+q\sin {\mathrm{\&theta;}}^{\&prime;}\sin {\mathrm{\&phi;}}^{\&prime;}+\cos {\mathrm{\&theta;}}^{\&prime;})$

I_eff＝I₀·IN

其中I₀是步骤（1-1）获得的月表接收的太阳辐照度。

可以看到，本地太阳入射角、本地太阳方位角、二维粗糙度以及遮蔽函数是影响有效太阳辐照度的重要因素。

本步骤（1）的优点在于考虑了地形起伏（包括大尺度和小尺度）对有效太阳辐照度的影响。具体而言，大尺度的地形起伏是通过基准面的倾斜来改变太阳入射角和太阳方位角，从而改变有效太阳辐照度；小尺度的地形起伏（如辐射计探测范围内的点相对于基准面的高度起伏）则是通过遮蔽效应来改变太阳光的照射情况，从而改变有效太阳辐照度。

（2）根据月壤的物理参数和步骤（1）中获得的有效太阳辐照度并使用热传导理论和月壤分层模型获得月壤的温度剖面，即某一时刻月壤不同深度（每一层）的温度T；月壤的物理参数包含月壤密度ρ，月壤的介电常数ε^*、月壤的比热C，月壤的热导率K。

（2-1）将月球表面的月壤层划分成N层，如图8所示，其中第一层的厚度是2厘米，为月尘层，第二层到第N层的厚度满足以下等式：

d_m＝A·e^B·m

其中m(2≤m≤N)表示层数，e为自然对数函数的底数，A和B均为常数，其取值与月壤的厚度d成正比，并且

其中d是所求点处月壤厚度，是不可知的，但是根据实测和各种实验数据，其值一般有一个大概范围。研究表明，月球高地的月壤厚度变化范围在1.0～18.0米之间，平均厚度为12米，月海的月壤厚度在1.5～10.0米之间，平均厚度为5米。考虑到计算速度和运算复杂度，N一般取1000～2000之间的数，mindelta是月壤的最小分层间隔，一般取0.0001。

（2-2）获取月球表面物理参数的模型，月壤的物理参数包括月壤密度ρ、月壤的复介电常数ε^*（分为实部ε′和虚部ε″）、月壤的比热C和月壤的导热率K；

（2-2-1）根据登月点实测和经验公式获得月壤密度ρ（g/cm³），其为月壤深度z（cm）的双曲线模型，具体采用以下公式：

$\mathrm{\&rho;}=1.92\frac{z+12.2}{z+18}$

其中ρ是月壤密度（以克每立方厘米为单位），z是月壤深度(以厘米为单位)。月壤密度的范围为1.3g/cm³～1.92g/cm³。

（2-2-2）根据月球样品的回归分析获得复介电常数ε^*，其余与月壤密度ρ和氧化铁与二氧化钛的含量S的关系如下：

ε^*＝ε₀（ε′-jε″）

ε′＝1.919^ρ

ε″＝ε′·10^{0.038S+0.312ρ-3.260}

其中，ε₀＝0.8854F/m，ε′表示复介电常数的实部，ε″表示复介电常数的虚部，ρ表示月壤密度（g/cm³），S表示氧化铁与二氧化钛的含量（%）。

（2-2-3）根据实验获得月壤比热与月壤温度关系为：

C＝c₁T³+c₂T²+c₃T+c₄

其中c₁＝0.03142×10^-7，c₂＝-0.033662×10^-4，c₃＝0.15899×10^-2，c₄＝-0.05277，T为月壤的温度。

（2-2-4）根据Vasavada模型获得月壤热导率与月壤温度的关系为：

$K=K_c\&CenterDot;(1+\mathrm{\&chi;}\&CenterDot;{\left(\frac{T}{T_{350}}\right)}^3)$

其中K_c表示固体导热率，χ是辐射热导率与固体热导率的比值，T₃₅₀＝350K，在步骤（2-1）的分层模型中：在顶层2厘米，K_c＝9.22×10^-4Wm^-1·K^-1,χ＝1.48；在第二层到第N层,K_c＝4×10^-2Wm^-1·K^-1，χ＝1.48。

（2-3）根据步骤（2-2）获得的参数并根据热传导理论求解热传导方程，得到月壤不同深度的温度。具体为：根据半有限固体的热传导理论，月球表面物质同时受到来自外部的太阳辐射、地球反照和月球内部热流三个热源的作用。根据能量守恒定律导出的一维热传导方程为：

$\mathrm{\&rho;}(z,T)C(z,T)\frac{\&PartialD;T}{\&PartialD;t}=\frac{\&PartialD;}{\&PartialD;z}\left[K(z,T)\frac{\&PartialD;T}{\&PartialD;z}\right]+Q(z,t)$

上式指出了月壤内部中任意微元体内温度与坐标位置和时间的关系。其中，ρ(z,T)表示密度（kg/m²），C(z,T)表示比热(J/(kg·K))，K(z,T)表示热导率(W/(m·K))，Q(z,t)表示部分透明介质由于吸收太阳的辐射而产生的源项,这里不考虑这个辐射,因此Q(z,t)＝0。根据步骤（2-2）获得的参数模型，月壤密度ρ(z,T)是月壤深度z的函数，月壤比热C(z,T)和月壤热导率K(z,T)都是月壤温度T的函数。

要求解上述一维热传导方程，还需要如下两个边界条件：

（2-3-1）在月壤表层：

$K_s\frac{\&PartialD;T}{\&PartialD;z}{|}_s=\mathrm{\&epsiv;}{\mathrm{\&sigma;}}_B{T}_s^4-(1-A_b)[I_{\mathrm{eff}}+E]+J_0$

其中,

表示表面处的温度梯度，K_s为表面热导率，

表示表示传入次表面的能量；ε为红外表面发射率(一般设定为0.90-1.0之间)，σ_B为Stefan-Boltzman常数5.6703×10^-8Wm^-2·K^-4，T_s为月壤表层温度，

表示月球表面辐射的红外能量；A_b为月表热辐射反照度0.12，E表示地球反射的太阳辐照度；J₀表示月球内部发射的热通亮,大约在0.02-0.04Wm^-2之间，比月表的太阳辐照度低5个数量级，忽略不计。

（2-3-2）在热平衡深度Z₀

$\frac{\&PartialD;T}{\&PartialD;z}{|}_{\mathrm{depth}}=\frac{-J_0}{K_{\mathrm{depth}}}<<1$

其中，K_depth为在Z₀时的温度梯度，J₀是远小于1的常数，它对亮温变化没有影响，所以可以忽略。可以选择大于这个深度的任意深度作为下边界,此时 $(\&PartialD;T/\&PartialD;z){|}_{Z_0}=0.$

由以上热传导方程和边界条件，利用数值求解方法获得热传导方程的数值解，如利用前向差分法可以计算任意时刻，表面（或剖面）温度。该温度是考虑了地形起伏引起的遮蔽效应后，一个辐射计分辨率范围内的平均表面（或剖面）温度。

通过以上分析可知，月表的地形起伏影响有效太阳辐照度，而有效太阳辐照度可以代入热传导方程进行迭代求解，最终得到月壤在分层模型中每一层的物理温度。

（3）根据步骤（2）获得的温度剖面并利用burke多层平面分层亮温模型以及电磁波极化理论计算月球表面亮温。本步骤包括以下子步骤：

（3-1）根据步骤（2）中获得月壤第i（1≤i≤N）层的物理温度为T_i，复介电常数为

获得月壤第i层的复波阻抗η_i和复波数k_i，计算公式如下：

${\mathrm{\&eta;}}_i=\sqrt{\frac{{\mathrm{\&mu;}}_0}{{\mathrm{\&epsiv;}}_i^{*}}},$ $k_i=2\mathrm{\&pi;f}\sqrt{{\mathrm{\&mu;}}_0{\mathrm{\&epsiv;}}_i^{*}}$

其中μ₀＝4π×10^-7H/m。

（3-2）根据Burke推导出的具有不同温度、介电常数的多层平面分层结构辐射亮温的一阶近似解，获得辐射计在观测角θ₁接收到的亮温如下：

${\mathrm{TB}}_r=(1-{\mathrm{\&Gamma;}}_{1,r}\left({\mathrm{\&theta;}}_1\right))\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{T}_{i+1}(1-\frac{1}{L_{i+1}})(1+\frac{{\mathrm{\&Gamma;}}_{i+1,r}}{L_{i+1}})\underset{m_3=2}{\overset{i}{\mathrm{\&Pi;}}}\left(\frac{1-{\mathrm{\&Gamma;}}_{m_3,r}}{L_{m_3}}\right)$

其中，r表示极化状态，可取h或v极化，即由上式可以得到h极化的亮温TB_h和v极化的亮温TB_v，T_i+1是第i+1层的物理温度，

是辐射在第i+1层的功率损耗因子，

为功率吸收系数，为第i+1层的复介电常数，d_i+1为第i+1层的厚度，θ_i+1是传播光线在第i+1层的方向角，满足Γ_i+1是第i+1层边界上的菲涅尔反射率，Γ_i+1,r＝|R_i+1,r|²，在水平极化和铅垂极化下，分别为：

$R_{i+1,h}=\frac{{\mathrm{\&eta;}}_{i+2}\cos {\mathrm{\&theta;}}_i-{\mathrm{\&eta;}}_{i+1}\cos {\mathrm{\&theta;}}_{i+1}}{{\mathrm{\&eta;}}_{i+2}\cos {\mathrm{\&theta;}}_i+{\mathrm{\&eta;}}_{i+1}\cos {\mathrm{\&theta;}}_{i+1}},$ $R_{i+1,v}=\frac{{\mathrm{\&eta;}}_{i+1}\cos {\mathrm{\&theta;}}_i-{\mathrm{\&eta;}}_{i+2}\cos {\mathrm{\&theta;}}_{i+1}}{{\mathrm{\&eta;}}_{i+1}\cos {\mathrm{\&theta;}}_i+{\mathrm{\&eta;}}_{i+2}\cos {\mathrm{\&theta;}}_{i+1}}$

其中η_i+1为第i+1层的复波阻抗。

（3-3）根据电磁波极化理论获得最终的中的亮温。具体为：根据步骤（3-2）获得的h极化的亮温TB_h和v极化的亮温TB_v，得到天线接收到的微波亮温TB，计算公式如下：

TB＝TB_v·cos²α+TB_h·sin²α

其中α是拟合倾斜平面的倾斜角。

A15地区是月陆和月海的交界，有地形起伏较大的区域也有地势较平坦的区域。本发明选取了CE-2卫星于当地时间13.64点扫过A15地区的一条探测轨迹，在轨迹上取14个点，分别求出这14个天线足印的粗糙分层月壤微波辐射亮温，并将计算值与CE-2的多通道辐射计实测微波亮温数据进行比较，以验证本发明方法的有效性。

CE-2微波辐射天线在对月观测时，天线指向月球中心方向，即天顶角观测，所以对于可以忽略本地倾斜的被测点的天线足印区域而言，其天线观测角θ₁＝0。但对于考虑本地倾斜的区域而言，其天线观测角不为0。如图3所示，本地倾斜平面的法向偏离了天顶角方向，使得天线观测方向存在一定的观测角，这个观测角正好等于本地倾斜平面的倾斜角α。

表1列出了A15地区所选探测轨迹上14个点的相应实际参数。表中α为该点的本地倾斜基准面的倾斜角，β为本地倾斜基准面的方位角，φ为太阳全局方位角，w为二维粗糙度。

表1 A15地区13.64时刻沿轨各点的相应参数

该地区纬度较低的部分亮温变化比较明显，而纬度较高的部分亮温变化不大，因为前者属于月陆地区，地形起伏较大，对有效太阳辐照度的影响比较明显，使得温度和最终的亮温变化较大，而后者属于月海地区，地势比较平坦。从图9所示的模拟结果来看，无论是19.35GHz还是37GHz，平面模拟亮温的变化趋势与实测相差较远，所以如果忽略了月球表面地形起伏，而将月球表面视为平面，则模拟亮温会产生较大的误差。考虑了地形起伏之后，模拟亮温与实测亮温的变化趋势明显的在总体上保持一致。并且考虑地形起伏之后，点与点之间的亮温差基本与实测亮温差相近。

从仿真结果来看，本发明所建立的起伏月面微波辐射亮温正向模型比平面模型有很大的改进，其随纬度变化的趋势与实测值总体上一致，这也说明本发明模型的必要性和准确性。

通过A15地区不同频率亮温的计算及与实测亮温的对比分析可以得到以下结论：地形起伏模型亮温随纬度的变化规律与实测亮温保持一致，这与平面模型相比有很大的改善，从而说明了在计算有地形起伏亮温时必须地形起伏的影响。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种起伏月面微波辐射亮温的计算方法，其特征在于，包括以下步骤： 

（1）获取月球表面的微波数据参数，并利用平面拟合、坐标转换以及遮蔽函数对月球表面的微波数据参数进行计算，以获得月球表面有效太阳辐照度；本步骤包括以下子步骤： 

（1-1）通过微波辐射计获取月球表面被测点的微波数据参数，包括经度lgt、纬度lat、高程H、太阳入射角θ、太阳方位角φ、对应的月球时间t、月表接收的太阳辐照度I₀、月壤中氧化铁和二氧化钛的含量S以及微波辐射计的工作频率f，这里的高程H和太阳入射角θ、太阳方位角φ都是在全局坐标系下得到的； 

（1-2）根据步骤（1-1）中获得的微波数据参数并使用最小二乘法将微波辐射计天线的月面足印内小面元中所获取的所有点拟合成一个平面，并获取该平面的倾斜角和方位角； 

（1-3）依据步骤（1-1）和（1-2）中获得的参数，以步骤（1-2）中的拟合平面作为基准面建立本地坐标系，根据坐标转换获得太阳本地入射角θ′和太阳本地方位角φ′，其计算公式为： 

cosθ′＝cosαcosθ+sinαsinθcos(φ-β) 

sinθ′＝(1-cos²θ′)^1/2

sinφ′＝sinθsin(φ-β)/sinθ′ 

cosφ′＝(cosαsinθcos(φ-β)-sinαcosθ)/sinθ′ 

（1-4）依据步骤（1-1）至（1-3）中获得的参数获得微波辐射计探测范围内的月球表面斜率p、q以及月球表面的局部法向矢量n表面均方根高度σ、表面相关长度l_x和l_y、表面均方根斜率w； 

(1-5)依据步骤（1-3）和（1-4）中获得的参数并根据遮蔽函数理论获得微波辐射计探测范围内的月球表面的二维双向遮蔽函数S(p,q,θ′,φ′,θ₁)，其中观测角θ₁＝α； 

（1-6）依据步骤（1-1）至（1-5）获得的参数计算太阳辐照度的衰减因子IN和太阳有效辐照度I_eff，其计算公式为： 

I_eff＝I₀·IN 

（2）根据月壤的物理参数和步骤（1）中获得的有效太阳辐照度并使用热传导理论和月壤分层模型获得月壤不同深度的温度T； 

（2-1）将月球表面的月壤层划分成N层，第二层到第N层的厚度满足以下等式： 

d_m＝A·e^B·m

其中m表示层数，e为自然对数函数的底数，A和B均为常数，其取值与月壤的厚度d成正比； 

（2-2）获取月球表面物理参数的模型，月壤的物理参数包括月壤密度ρ、月壤的复介电常数ε^*、月壤的比热C和月壤的导热率K； 

（2-3）根据步骤（2-2）获得的参数并根据热传导理论求解热传导方程，得到月壤不同深度的温度； 

（3）根据步骤（2）获得的温度剖面并利用burke多层平面分层亮温模型以及电磁波极化理论计算月球表面亮温。 

2.根据权利要求1所述的起伏月面微波辐射亮温的计算方法，其特征在于，步骤（1-2）中，曲面拟合平面的具体方法如下： 

（1-2-1）将小面元中n个点拟合成的平面z的表达式为 

z＝b₁+b₂x+b₃y 

其中，b₁、b₂和b₃为待求系数，则有 

（1-2-2）拟合平面z＝b₁+b₂x+b₃y与x轴、y轴、z轴的交点分别为 

a₃(0,0,b₁)，

为拟合平面的外法向矢量，则有 

（1-2-3）解得

即

（1-2-4）获得拟合平面的倾斜角α和方位角β为： 

。

3.根据权利要求2所述的起伏月面微波辐射亮温的计算方法，其特征在于，步骤（1-4）的计算公式为： 

其中

分别为x轴、y轴和z轴的单位方向矢量； 

表面均方根高度σ为： 

其中 

M₁，M₂分布为微波辐射计探测范围内的x和y方向的高程数目； 

对于微波辐射计探测范围内的二维表面，方向和y方向的均方根斜率分别为： 

二维均方根斜率为： 

4.根据权利要求3所述的起伏月面微波辐射亮温的计算方法，其特征在于，步骤（1-5）的计算公式为： 

当θ₁＞θ′时， 

其中

erfc是误差函数的余集函数，定义为： 

当0＜θ₁＜θ′时： 

其中

当θ₁＜0时： 

其中

μ＝cotθ′，

5.根据权利要求1所述的起伏月面微波辐射亮温的计算方法，其特征在于，步骤（2-2）包括以下子步骤： 

（2-2-1）根据登月点实测和经验公式获得月壤密度ρ，具体采用以下公式： 

其中ρ是月壤密度，z是月壤深度； 

（2-2-2）根据月球样品的回归分析获得复介电常数ε^*，其余与月壤密度ρ和氧化铁与二氧化钛的含量S的关系如下： 

ε^*＝ε₀（ε′-jε″） 

ε′＝1.919^ρ

ε″＝ε′·10^{0.038S+0.312ρ-3.260}

其中，ε₀＝0.8854F/m，ε′表示复介电常数的实部，ε″表示复介电常数的虚部，S表示氧化铁与二氧化钛的含量； 

（2-2-3）根据实验获得月壤比热与月壤温度关系为： 

C＝c₁T³+c₂T²+c₃T+c₄

其中T为月壤的温度。 

（2-2-4）根据Vasavada模型获得月壤热导率与月壤温度的关系为： 

其中K_c表示固体导热率，χ是辐射热导率与固体热导率的比值。 

6.根据权利要求1所述的起伏月面微波辐射亮温的计算方法，其特征在于，步骤（2-3）具体为：根据能量守恒定律导出的一维热传导方程为： 

其中，ρ(z,T)表示密度（kg/m²），C(z,T)表示比热(J/(kg·K))，K(z,T)表示热导率(W/(m·K))，Q(z,t)表示部分透明介质由于吸收太阳的辐射而产生的源项，且上述方程满足以下边界条件： 

（2-3-1）在月壤表层： 

其中,

表示表面处的温度梯度，K_s为表面热导率，

表示表示传入次表面的能量；ε为红外表面发射率(一般设定为0.90-1.0之间)，σ_B为Stefan-Boltzman常数5.6703×10^-8Wm^-2·K^-4，T_s为月壤表层温度，表示月球表面辐射的红外能量；A_b为月表热辐射反照度0.12，E表示地球反射的太阳辐照度；J₀表示月球内部发射的热通亮； 

（2-3-2）在热平衡深度Z₀

其中，K_depth为在Z₀时的温度梯度，J₀是远小于1的常数。 

7.根据权利要求1所述的起伏月面微波辐射亮温的计算方法，其特征在于，步骤（3）包括以下子步骤： 

（3-1）根据步骤（2）中获得月壤第i层的物理温度为T_i，复介电常数为 

获得月壤第i层的复波阻抗η_i和复波数k_i，计算公式如下： 

其中μ₀＝4π×10^-7H/m； 

（3-2）根据Burke推导出的具有不同温度、介电常数的多层平面分层结构辐射亮温的一阶近似解，获得辐射计在观测角θ₁接收到的亮温如下： 

其中，r表示极化状态，可取h或v极化，即由上式可以得到h极化的亮温TB_h和v极化的亮温TB_v，T_i+1是第i+1层的物理温度，

是辐射在第i+1层的功率损耗因子，

为功率吸收系数，

为第i+1层的复介电常数，d_i+1为第i+1层的厚度，θ_i+1是传播光线在第i+1层的方向角，满足

Γ_i+1是第i+1层边界上的菲涅尔反射率，Γ_i+1,r＝|R_i+1,r|²，在水平极化和铅垂极化下，分别为： 

其中η_i+1为第i+1层的复波阻抗； 

（3-3）根据电磁波极化理论获得最终的中的亮温；具体为：根据步骤（3-2）获得的h极化的亮温TB_h和v极化的亮温TB_v，得到天线接收到的微波亮温TB，计算公式如下： 

TB＝TB_v·cos²α+TB_h·sin²α。 

Images