在雷達仿真過程中,運動仿真的必要性,以及運動仿真可以實現哪些功能,在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