运动仿真——phased.Platform

在雷达仿真过程中，运动仿真的必要性，以及运动仿真可以实现哪些功能，在matlab对应的user guide中已经讲的很清楚了，这里不再赘述。

本文主要介绍phased.Platform的一些“坑”，和典型的用法。

第一坑：系统对象机制

系统对象（system object）在调用的时候，返回当前的状态值，并计算下一状态值存储在系统对象中，直到调用release函数复位。假如仿真的时间步长为T，第一次调用系统对象返回零时刻的状态值。如果我们想要知道nT时刻的返回值，则需要调用系统对象函数n+1次。

第二坑：InitialVelocity属性

帮助文档是这么描述该属性

红框中说：当MotionModel设置为Velocity且VelocitySource设置为Input port的时候，InitialVelocity起作用。但是问题来了，当VelocitySource设置为Input port的时候，系统对象调用的形式必须满足[Pos,Vel] = platform(T,V)，那么在调用的时候必须还要输入V参数，可以推断，第一次调用系统对象时，返回的Vel值应该是Initialvelocity，而输入参数V则是为了更新下一次数据使用。更新原则为：Pos = Pos + T*V和Vel = V。

验证代码如下：

clear;clc;v0 = [1, 2, 3]';
v1 = [2, 3, 4]';
v2 = [3, 4, 5]';rdPlatform = phased.Platform();
rdPlatform.MotionModel = 'Velocity';
rdPlatform.InitialPosition = [0;0;0];
rdPlatform.VelocitySource = 'Input port';
rdPlatform.InitialVelocity = v0;T = 1;% 输出pos=initialPosition=[0,0,0],vel=[1,2,3]，
% 更新vel = [2,3,4]
% 更新pos = initialPosition + T*v1 =[2,3,4]
[pos, vel] = rdPlatform(T, v1)% 输出pos = [2,3,4]和vel=[2,3,4]
% 更新vel = [3,4,5]
% 更新pos=[2,3,4]+[3,4,5]=[5,7,9]
[pos, vel] = rdPlatform(T, v2)% 输出pos = [5,7,9]和vel=[1,2,3]
% 更新vel = [1,2,3]
% 更新pos=[5,7,9]+[1,2,3]=[6,9,12]
[pos, vel] = rdPlatform(T, v0)

第三坑：Velocity属性

每次调用函数的时候，都不需要更新velocity，这是最基本的用法。

由于velocity这个属性tunable，所以在仿真过程中可以更新velocity的值，。

clear;clc;rdPlatform = phased.Platform();
rdPlatform.MotionModel = 'Velocity';
rdPlatform.InitialPosition = [0;0;0];
rdPlatform.VelocitySource = 'Property';
rdPlatform.Velocity = [10;20;-10];T = 1;[pos, vel] = rdPlatform(T)
rdPlatform.Velocity = [10;10;-10];
[pos, vel] = rdPlatform(T)rdPlatform.Velocity = [0;10;-10];
[pos, vel] = rdPlatform(T)

第四坑：旋转模式

旋转模式比较复杂，需要搞清楚两个问题：1）旋转什么东西？2）绕着什么旋转？

首先看旋转模式的定义：

这里说，在Circulor模式下，沿着平台方向坐标系中方位指向的顺时针方向旋转，这个说法非常难懂。因此是orientation axes绕着aizmuth direction旋转。这里出现了第一个问题：orientation axes是如何定义的？azimuth direction如何定义？

帮助文档中关于旋转模式的解释还出现在InitialOrientationAxes的定义中：

意思是orientation axes（简称oa）就是局部坐标系，这解释了第1个问题的一半疑问。紧接着第2个问题，平台有transitional和rotational的运动分量，这两种运动分量是如何影响oa的值？

根据系统对象调用的规则：[Pos,Vel,Laxes] = step(___)可以返回局部坐标系，而局部坐标系Laxes定义如下：

可以知道Laxes就是platform orientation axes，也叫platform axes，Laxes是绕着运动轨迹的方向进行旋转，这句话解释了第一个问题的后半个问题，aizmuth direction就是运动轨迹的法向量。

总结下，仿真模型中定义了平台的局部坐标系，该坐标系绕着运动轨迹的法向量旋转。这里出现了第三个问题，运动轨迹的法向量可能是随时间变化的，如何描述绕着非固定旋转轴的旋转？为了解答这个问题，我对分别只有transitional和rotational运动分量的oa变化做了研究。以下是研究过程是结论。

仅有Transitional运动分量

首先，编写代码，可视化oa和当前运动之间的关系。仿真一个平抛运动，仿真代码如下：

clear;clc;
close all;% 非零初始速度
v = [10; 10; 0];% 初始局部坐标系
initaxes = azelaxes(0, 0);% 仿真时间步长
T = 0.1% 选择加速运动模式
rdPlatform = phased.Platform(...
'MotionModel','Acceleration', ...
'InitialPosition',[0,0,0]',...
'InitialVelocity',v, ...
'AccelerationSource','Property', ...
'Acceleration',[0,0,-9.8]');% 选择非旋转模式
% 初始的局部坐标系设置为initaxes；
rdPlatform.ScanMode = 'None';
rdPlatform.InitialOrientationAxes = initaxes;
rdPlatform.OrientationAxesOutputPort = true;% 创建绘图窗口
figure;
ha = axes("Parent", gcf);
hold(ha, 'on');
xlabel('x')
ylabel('y')
zlabel('z')% 循环仿真
for i = 1:50 [pos, vel, lax] = rdPlatform(T);xaxes(vel)*transpose(xaxes(v))*initaxeslaxpltax(ha, pos, vel, lax);
endaxis(ha, 'equal')

仿真结果如下：

可以看出，oa随着速度矢量的变化而变化。假如每一时刻的速度矢量为$v_i$，由速度矢量确定的坐标系（我自己写一个函数，来确定平台局部坐标系，该函数根据$x$轴来确定局部坐标系，因此叫xaxes)为$P_i$，那么$P_i$可以由$P_1$通过线性变换$T_i$得到，写作

$$P_i=T_i\times P_1$$

由此可以得到线性变换为$T_i=P_i\times P_1^T$，假设平台没有旋转运动，那么平台初始局部坐标系为$L_1$，第$i$时刻的局部坐标系$L_i$就是$L_i=T_i\times L_1=P_i\times P_1^T\times L_1$。

为了验证我的结论，可以修改上述代码的初始速度、加速度矢量为任意值，比较每一次仿真输出的局部坐标系的值和我们自己计算的局部坐标系的值是否相等。经过我的测试验证，满足结论。

仅有Rotational运动分量

假设旋转坐标系的初始为$R_1$，旋转角度为$\theta$，那么旋转变换可以表示为$\mathrm{rotz}(\theta )$，那么加旋转的局部坐标系就是$\mathrm{rotz}(\theta )\times {\mathrm{R}}_1$

clear;clc;close all;v = [0; 0; 0];scanrate = 10;rdPlatform = phased.Platform(...
'MotionModel','Acceleration', ...
'InitialPosition',[0,0,0]',...
'InitialVelocity',v, ...
'AccelerationSource','Property', ...
'Acceleration',[0, 0, 0]');rdPlatform.ScanMode = 'Circular';
rdPlatform.InitialScanAngle = 0;
rdPlatform.AzimuthScanRate = scanrate;
rdPlatform.InitialOrientationAxes = azelaxes(0, -20);
rdPlatform.OrientationAxesOutputPort = true;T = 1;refpos = [0; 0; 0];figure;ha = axes("Parent", gcf);
hold(ha, 'on');
xlabel('x')
ylabel('y')
zlabel('z')v0 = xaxes(v);R1 = azelaxes(0, -20);for i = 1:5[pos, vel, lax] = rdPlatform(T);% [~, ang] = rangeangle(lax(:,1))% R = rotz(scanrate*(i-1)*T);rotz(-scanrate*(i-1)*T) * R1laxpltax(ha, pos, vel, lax);
endaxis(ha, 'equal')

含有Transitional和Rotational运动分量

该模式下的建模机理待讨论补充

附件

function lax = xaxes(x)x = x / norm(x);
% define vector v in x-y plane of global coordinate system
v = [0; 0; 1];
% compute y;
y = cross(v, x);
y = y / norm(y);
% compute z;
z = cross(x, y);
% construct local axes;
lax = [x, y, z];
end

function pltax(ha, pos, vel, ax)vel = vel / norm(vel);
scatter3(ha, pos(1), pos(2), pos(3), 'k');
quiver3(ha, pos(1), pos(2), pos(3), vel(1), vel(2), vel(3), 'k', 'MaxHeadSize', 2, 'AutoScale', 'off');
quiver3(ha, pos(1), pos(2), pos(3), ax(1, 1), ax(2, 1), ax(3, 1), 'b', 'MaxHeadSize', 2, 'AutoScale', 'off');
quiver3(ha, pos(1), pos(2), pos(3), ax(1, 2), ax(2, 2), ax(3, 2), 'r', 'MaxHeadSize', 2, 'AutoScale', 'off');
quiver3(ha, pos(1), pos(2), pos(3), ax(1, 3), ax(2, 3), ax(3, 3), 'r', 'MaxHeadSize', 2, 'AutoScale', 'off');end