一种基于最大吸收功率的卫星太阳能电池板安装方案

问题的引入

卫星上的电源有太阳能电池电源、化学电源和核电源等多种。太阳能是一种取之不尽、用之不竭的能源。因此，目前绝大多数卫星都采用太阳能电池。这种电池利用半导体材料的光电效应将太阳能直接转换成电能，能工作几年甚至几十年。
早期的太阳能电池片是用半导体硅做成的，现在开始广泛采用砷化镓材料，因为它的光电转换效率高，能达到 20%以上。另外，为了进一步增加太阳能电池片数量，中大型卫星多采用由数块太阳能电池板连接而成的太阳翼。由于太大，火箭装不下，所以在发射时太阳翼处于折叠状态，待卫星与火箭分离后再展开。为了使卫星上的太阳翼总是朝向太阳，以获取最大的电能，不少卫星上的太阳翼采用了一些先进技术：一是装有驱动机构，带着太阳翼转动；二是用太阳敏感器来捕获太阳的方位，从而控制驱动机构转动太阳翼，使太阳光尽可能垂直于太阳翼，为卫星提供充足的能源。
假设有一颗立方体卫星，我们可以在它的各个面安装太阳能电池板。卫星环绕地球运动，但轨道任意。怎样安装太阳能电池板，才能使卫星吸收最大的太阳辐射功率？本文通过Matlab平台仿真，对这个问题进行了讨论。

分析思路

首先获取一颗卫星的两行轨道根数（TLE）数据，提取开普勒轨道根数，以及卫星在近焦点坐标系中的坐标。在此基础上，我们将以下列步骤推进这个问题：

1. 利用普朗克黑体辐射公式，计算可见光波段（400nm~800nm）太阳的辐射强度总和M、辐射通量∅、太阳的发光强度I；
2. 利用已知参数，得到太阳在近焦点坐标系中的坐标；
3. 将太阳在近焦点坐标系中的坐标转换为地心惯性坐标系中的坐标；
4. 再将太阳在地心惯性坐标系中的坐标转换到卫星的近焦点坐标系中的坐标；
5. 利用太阳与卫星的几何关系和已求出的发光强度I，得到太阳在卫星处的照度；
6. 通过已求出的照度，结合太阳光照方向与卫星各面法向的几何关系，计算卫星各面太阳能电池板实际最大输出功率；
7. 分析各方向电池板最大输出功率数据，提出能够优化太阳能电池功率 输出的立方体卫星太阳能电池板安装或工作方案。

太阳光照强度参数的计算

由普朗克黑体辐射公式，黑体辐射的单色辐射强度为
$$M(\lambda ,T_0)=10^{12}\frac{c_1}{{\lambda }^5(e^{c_2/\lambda T_0}-1)}$$
式中λ为波长(um)，$T_0$为黑体的温度(K)，通常太阳辐射可以认为是温度为 5900K 的黑体辐射，$c_1=3.742\times 10^{(}-4)W\bullet u{m}^2$为第一辐射常数，$c_2=14388um\bullet K$为第二辐射常数。单色辐射强度$M(\lambda ,T_0)$单位为$W/(m^2\bullet um)$，在${\lambda }_1-{\lambda }_2$波长范围内太阳的发光强度总和$M$与辐射通量$\varnothing$分别为
$$M={\int }_{{\lambda }_1}^{{\lambda }_2}M(\lambda ,T_0)\mathrm{d}\lambda ={\int }_{{\lambda }_1}^{{\lambda }_2}10^{12}\frac{c_1}{{\lambda }^5(e^{c_2/\lambda T_0}-1)}\mathrm{d}\lambda$$
$$\varnothing =4\pi {R_S}^{2}M$$
式中$R_s=6.9599\times 10^{8}m$为太阳半径，$M$的单位为$W/(m^2)$。假设太阳所发出的总辐射通量在空间方向上的分布是均匀的，则在${\lambda }_1-{\lambda }_2$波长范围内太阳的发光强度为
$$I=\varnothing /4\pi =R_s^{2}M$$
根据距离平方反比定律，在${\lambda }_1-{\lambda }_2$波长范围内太阳在卫星处的照度(即大气层外的照度)$E_s$为
$$E_s=I/D^2=R_s^{2}M/D^2$$
根据上式，若知道太阳与卫星的距离，便可得到太阳在卫星处的照度。但根据实际情况，卫星可能运行到地球阴影面，导致光照消失。在后续计算照度时应注意这一点。

太阳在各坐标系中坐标的转换

考虑到计算照度时需要得到太阳到卫星的距离，将太阳坐标转换到卫星的近焦点坐标系将会便于求解。
太阳轨道参数：

已知以上参数后，便能够求解太阳在其自身的近焦点坐标系中的坐标。太阳在近焦点坐标系中的坐标为：
$$x_{sun,o}=a(cosE-e)$$
$$y_{sun,o}=a\surd (1-e^2)sinE$$
$$z_{sun,o}=0$$
近焦点坐标系到地心赤道惯性坐标系的转移矩阵可表示为：

因此太阳在地心赤道惯性坐标系中的坐标为：

用同样的方法，我们可以得到卫星的近焦点坐标到地心赤道惯性坐标系的转移矩阵${T_{ip}}^{\ast }$。设$Y_{isun}$为太阳在卫星近焦点坐标系中的坐标，$Y_{p}sun$为太阳在地心赤道惯性坐标系中的坐标，则
$$Y_{p}sun={T_{ip}}^{\ast }\ast Y_{i}sun$$
利用矩阵左除的方法，可由$Y_{p}sun$和${T_{ip}}^{\ast }$得到$Y_{i}sun$，从而达成太阳坐标转换的目标。

太阳在卫星处照度的计算

在照度的计算中，最关键的是确定当前时刻卫星所在位置是否能够收到光线，即卫星是否处于地球阴影面。
任意时刻，太阳、地球、卫星都处于同一个平面，它们的相对位置关系如下图：

图中，蓝色的圆表示地球，红色的圆表示太阳。卫星在围绕地球的轨道上运动，图中黑色的点A标出的是卫星有照度-无照度的临界位置：若卫星到太阳的距离大于该位置到太阳的距离，则卫星进入地球阴影区，无光照。
由于太阳到地球的距离远远大于地球半径，可近似认为太阳光是平行光，这样，点A的位置可由几何关系确定。设D是过A向地球圆所作切线与地球圆的焦点，则$AD\Vert OC$。
延长OC，过A向OC作垂线，垂足为B，则AB=OD=地球半径；设卫星在近焦点坐标系中的位置坐标为$(x,y,z)$，则$AO=\surd (x^2+y^2+z^2)$；由于ABO是直角三角形，可以根据勾股定理计算出BO的长度；OC表示太阳到地心的距离，可由太阳的坐标直接求出，得到BC=BO+OC；同时ABC也是直角三角形，AB、BC长度已知，根据勾股定理便能够计算出AC的长度，而AC就表示了临界点卫星到太阳的距离$D_0$。
在卫星运行的每一时刻，都对卫星到太阳的距离D进行计算，并将计算结果与临界点的距离值进行比较。若$D>D_0$，说明卫星处在阴影区，照度为0；若$D<D_0$，说明卫星处在光照区，依照公式$E_s=I/D^2=R_s^{2}M/D^2$计算照度即可。

卫星各面太阳能电池板实际最大输出功率的计算

立方体卫星的各面示意如下：

由于卫星在轨道坐标系统中俯仰、偏 航、翻滚角均为0度，故卫星六个面中，必有一面法线平行于轨道切向（A），一面法线指向地心（B）。与这两面方向不同的另一面（C），其法线可由这两面法线的叉乘得到。剩下的三个面与这三个面相对，法向刚好相反。也就是说，有了A面和B面的法线方向，便能求出全部6个面的方向。
A面的法向，即轨道切向，可由两相邻时刻卫星的运动距离矢量得到。当时间差选取地足够小，卫星在两相邻时刻运动的距离矢量，近似平行于卫星在该位置处的切线；由于计算在卫星的近焦点坐标系中进行，B面的法向就是卫星的位矢，由卫星的坐标直接确定。
太阳光线方向可由太阳与卫星的坐标差确定。设太阳光线矢量a,卫星某一面法向为b，a与b的夹角为θ,则

如果一块太阳能电池板面积为S，电能转化效率为η，最大输出功率为$W_{max}$，与太阳的距离为D，那么太阳能电池板实际最大输出功率为
$$W=max(0,min(\eta E_{s}cos\theta ,W_{max}))$$

Matlab代码讲解

1）计算太阳光照强度与位置参数
依照光照强度参数计算及坐标转换中的数据及公式，代入计算即可。
太阳位置参数计算：

asun=149597870;                                 %半长轴
esun=0.016709-0.4204*0.1^3;                     %偏心率
isun=0;                                         %轨道面倾角
Msun=pi*(357.5277233+35999.050)/180;            %真近地点角
omegasun=pi*(282.9404+1.7201)/180;              %近地点幅角
Esun=Msun+esun*(1+esun*cos(Msun))*sin(Msun);    %偏近点角
ibusasun=pi*(23.439291-0.0130042)/180;          %黄赤夹角

400-800nm太阳发光强度计算：

syms lamda
T0=5900;c1=3.742*10^(-4);c2=14388;
f=10^12*c1/(lamda^5*exp(c2/(lamda*T0))-1);
M=double(int(f,lamda,0.4,0.8))*10^(-6);
Rs=6.9599*10^8;
I=Rs^2*M;

需注意对符号函数进行定积分运算时，最好对计算结果进行double强制类型转换，否则可能因类型不匹配而报错。

对太阳进行坐标转换
首先根据公式计算太阳在它自己的近焦点坐标系中的坐标：

xsunp=asun*(cos(Esun)-esun);
ysunp=asun*sqrt(1-esun^2)*sin(Esun);
zsunp=0;

然后利用转移矩阵，转换到赤道惯性坐标系：

T1=[1 0 0;0 cos(ibusasun) -sin(ibusasun);0 sin(ibusasun) cos(ibusasun)];
T2=[cos(omegasun) -sin(omegasun) 0;sin(omegasun) cos(omegasun) 0;0 0 1];
Tsun=T1*T2;
[Ysuni]=Tsun*[xsunp;ysunp;zsunp];

自行定义了函数ECI2PF，用以完成从地心惯性坐标系到卫星的近焦点坐标系的坐标转换：

function [Y1] = ECI2PF (Y0, Kep)

% 开普勒参数
a = Kep(1);  Omega = Kep(4);
e = Kep(2);  omega = Kep(5);
i = Kep(3);  M0    = Kep(6);

%转换矩阵计算
Tip1=[cos(Omega),-sin(Omega),0;sin(Omega),cos(Omega),0;0,0,1];
Tip2=[1,0,0;0,cos(i),-sin(i);0,sin(i),cos(i)];
Tip3=[cos(omega),-sin(omega),0;sin(omega),cos(omega),0;0,0,1];
Tip=Tip1*Tip2*Tip3;

% 坐标系转换
Y1=Tip\Y0;
end

函数的输入为赤道惯性坐标系中的坐标Y0和卫星的开普勒参数Kep，输出所求点在卫星的近焦点坐标系中的坐标Y1。Tip原本是近焦点坐标系-赤道惯性坐标系的转移矩阵，这里利用MATLAB带有的“左除”这一特殊运算，完成了反向转换的功能。
通过该函数，就把太阳的坐标从赤道惯性坐标系转到了卫星的近焦点坐标系：

[Ysun]=ECI2PF(Ysuni,Kep);

照度的计算
太阳坐标已有，现在确定卫星的坐标：

[Y0] = State (gm, Kep, dt);
Ys=0.1^3*Y0(:,1);

State函数在上次实验中已有使用，其功能是确定卫星在它自己的近焦点坐标系中的位置和速度坐标，在这里不再详细阐述。第二行命令的作用是提取卫星的位置坐标，并将其存放在向量Ys中。
根据日-地-卫星图中说明的原理，卫星是否处于地球阴影区取决于它与太阳的距离：

ds_sun0=6371^2+(sqrt(sum(Ysun.^2))+sqrt(sum(Ys.^2)-6371^2))^2;
ds_sun=sum((Ys-Ysun).^2);

ds_sun0表述临界距离，ds_sun表示卫星与太阳的实际距离。当实际距离大于临界距离时，照度为0.否则按照照度计算公式计算卫星处的照度：

if ds_sun>ds_sun0
    Es(dt/60)=0;
else
    Es(dt/60)=I/ds_sun;
end

对上述步骤按一定时间步长进行迭代（控制总时长12000s），即可得到从规定时刻起200分钟内太阳在卫星处的照度曲线：

for dt=60:60:12000
…
end

计算卫星各面段太阳能板实际最大输出功率
这一步关键在于找到3个互相垂直的法线方向与太阳光照方向（太阳与卫星连线方向）的夹角余弦值。
首先是卫星前进方向，法线方向与轨道切线平行，利用很短时间内卫星的位移方向近似替代。
在前一时刻，把卫星坐标存放在临时变量temp中：

temp=Ys;

到下一时刻，用现在的卫星坐标与temp相减，就是短时间内卫星的位移：

Yds=Ys-temp;

利用矢量夹角公式，得到该方向的夹角余弦值：

costheta1=sum(Ysun.*Yds)/sqrt(sum(Ysun.^2)*sum(Yds.^2));

其次是法线指向地心的方向。由于各坐标均属于卫星的近焦点坐标系，卫星的坐标就表示了该法线方向，运算相对方便：

costheta2=sum(Ysun.*Ys)/sqrt(sum(Ysun.^2)*sum(Ys.^2));

最后一个法向用前两个法向的叉乘表示：

Yon=cross(Ys,Yds);
costheta3=sum(Ysun.*Yon)/sqrt(sum(Ysun.^2)*sum(Yon.^2));

利用三个余弦值，结合太阳能电池板实际最大输出功率的计算公式，便能得到各方向太阳能板实际最大输出功率（有效面积60平方厘米，效率30%，最大输出功率2.4W）：

W1(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*costheta1,2.4));
W2(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*(-costheta1),2.4));
W3(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*costheta2,2.4));
W4(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*(-costheta2),2.4));
W5(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*costheta3,2.4));
W6(dt/60)=max(0,min(0.3*60*0.1^4*Es(dt/60)*(-costheta3),2.4));

仿真结果分析

1）确定时刻起至之后 200 分钟内太阳在卫星处的照度曲线
卫星：amateur；起始时刻：2020.10.4.22:10

曲线形如矩形方波。由于卫星轨道的大小范围远小于日地距离，故卫星在非阴影区运动时，其到太阳的距离变化对照度的影响可以忽略不计，近似为恒定值；而一旦运动到阴影区，照度立马减小到0。这说明仿真结果是合理的，与实际情况能够较好契合。

2）立方体卫星各面太阳能板实际最大输出功率曲线
卫星各面如下图所示：

以下为各面最大输出功率在2020.10.4.22:10到此后200分钟的曲线图及各图的分析：
正面与背面

上图中有两个完全相同的脉冲波形，两波形之间输出功率恒为零。由于卫星的绕转具有周期性，可估算amateur卫星的绕转周期约7000s。
波形在靠近2000s、超出8000s处有两个陡增点，且这两个陡增点与照度曲线的陡增点完全匹配，说明卫星在这两个时刻从阴影区进入了光照区，并受到一定功率的太阳光照射；而后曲线缓慢上升，又逐渐减小到0，这也是与实际情况相符合的：

位置1：卫星进入光照区，照度突变为一定值；
位置2：θ=0，功率达到最大；位置1-位置2 θ减小，功率增大
位置3：θ=90，功率减小至0；位置2-位置3 θ增大，功率减小；位置3-位置1过程由于正面始终背向太阳光，功率保持0不变。

背面与正面的输出功率曲线具有对称性：位置3（图7中大于4000s处的起始点）是正面输出功率减小到0的点，却是背面输出功率从0开始逐渐增加的点；从位置3开始，背面输出功率逐渐增加到最高点（位置4），又逐渐下降，指导进入地球阴影区（位置5），功率突变减小至0。
右面与左面

背向地球的右面，最大输出功率曲线始终是平滑的，说明这一面未受到阴影区的影响：从位置2-位置3，右面法线与阳光方向夹角减小，功率增大；位置3-位置4，右面法线与阳光方向夹角增大，功率减小；而在受到阴影区影响的位置2-位置4过程，由于几何遮挡，功率始终为0。

假设没有地球的存在，左面的功率曲线应与右面在一个周期内实现对称。但因为地球的遮挡，左面本应到达最大功率的位置5-位置1区间，完全没有照度，在图像上体现为在原本的脉冲峰值处呈现矩形凹陷，且凹陷点的时刻与照度曲线相匹配。
上面与下面


对amateur卫星来说，上面最大输出功率始终为零，下面最大输出功率则同照度曲线完全一致。若更换不同的卫星，则上面、下面的功率曲线可互换，但呈现的规律保持不变。也就是说，对于这两个面，在非阴影区，必有一个面始终受到光照且其法线方向同阳光方向的夹角保持不变，而另一个面始终处于被几何遮挡的状态。这与实际情况也是相符合的。

总结–立方体卫星太阳能电池板安装方案的设计

对于立方体卫星太阳能电池板的安装方案，有以下假设：
① 为每个面都安装电池板是不现实的；
②电池板的安装应满足对称性，即某一面若安装电池板，其对应面也应安装电池板；
② 最优安装方案应使得所有电池板获得的总功率在可选方案中最大。
为选定最优安装方案，我们需要求出各面的最大输出功率，或者是最大输出功率的相对大小，并对结果进行排序。这一步的程序实现为：

S1=sum(W1);
S2=sum(W2);
S3=sum(W3);
S4=sum(W4);
S5=sum(W5);
S6=sum(W6);
[S,I]=sort([S1+S2 S3+S4 S5+S6],'descend')

W1-W6表示分别正、背、右、左、上、下面在各时刻的最大输出功率，则S1-S6就能够代表各面平均最大输出功率的相对。由于假设②，进行排序的是正+背、右+左、上+下面平均功率之和。一组面的平均功率之和越大，在其上安装太阳能电池板的效率越高。
使用不同卫星得到的结果如下：

卫星名称	正+背面平均功率相对大小	右+左面平均功率相对大小	上+下面平均功率相对大小	排序（降序）
amateur	119.2251	103.6662	190.2461	上+下＞正+背＞右+左
nnss	154.4545	115.8150	103.6619	正+背＞右+左＞上+下
noaa	132.4896	112.3542	143.6318	上+下＞正+背＞右+左
radar	91.2918	85.8307	290.1408	上+下＞正+背＞右+左
geo	68.3337	307.3655	49.9478	右+左＞正+背＞上+下
gps-ops	131.7696	270.6468	72.5789	右+左＞正+背＞上+下

选取多组卫星并进行排序比较后，可以得到这样的结果：

① 正+背面，即法线平行于轨道切线的两个面，其平均功率总不会最小；
②右+左与上+下面，其平均功率与具体的卫星有很大的关系：一般来说，高纬度运行的卫星，其上+下面获得的功率最大；赤道平面卫星或低纬度运行的卫星，其右+左面获得的功率最大。
据此，一个立方体卫星若有两组方向可以安装太阳能电池板，其最佳安装方案为：
①正+背面必安装电池板；
②若卫星运行轨道纬度较高，在上+下面安装电池板；若为赤道平面卫星或低纬度运行的卫星，在右+左面安装电池板。