基于Matlab模拟VLBI 基线天空覆盖

✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab仿真内容点击👇

智能优化算法  神经网络预测 雷达通信  无线传感器

信号处理 图像处理 路径规划 元胞自动机 无人机  电力系统

⛄ 内容介绍

基于甚长基线干涉测量技术(Very Long Base Interferometer,VLBI)的时延,时延率以及USB(United S-band)/UXB(United X-band)的测距特点,采用联合统计定位及月面高程约束策略,实现了高精度的探测器的实时单点定位和准实时联合统计定位,且实时单点定位不受力学约束能够快速准确地给出三维位置信息

⛄ 部分代码

%

% function vlbiCoverageMap([{stations}, 'sourcename'])

%

% Plot 2D maps of baseline-elevation against (RA x Dec) for every baseline 

% formed by {stations}, and adds the track of source 'sourcename'. 

% Allowed station names are those returned by [~,names]=get_vlbi_stations_xyz()

% A NED database lookup is performed to get the source RA,Dec coodinates. 

% The "baseline-elevation" is the minimum of the source elevations at 

% either antenna on a given baseline.

%

% Requires:  Aerospace Toolbox and Mapping Toolbox

%

% Example plots:

%   vlbiCoverageMap({'KT','KU','KY'}, 'Sgr A*')

% The plots are exported to files that are

%   vlbi-cover-<plotNr>-<station1>-<station2>.pdf

%

function [figfiles] = vlbiSkyCoverage(varargin)

    min_elev = 10.0;    % minimum elevation (degrees)

    freq = 86e9;        % observing frequency (Hz) used to estimate baseline length

    %% Helpers

    fnt = struct();

    fnt.FontSize = 10; 

    fnt.FontWeight = 'bold'; 

    fnt.FontName = 'Times New Roman';

    figfiles = {};

    

    %% Args

    if ~(length(varargin) == 0 || (length(varargin) == 2))

        fprintf(1, 'vlbiCoverageMap() : specify either 2 or zero arguments!\n');

        return

    end

    % Defaults

    if (length(varargin) == 0),

        stations = {'KU', 'KT', 'Eb', 'Pv', 'CARMA', 'SMA', 'APEX', 'JCMT'};

        target = 'NRAO530';

    else

        stations = varargin{1};

        target = varargin{2};

    end

    

    %% Geometry

    [xyz,~] = get_vlbi_stations_xyz(stations);

    lla = ecef2lla(xyz, 'WGS84');

    fprintf(1, 'Looking up %s on NED...\n',target);

    [~,src] = m_getNED(target);

    %% Make grid of sky RA & Dec positions

    Nst = size(lla,1);

    if (Nst <= 1),

        fprintf(1, 'Error: need 2 or more stations, got %d (%s)\n', Nst, char(stations));

        return;

    end

    

    Npt = 100;

    decs = deg2rad(linspace(-90,+90,Npt)); %decs(end+1) = decs(1);

    ras  = deg2rad(linspace(0,360,Npt)); %ras(end+1) = ras(1);

    lats = deg2rad(lla(:,1));

    lons = deg2rad(lla(:,2));

    [decsMN,rasMN] = meshgrid(decs,ras);

    %% Evaluate Az,El across sky RA&Dec grid at site of every antenna/baseline

    baselines = nchoosek(1:Nst, 2);

    Nbl = size(baselines,1);

    for bb=1:Nbl,

        s1 = baselines(bb,1);

        s2 = baselines(bb,2);

        % station 1: alt, az 

        HA  = rasMN + lons(s1);

        lat = lats(s1);

        [alt1,az1] = m_pos2altaz(lat,HA,decsMN);

        % station 2: alt, az 

        HA  = rasMN + lons(s2);

        lat = lats(s2);

        [alt2,az2] = m_pos2altaz(lat,HA,decsMN);

        % baseline: alt clip

        if 0,

            altbl = (alt1 + alt2) / 2;

        else

            altbl = min(alt1,alt2); 

        end

        % altbl(altbl < min_elev) = NaN;

        % sky map of baseline

        figure(1), close(gcf);

        figure(1), clf;

        axesm('MapProjection','hammer', 'AngleUnits','degrees', 'grid','on', ...

            'frame','on', 'meridianlabel','on', 'labelformat','signed', ...

            'parallellabel','on', 'mlabelparallel',0);

        geoshow(rad2deg(decsMN), rad2deg(-rasMN), altbl, 'DisplayType', 'contour');

        hgeo = gca();

        box off; hold on;

        % overlay the station names

        x = rad2deg(lats); 

        y = rad2deg(lons);

        plotm(x(s1), y(s1), 'bo', 'MarkerFaceColor', 'r');

        plotm(x(s2), y(s2), 'bo', 'MarkerFaceColor', 'r');

        h1 = textm(x(s1), y(s1)+180/50, stations{s1});

        h2 = textm(x(s2), y(s2)+180/50, stations{s2});

        set(h1,fnt);

        set(h2,fnt);

        % overlay source track

        x = ones(size(ras)) * src.dec;

        y = rad2deg(ras);

        x(end+1) = x(1);

        y(end+1) = y(1);

        plotm(x,y, 'k-.');

        plotm(src.dec, src.RA, 'bo', 'MarkerFaceColor', 'g');

        h = textm(src.dec+5*sign(src.dec), src.RA, target);

        set(h,fnt);

        % add color scale of elevation

        axh = colorbar('Peer',gca());

        set(get(axh,'YLabel'),'String','Elevation on baseline (deg)','FontSize',12);

        % add title, add axis labels

        dlonH = (12/pi)*abs(lons(s1)-lons(s2));

        dlatD = (180/pi)*abs(lats(s1)-lats(s2));

        lambda = 299792458.0/freq;

        dlambda = norm(xyz(s1,:)-xyz(s2,:), 2) / lambda;

        [maxEl,maxElIdx] = max(altbl(:)); 

        [M1,M2] = ind2sub(size(altbl),maxElIdx);        

        s={};

        s{end+1} = ['\makebox[4in][c]{', sprintf('%s -- %s', stations{s1}, stations{s2}), '}'];

        s{end+1} = ['\makebox[4in][c]{', sprintf('dHA=%.2f\\,h dDec=%.0f\\,deg B=%.0f\\,M$\\lambda$ @ %.0f\\,GHz', dlonH,dlatD,dlambda*1e-6,freq*1e-9), '}'];

        s{end+1} = ['\makebox[4in][c]{', sprintf('Maximum El=%.0f\\,deg for Dec=%.0f\\,deg', maxEl,rad2deg(decsMN(M1,M2))), '}'];

    

        title(s, 'Interpreter','latex','FontSize',12,'HorizontalAlignment','center');

        xlabel('Celestial RA (deg)');

        ylabel('Celestial Dec (deg)');

        

        % fix the declination color scale

        set(hgeo,'clim',[-90 90]);

        

        % export to graphics file (PNG: low quality, SVG: not supported in Matlab R2014)

        figfiles{end+1} = sprintf('vlbi-cover-%02d-%s-%s', numel(figfiles)+1, stations{s1}, stations{s2});

        saveas(gcf(),[figfiles{end} '.pdf'],'pdf');

        W = 29.7; H = 21.0; Z = 3.0;

        set(gcf,'PaperPositionMode','auto', 'PaperOrientation','portrait', 'PaperUnits','centimeters', 'Units','centimeters');

        set(gcf,'PaperSize',[W H], 'PaperPosition',[-Z, -Z, W+Z, H+Z], 'InvertHardcopy','off', 'Resize','on');

        print(gcf(),[figfiles{end} '.png'],'-dpng','-r300');

        if (bb ~= Nbl), input('next'); end

    end

end

function [alt1,az1] = m_pos2altaz(lat,HA,dec)

    % http://www.stargazing.net/kepler/altaz.html

    % note: input args are in radians, and can be matrices

    alt = asin(sin(dec).*sin(lat) + cos(dec).*cos(lat).*cos(HA));

    az  = (sin(dec) - sin(alt).*sin(lat)) ./ (cos(alt).*cos(lat));

    

    % remove roundoff errors

    [i,j] = find(abs(az)>1.0); 

    az(i,j) = sign(az(i,j)); 

    % change hemisphere

    az = acos(az);

    az(sin(HA)>0) = deg2rad(360) - az(sin(HA)>0); 

    

    alt1 = rad2deg(alt);

    az1 = rad2deg(az);

end

function [entry,sentry] = m_getNED(srcname)

    % Query that returns ASCII:

    % http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=MRK+1496&extend=no&out_csys=Equatorial&out_equinox=J2000.0&obj_sort=RA+or+Longitude&of=ascii_bar&img_stamp=NO    

    nedURL = 'http://ned.ipac.caltech.edu/cgi-bin/objsearch';

    nedopts = { 'objname', srcname', ...

        'extend','no', ...

        'out_csys','Equatorial', ...

        'out_equinox','J2000.0', ...

        'obj_sort','RA or Longitude', ...

        'of','ascii_bar', ...

        'img_stamp','NO' };

    [str,status] = urlread(nedURL, 'get',nedopts);

% fprintf(1, 'NED returned: %d - %s\n', status, str);

    str = strtrim(str);

    % Note:

    % When a source is found in NED, status==1 and NED's reply looks like:

    % No.|Object Name|RA(deg)|DEC(deg)|Type|Velocity|Redshift|Redshift Flag|Magnitude and Filter|Distance (arcmin)|References|Notes|Photometry Points|Positions|Redshift Points|Diameter Points|Associations

    % 1|MESSIER 081|148.88822 |  69.06529 |G|   -34|-0.000113 |    | 7.89||1879|32|186|9|18|8|1

    entry = {'', 0, 0, 0, 'xxxx-xx', {}, {}};

    sentry = struct();

    sentry.name = '';

    sentry.RA = 0.0;

    sentry.dec = 0.0;

    

    if (status==1),

        % TODO: the NED reply *can* contain several lines. For example a

        % look-up of NGC 4649 returns 536 source names. How to treat this...

        obj = regexp(str,'\n','split');

        obj = obj{end};

        data = regexp(char(obj),'\|','split');

        if (numel(data)<7),

            fprintf(1, 'surfDb_getNED: unexpected reply from NED:\n');

            data,

            return;

        end

        name = data{2};

        name = upper(name);

        name = strrep(name,'GROUP','');        

        vsys = floor(str2double(data{6}));

        sra  = str2double(data{3});

        sdec = str2double(data{4});

        sz   = str2double(data{7});

        

        entry = {name, vsys, sz, 0.0, 'xxxx-xx', {}, {}, sra, sdec};

        sentry.name = name;

        sentry.RA = sra;

        sentry.dec = sdec;

    else

        entry = {'', 0, 0.0, 0.0, 'xxxx-xx', {}, {}, 0.0, 0.0};

    end

end

⛄ 运行结果

⛄ 参考文献

[1]王霄, 雷辉, 杨旭海,等. 基于VLBI观测的基线改正方法[J]. 测绘科学, 2019, 44(10):7.

⛄ 完整代码

❤️部分理论引用网络文献，若有侵权联系博主删除

❤️ 关注我领取海量matlab电子书和数学建模资料