Matlab连接STK获取可见性、距离数据

1. 常用命令
   用Matlab创建STK对象
   创建场景：

uiap = actxserver(‘STK11.application’);
root = uiap.Personality2;
root.NewScenario(‘myscenario1’);
sc = root.CurrentScenario;

创建卫星：

sat = sc.Children.New(‘eSatellite’,‘mysat’);
sat.Propagator.Propagate;#显示卫星轨迹

STK创建好对象后，用Matlab获取
获取stk中所有对象的名称（包括卫星、地面站、星座等对象）：

stkInit; 
remMachine = stkDefaultHost;
conid = stkOpen(remMachine);
objNames = stkObjNames;

获取可见性数据：

% 参数依次是：参考卫星，目标卫星，报告样式，开始时间，结束时间，时间间隔
[secData, secNames] =stkAccReport(char(objNames(int32(i))),char(objNames(int32(j))), style,startTime,endTime,dt);

secData是元胞数据(相当于python中的字典dict)，secData{1}获取第一个单元的数据，secData{1}.data返回可见次数、开始时间集、结束时间集、每次持续时间。

Matlab得到可见性数据及可见性矩阵（01矩阵）

#每隔10sec获取一次星间可见性
stkInit;
remMachine = stkDefaultHost;
conid = stkOpen(remMachine);

objNames = stkObjNames;
dt = 10;
style = 'Access';
startTime = 0;
endTime = 10;
result = zeros(24);

for k = 1:60
    for i = 4:35
        disp(strcat(num2str(k),' ---------- ',num2str(i-3)))
        result(i-3,i-3) = 1;
        for j = (i+1):35
            [secData, secNames] =stkAccReport(char(objNames(int32(i))),char(objNames(int32(j))), style,startTime,endTime,dt);
            [temp1,temp2] = size(secData{1});
            if temp2 == 4
                [a,b,c,d] = secData{1}.data;
                % a: Access b:Start Time c:End Time d:Duration
                if d == dt
                    result(i-3,j-3) = 1;
                    result(j-3,i-3) = 1;
                end
            end
        end
    end
    csvwrite(strcat('C:\\Users\\Administrator\\Desktop\\stk\\stk_files\\result_',num2str(k),'.csv'),result);
    startTime = startTime + 10;
    endTime = endTime+10;
end
stkClose(conid);
stkClose;

（可跳过此代码！）
使用for+eval方便计算多卫星或多地面站之间的可见性：

'''
使用此方法连接stk需用matlab代码创建完整场景，包括卫星等，然后操纵stk对象
使用上面conid的方法连接stk，可以在stk创建好场景，然后matlab获得对象名称
'''
    uiApplication = actxGetRunningServer('STK11.application');
    root = uiApplication.Personality2;
    scenario = root.CurrentScenario;
    
    %waitbar
    hWait = waitbar(0,'Please Wait~');
    
    %scenario start time
    startTime = scenario.StartTime;
    formatIn = 'dd mmm yyyy HH:MM:SS';
    startTimeNum = datenum(startTime,formatIn);
    
    
    % allChildren includes satellites and facilities
    allChildren = scenario.Children;
    
    allSatellites = allChildren.GetElements('eSatellite');
    allFacilities = allChildren.GetElements('eFacility');
    
    
    satNum = allSatellites.Count;
    for i=0:1:satNum-1
        % eg: sat0=allSatellites.Item(cast(0,'int32'));
        eval(['sat',num2str(i),'=allSatellites.Item(cast(i,''int32''));'])
    end
    
    facNum = allFacilities.Count;
    for i=0:1:facNum-1
        eval(['fac',num2str(i),'=allFacilities.Item(cast(i,''int32''));'])
    end
    
    % all access in Time(day month year...)
    allAccessIntervals=[];
    % all access in Seconds
    allAccessIntervalsSec=[];
    
    % only ISL
    if(var1==1)
        
        for from=0:1:satNum-1
            waitbar(from/satNum,hWait);
            for to=from+1:1:satNum-1
                eval(['fromSat=sat',num2str(from),';']);
                eval(['toSat=sat',num2str(to),';']);
                
                access = fromSat.GetAccessToObject(toSat);
                access.ComputeAccess;
                
                accessIntervals = access.ComputedAccessIntervalTimes;
                if(accessIntervals.Count~=0)
                    
                    computedIntervals = accessIntervals.ToArray(0, -1);
                    
                    temp=cell(accessIntervals.Count,2);
                    for i=1:1:accessIntervals.Count
                        temp{i,1}=from+1;
                        temp{i,2}=to+1;
                    end
                    allAccessIntervals = [allAccessIntervals;temp,computedIntervals];
                    
                end
                
            end
        end
    end
 
    % transfer Time(day month year time) to Seconds
    [rows,cols] = size( allAccessIntervals);
    for row=1:1:rows
        intervalStart =  allAccessIntervals{row,3};
        intervalEnd =  allAccessIntervals{row,4};
        
        %in day
        intervalStartNum = datenum(intervalStart,formatIn);
        intervalEndNum = datenum(intervalEnd,formatIn);
        
        %in sceconds
        intervalSatrtSec = (intervalStartNum-startTimeNum)*24*60*60;
        intervalEndSec = (intervalEndNum-startTimeNum)*24*60*60;
        
        allAccessIntervalsSec =    [allAccessIntervalsSec;allAccessIntervals{row,1},allAccessIntervals{row,2},intervalSatrtSec,intervalEndSec];
    end

Matlab连接STK后，获取星间可见性数据及距离数据完整代码：
dt为时间间隔，由于卫星的移动性特点，在此设dt为60sec获取每隔60sec的星间可见性矩阵及距离矩阵：
主文件：

[conid,objNames] = initConn();%运行一次即可，运行后跑其他代码，将这行注释！
bias = 3;%偏移量
satnum = 32;%卫星总数
startobj = 1 + bias;
endobj = satnum + bias;
startTime = 0;
endTime = 7200;
dt = 60;
result = getVisibility(objNames,startTime,endTime,satnum,bias);
analysis(result,startTime,endTime,satnum,dt);
result_range = getRange(objNames,satnum,bias);
analysisRange(result_range,startTime,endTime,satnum,dt);
closeConn(conid);

函数文件：
getVisibility.m

% 获取从startTime到endTime之间所有的可见性数据
function [result] = getVisibility(objNames,startTime,endTime,satnum,bias)
    style = 'Access';
    dt = 1;
    result = {};
    for i = 1:satnum
        disp(strcat(num2str(i)))
        result{i,i} = 1; % 自己对自己设置为可见 1
        for j = (i+1):satnum
            [secData, ~] =stkAccReport(char(objNames(int32(i+bias))),char(objNames(int32(j+bias))),style,startTime,endTime,dt);
            result{i,j} = secData;
            result{j,i} = secData;
        end
    end
end

analysis.m

% 解析可见性数据
function [] = analysis(result,startTime,endTime,satnum,dt)
    for k = startTime:dt:endTime%每隔dt(sec)输出可见性矩阵
        disp(strcat(num2str(k)))
        visibility = zeros(satnum);%可见性矩阵01矩阵
        for i = 1:satnum
            visibility(i,i) = 1;%自己和自己都是可见的
            for j = (i+1):satnum
                if iscell(result{i,j})  % 自己和其他卫星%判断是否为元胞数据结构
                    [~,tmp2] = size(result{i,j}{1,1});
                    if tmp2 == 4 % 可能存在可见性
                        [access,startTime,stopTime,durationTime] = result{i,j}{1,1}.data;
                        [accessLength1,~] = size(access);
                        % 检查是否在可见时间段内
                        for z = 1:accessLength1
                            if k >= startTime(z) && k <=stopTime(z) % 存在可见性
                                visibility(i,j) = 1;
                                visibility(j,i) = 1;
                                break;
                            end
                        end
                    end
                end
            end
        end
        csvwrite(strcat('C:\\Users\\Administrator\\Desktop\\stk\\stk_files\\visibility_',num2str(k),'.csv'),visibility)
    end
end

getRange.m

% 获取距离
function [result] = getRange(objNames,satnum,bias)
    style = 'AER';
    result = {};
    for i = 1:satnum
        disp(strcat(num2str(i)))
        result{i,i} = 0; % 自己与自己的距离为 0
        for j = (i+1):satnum
            [rangeData, ~] =stkAccReport(char(objNames(int32(i+bias))),char(objNames(int32(j+bias))),style);
            result{i,j} = rangeData;
            result{j,i} = rangeData;
        end
    end
end

analysisRange.m

% 解析数据
function [] = analysisRange(result,startTime,endTime,satnum,dt)
    for k = startTime:dt:endTime
        visibility = load(strcat('C:\\Users\\Administrator\\Desktop\\stk\\stk_files\\visibility_',num2str(k),'.csv'));
        disp(strcat(num2str(k)))
        range = zeros(satnum);
        for i = 1:satnum
            range(i,i) = 0;% 对自身距离为 0
            for j = (i+1):satnum
                if iscell(result{i,j})% 其他卫星
                    [~,tmp2] = size(result{i,j}{1,1});
                    if visibility(i,j) == 1 % 存在距离
                        [time,azimuth,elevation,rangeTmp] = result{i,j}{1,1}.data;
                        [rangeLength,~] = size(rangeTmp);
                        % 检查是否在可见时间段内
                        for z = 1:rangeLength
                            if z <= rangeLength - 1
                                if k >= time(z) && k <time(z+1)
                                    range(i,j) = rangeTmp(z);
                                    range(j,i) = rangeTmp(z);
                                    break;
                                end
                            else
                                range(i,j) = rangeTmp(z);
                                range(j,i) = rangeTmp(z);
                            end
                        end
                    else% 不存在距离  
                        range(i,j) = 999;
                        range(j,i) = 999;
                    end
                end
            end
        end
        csvwrite(strcat('C:\\Users\\Administrator\\Desktop\\stk\\stk_files\\range_',num2str(k),'.csv'),range)
    end
end