在雷达仿真过程中,运动仿真的必要性,以及运动仿真可以实现哪些功能,在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 v_i vi,由速度矢量确定的坐标系(我自己写一个函数,来确定平台局部坐标系,该函数根据 x x x轴来确定局部坐标系,因此叫xaxes)为 P i P_i Pi,那么 P i P_i Pi可以由 P 1 P_1 P1通过线性变换 T i T_i Ti得到,写作
P i = T i × P 1 P_i=T_i\times P_1 Pi=Ti×P1
由此可以得到线性变换为 T i = P i × P 1 T T_i=P_i\times P_1^T Ti=Pi×P1T,假设平台没有旋转运动,那么平台初始局部坐标系为 L 1 L_1 L1,第 i i i时刻的局部坐标系 L i L_i Li就是 L i = T i × L 1 = P i × P 1 T × L 1 L_i=T_i\times L_1=P_i\times P_1^T\times L_1 Li=Ti×L1=Pi×P1T×L1。
为了验证我的结论,可以修改上述代码的初始速度、加速度矢量为任意值,比较每一次仿真输出的局部坐标系的值和我们自己计算的局部坐标系的值是否相等。经过我的测试验证,满足结论。
仅有Rotational运动分量
在这个模式下,就很容易理解The current platform axes rotate around the normal vector to the path of the platform.
。具体含义为:平台有一个和全局坐标系重合的局部坐标系,在局部坐标系中构建旋转坐标系,记为 R 1 R_1 R1。沿着方位角方向旋转角度 θ \theta θ,就是沿着z轴看向局部坐标系的原点,绕着局部坐标系的z轴,顺时针旋转。那么旋转变换可以表示为 r o t z ( − θ ) \rm{rotz}(-\theta) rotz(−θ),那么加旋转的局部坐标系就是 r o t z ( − θ ) × R 1 \rm{rotz}(-\theta)\times R_1 rotz(−θ)×R1。
另外还有一个参数是InitialScanAngle
,该参数决定了初始角度。由于我们已经定义了和全局坐标系重合局部坐标系,那么考虑初始角度 r o t z ( − θ ) × r o t z ( i n i t A n g ) × R 1 \rm{rotz}(-\theta)\times \rm{rotz}(initAng)\times R_1 rotz(−θ)×rotz(initAng)×R1。
用代码实现如下:
clear;clc;close all;% 初始速度为零
v = [0; 0; 0];T = 1;% 初始orientation axes
initOA = azelaxes(0, 10);% 初始旋转角度
initAng = 90;% 扫描速度为10deg/sec
scanrate = 45;% 无加速度
rdPlatform = phased.Platform(...
'MotionModel','Acceleration', ...
'InitialPosition',[0,0,0]',...
'InitialVelocity',v, ...
'AccelerationSource','Property', ...
'Acceleration',[0, 0, 0]');% 圆周扫描,初始扫描角度为45deg。
rdPlatform.ScanMode = 'Circular';
rdPlatform.InitialScanAngle = initAng;
rdPlatform.AzimuthScanRate = scanrate;
rdPlatform.InitialOrientationAxes = initOA;
rdPlatform.OrientationAxesOutputPort = true;figure;
ha = axes("Parent", gcf);
hold(ha, 'on');
xlabel('x')
ylabel('y')
zlabel('z')for i = 1:3[pos, vel, lax] = rdPlatform(T);rotz(-scanrate*(i-1)*T) * rotz(initAng) * initOAlaxpltax(ha, pos, vel, lax);
endaxis(ha, 'equal')
仿真结果如下,可以看出仿真结果和我的猜想一致。
含有Transitional和Rotational运动分量
包含平动和转动的运动状态,该状态下返回的Laxes值目前机理还不清楚,待讨论。
附件
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