본문 바로가기

작은 공부방

Project - 위성의 궤도 plot 프로그램

반응형

Introduction

  ,TLE 파일은 위성의 상대적인 위치, 궤도 정보 등 다양한 데이터를 가지고 있다. TLE 파일이 있으면 특정 시간의 위성 위치를 알 수 있다. 따라서 이러한 일련의 데이터를 바탕으로 위성을 Visualization하는 프로그램을 제작하고자 한다. 본 프로그램은 Pyvista를 이용해서 제작되었다. 본 프로그램은 독립적인 프로그램으로 gRPC를 이용해서 구현하였다. 따라서 다른 외부 프로그램에서 이를 쉽게 이용할 수 있다.

 

  본 프로그램이 수집하는 데이터는 다음과 같다.

  • 현재 시간
  • 위성 정보(Serial Number로 제공)
  • 위성의 현재 자세 정보(quaternion으로 제공)
  • 기준 벡터 정보

  본 프로그램이 나타내는 정보는 다음과 같다.

  • 지구 표현
  • 위성의 궤도 표시
  • 위성의 위치 표시
  • 지상국 위치 표시
  • 지상국과 통신 가능 위치 표시

  본 프로그램이 요구하는 데이터를 넣으면 3차원으로 자유롭게 상호작용할 수 있는 인터페이스가 나타난다. 따라서 본 프로그램이 어떻게 제작되었는지 아래에 설명하고자 한다.

 

Implementation

기본 가정

  결국 지구와 위성을 표현하는데 가장 중요한 것은 두 물체의 상대적인 위치를 표현하는 것이다. 이때, ECI나 ECEF와 같이 다양한 좌표계가 존재한다. 하지만 어떤 방식을 쓰던지 문제는 없는데, 결국 상대적인 위치 표현이기 때문이다.(별이나 절대적인 대상으로 그리는게 아니다.). 본 구현에서는 skyfiled.api의 wgs84에서 제공하는 satellite 객체를 이용했다. 이 객체는 위성의 latitude, longitude, height 정보를 제공하는데, 이를 이용하면 위도와 경도가 0인 지점에서의 상대적인 위치를 알 수 있기 때문이다. Frame은 다음과 같다.

X axis Y axis Z axis
attitude = 0, longitude = 0인 점과 지구 중심 사이를 연결한 선 X와 Z의 외적 지구의 자전축

 

입력파트

  위 프로그램이 동작하는데 필요한 모든 인자들을 넘겨주지 않으면 프로그램이 동작하지 않는다. 해당 프로그램에 인자를 넘기는 method는 gRPC를 이용해서 구현하였다. 특히 현재 시간, 위성 정보, 위성의 자세 정보 총 3가지의 필수 인자가 입력되지 않으면 동작하지 않는다. 이를 위해서 gRPC를 이용했는데, gRPC 문법은 기초 문법을 그대로 이용했으며, gRPC를 위해서 만든 통신 규약만 설명하고자 한다.

syntax = "proto3";

service Satrac{
  rpc SendTrajectoryInfo (SatTrajectoryInfo) returns (Response_code) {}
  rpc SendAttitudeInfo (SatAttitudeInfo) returns (Response_code) {}
  rpc SendGroundStationInfo (GroundStationInfo) returns (Response_code) {}
  rpc SendSatVecInfo (stream SatVectorInfo) returns (Response_code) {}
  rpc ChangeUrl (Url) returns (Response_code) {}
  rpc StartDraw (Empty) returns (Response_code) {}
}

  서비스 이름은 Satrac으로 정의하였으며, 궤도 정보 전송, 자세 정보 전송, 지상국 정보 전송, 위성 기타 벡터 정보 전송은 위성에 관련된 정보이며, TLE 파일을 불러오는 주소를 바꿀 수 있도록 하였다. 각각은 종료 후 Response_code를 반환하며, SendSatVecInfo에서는 여러 개의 백터를 한 번에 전송할 수 있으므로 stream 형태로 구현하였다.

message SatTrajectoryInfo {
  SatDate t = 1;
  int32 sat_num = 2;
}

message SatDate {
  int32 y = 1;
  int32 mo = 2;
  int32 d = 3;
  int32 h = 4;
  int32 m = 5;
  int32 s = 6;
}

 위성을 세팅하는데 TLE를 불러올 수 있는 정보와 현재 위치를 알 수 있는 date 정보를 포함하고 있다. 나머지는 업로드한 파일을 참고하면 된다.

 

출력파트

지구 표현

  지구를 표현하는 것은 생각보다 간단하다. pyvista는 pyplot과 비슷하게 시각화 도구이다. 대신 pyplot이 3d와 관련해서 조금 부족한 부분이 있다면 pyvista는 이를 간편한 방법으로 제공한다. 먼저 지구를 구형으로 가정하고 반지름이 지구와 같은 구를 만들었다. 이후에 지구를 펼처놓은 지도를 이 구에 mapping하는 식으로 지구를 표현하였다. 지구를 표현하는 부분은 다음과 같다.

earth_radius = 6378.1

sphere = pv.Sphere(radius=earth_radius, theta_resolution=240, phi_resolution=240, start_theta=270.0001, end_theta=270)
sphere.active_t_coords = np.zeros((sphere.points.shape[0], 2))


for i in range(sphere.points.shape[0]):
    x, y, z = sphere.points[i, 0]/earth_radius, sphere.points[i, 1]/earth_radius, sphere.points[i, 2]/earth_radius
    x, y, z = normalize(x), normalize(y), normalize(z)
     
     sphere.active_t_coords[i] = [0.5 + math.atan2(-x, y)/(2 * math.pi), 0.5 + math.asin(z)/math.pi]
  #cordinate chage
 
sphere.rotate_z(-90)

#load satellite

pl = pv.Plotter()
pl.add_background_image('starry-night-sky-fit.jpg', scale=1.001)
pl.add_mesh(sphere, texture=pv.read_texture("earth.jpg"), smooth_shading=False)

  코드의 핵심적인 부분을 이야기 하자면 2D인 지도를 3D로 매핑하는 기술이 필요한데 이는 UV mapping이라고 한다. 그렇게 복잡하지는 않으며, 변환식을 이용해서 계산한 2D 결과를 sphere 객체의 active_t_coords에 넣어주는 것으로 구현할 수 있다. 이후에 위의 기본 가정에서의 Frame과 지구를 맞추기 위해서 -90도를 돌려서 지구를 표현해 주었다.

 

위성의 궤도 표시 & 위성의 위치 표시

  위성의 궤도는 TLE 정보를 이용하면 계산할 수 있다. 이러한 method를 제공하는 것은 skyfiled.api의 wgs84이다. load를 이용해서 불러오면 List 형태로 satellite 객체를 만든다. 이후 원하는 위성을 선택하면 satellite 객체를 이용해서 utc 시간에서의 궤도의 위치를 찾을 수 있다. 이후 이렇게 받은 위치에 대한 latitude, longitude, height 자료를 이용해서 (height,0,0) 벡터를 y-axis에 대해서 longitude만큼 회전하고, lattitude만큼  z축을 회전하여서 좌표를 나타낸다. 위성의 위치는 위성의 궤도 중에서 현재 시간의 궤도 정보를 갖고 있으면 바로 알 수 있다. 

  위성의 궤도와 더불어서 위성이 해당 지점에 있을 때, 태양에 보이는지 보이지 않는지를 표현하기 sunlit 메소드를 이용하였다. eph파일을 불러오고 이를 바탕으로 태양에 보이는지 보이지 않는지도 표시하였다.

def lat_lon_rotation(lat, lon, vector):
    s = math.sin(lat * math.pi / 180)
    c = math.cos(lat * math.pi / 180)
    dcm_y = np.array([[c, 0, s],
                      [0, 1, 0],
                      [-s, 0, c]])

    s = math.sin(lon * math.pi / 180)
    c = math.cos(lon * math.pi / 180)
    dcm_z = np.array([[c, -s, 0],
                      [s, c, 0],
                      [0, 0, 1]])

    res = (dcm_z @ (dcm_y @ vector))

    return res.T

def lat_lon_at(satellite, ti, earth_radius):
    geo = satellite.at(ti)

    lat, lon = wgs84.latlon_of(geo)
    lat, lon = lat.arcminutes()/60, lon.arcminutes()/60
    height = wgs84.height_of(geo).km

    return lat_lon_rotation(lat, lon, np.array([earth_radius + height, 0, 0]).T)
    
for ti in self.time_gap:
    x, y, z = lat_lon_at(self.satellite, ti, self.earth_radius)

    sunlit = self.satellite.at(ti).is_sunlit(self.eph)

 

위성의 자세 표시 & 지상국 위치 표현

  위성의 자세는 일반적으로 quaternion으로 표현되므로, 해당 위성 자세정보가 들어오게 되면 이 quaternion을 dcm으로 바꾸어서 벡터를 바로 변환해줄 수 있도록 한다. 이 dcm은 위성체 stl파일과 위성의 x,y,z 및 기준 벡터를 표현하는데 이용된다. 따라서 이 dcm으로 적절한 벡터 연산을 이용해서 위성의 자세롤 표현한다.

  지상국도 dcm을 이용한다는 점은 똑같은데 마찬가지로 지상국 stl파일이 초기에 가지고 있는 위치 정보를 이용해서 적당한 회전변환을 해줌으로써 위치를 올바르게 만들고, 표현하였다.

for i in range(3):
    pl.add_lines(np.array([sat_location + sat_att[i],sat_location]), color = color[i], width = 1.5)

for i in range(len(sat_vec)):
    pl.add_lines(np.array([sat_location , sat_location + 500*(sat_att_dcm@((sat_vec[i]).T)).T]), color = color[i+3], width = 1.5)

 

지상국과 통신 가능 위치 표현

  지상국의 위도와 경도를 제공하면 지상국이 볼 때 위성의 위치가 위도 몇 도인지 알려주는 find_event method가 있다. 따라서 이 method를 이용하게 된다면. 어느 시점에 위성이 떠오르고, 지는지를 알 수 있다. 실제 통신 가능 범위도 이와 유사한 방식으로 계산되기 때문에 이를 이용해서 표현할 수 있다.

lat, lon = math.asin(normalize(self.gs_location[2]/self.earth_radius)), math.atan2(self.gs_location[1], self.gs_location[0])
        
crisp = wgs84.latlon(lat, lon)
t_deg, events_deg = self.satellite.find_events(crisp, self.now_time, self.end_time, altitude_degrees=self.view_degree)

two_time_set = []

taos = self.now_time
aos_flag = False
for ti, event in zip(t_deg, events_deg):
    if event == 0:
        taos = ti

    elif event == 2:
        two_time_set.append([taos, ti])

if aos_flag == True:
    two_time_set.append([taos, self.end_time])

 

Result

  이러한 여러가지 것을 고려하여서 프로그램 제작을 완료하였다. 완료한 프로그램은 다음과 같다. gRPC로 제작하여서 서버를 실행 시킨 후 client도 실행 시켜야 프로그램이 동작하는 것을 볼 수 있다. 실제 외부 프로그램에서 gRPC를 이용하면 이를 따로 열 수도 있다. 실행 화면은 다음과 같다.

simulation.zip
5.52MB

 

 

반응형