Luquiescene

站点统计

文章

57

分类

5

标签

14

总字数

58,715

运行天数

0 天

最后活动

0 天前

公告

欢迎来到Luquiesceneの小屋ヾ(≧▽≦*)o

了解更多

标签

Luquiescene

公告

欢迎来到Luquiesceneの小屋ヾ(≧▽≦*)o

了解更多

标签

AI coding Arduino C++ Git GitHub Python 博客天文学学习笔记工具教程数据科学日常硬件

Luquiescene

公告

欢迎来到Luquiesceneの小屋ヾ(≧▽≦*)o

了解更多

标签

AI coding Arduino C++ Git GitHub Python 博客天文学学习笔记工具教程数据科学日常硬件

站点统计

文章

57

分类

5

标签

14

总字数

58,715

运行天数

0 天

最后活动

0 天前

2176 字

6 分钟

笔记｜天文数据科学02｜周期侦探

2026-07-08

理工研习

Python

/

数据科学

/

天文学

在上一篇文章中，我们通过肉眼检查了光变曲线中的凌星信号。但真正的数据分析需要量化方法——用数学来比较不同模型的优劣，并用计算机算法自动搜索周期性信号。

本文将逐步实现这些方法，最终用真实的 NASA Kepler 太空望远镜数据验证我们的技术。

1. 模型比较：哪个模型更好？#

1.1 绘制观测数据#

首先创建一条简单的光变曲线——5个亮度的观测值随时间变化：

1
import numpy as np
2
import matplotlib.pyplot as plt
3

4
time = np.array([1, 2, 3, 4, 5])
5
brightness = np.array([1.02, 0.98, 1.01, 0.99, 1.00])
6

7
plt.scatter(time, brightness)
8
plt.plot(time, brightness)
9

10
plt.xlabel('Time measurements')
11
plt.ylabel('Brightness')

plt.scatter() 绘制散点（实际观测值），plt.plot() 用折线连接它们以展示变化趋势。亮度在 1.0 附近上下波动。

1.2 最简单的模型：恒定亮度#

如果我们假设恒星亮度始终不变，那么可以构建一个”恒定模型”——将所有时间的预测值设为 1.00：

1
model_brightness = 1.00
2
predicted = np.array([model_brightness] * len(brightness))
3

4
print(predicted)
5
plt.scatter(time, brightness)
6
plt.plot(time, predicted)

输出：

1
[1. 1. 1. 1. 1.]

[model_brightness] * len(brightness) 将标量值重复 5 次，创建一个长度匹配的数组。

1.3 计算残差（误差）#

残差 = 预测值 − 观测值。正残差表示模型高估了亮度，负残差表示模型低估：

1
error = predicted - brightness
2

3
print('Observed brightness', brightness)
4
print('Predicted brightness', predicted)
5
print('Errors', error)

输出：

1
Observed brightness [1.02 0.98 1.01 0.99 1.  ]
2
Predicted brightness [1. 1. 1. 1. 1.]
3
Errors [-0.02  0.02 -0.01  0.01  0.  ]

1.4 残差总和的问题#

将残差直接相加看似合理，但正负误差会相互抵消：

1
total_error = np.sum(error)
2

3
print('Errors', error)
4
print('Sum of errors', total_error)

输出：

1
Errors [-0.02  0.02 -0.01  0.01  0.  ]
2
Sum of errors 0.0

残差总和为 0，但模型显然不完美！因为正的 +0.02 和负的 −0.02 相互抵消了。因此不能直接使用残差总和来衡量模型好坏。

1.5 残差平方和（SSE）#

解决方案：将误差平方后再求和。平方后所有值都为正，不会相互抵消：

$\text{SSE} = \sum (\text{predicted} - \text{observed})^2$

1
squared_error = (predicted - brightness) ** 2
2
print('Errors:', error)
3
print('Squared errors:', squared_error)
4

5
sum_squared_error = np.sum(squared_error)
6
print('Sum of squared error:', sum_squared_error)

输出：

1
Errors: [-0.02  0.02 -0.01  0.01  0.  ]
2
Squared errors: [0.0004 0.0004 0.0001 0.0001 0.    ]
3
Sum of squared error: 0.0010000000000000018

SSE 越小，模型越贴合数据。平方操作放大了大误差的影响，惩罚了偏离程度较大的点。

1.6 比较两个模型#

现在比较两个恒定亮度模型：model_1 = 1.00、model_2 = 0.97：

1
model_1 = 1.00
2
model_2 = 0.97
3

4
predicted_1 = np.array([model_1] * len(brightness))
5
predicted_2 = np.array([model_2] * len(brightness))
6
errors_1 = predicted_1 - brightness
7
errors_2 = predicted_2 - brightness
8

9
squared_error_1 = errors_1 ** 2
10
squared_error_2 = errors_2 ** 2
11

12
total_error_1 = np.sum(squared_error_1)
13
total_error_2 = np.sum(squared_error_2)
14

15
print('Errors_1:', total_error_1)
16
print('Errors_2:', total_error_2)
17

18
# model 1 is better

输出：

1
Errors_1: 0.0010000000000000018
2
Errors_2: 0.00550000000000001

SSE₁ (0.001) < SSE₂ (0.0055)，模型 1 更好。将两个模型可视化：

1
plt.scatter(time, brightness)
2
plt.plot(time, predicted_1)
3
plt.plot(time, predicted_2)
4

5
plt.xlabel("Time")
6
plt.ylabel("Brightness")
7
plt.title("Which model is better?")

橙色线（model_2 = 0.97）整体低于数据点，SSE 更大。蓝色线（model_1 = 1.00）更接近数据中心。

1.7 自动搜索最佳模型#

手动测试两个模型效率太低。可以用 np.linspace 在可能范围内生成 101 个候选值，遍历计算每个的 SSE，最后用 np.argmin 找到 SSE 最小的模型：

1
possible_models = np.linspace(0.95, 1.05, 101)
2

3
sse_values = []
4

5
for model in possible_models:
6
    predicted = np.array([model] * len(brightness))
7
    errors = brightness - predicted
8
    sse = np.sum(errors ** 2)
9
    sse_values.append(sse)
10

11
sse_values = np.array(sse_values)
12

13
best_index = np.argmin(sse_values)
14
best_model = possible_models[best_index]
15
best_sse = sse_values[best_index]
16

17
print("Best model brightness:", best_model)
18
print("Smallest squared error:", best_sse)

输出：

1
Best model brightness: 1.0
2
Smallest squared error: 0.0010000000000000018

最佳模型亮度 = 1.00，与我们之前的直观判断完全一致。这种”暴力搜索”方法虽然简单，但非常有效——它的核心思想与后续 BLS 算法搜索最佳周期一脉相承。

2. 比较含凌星与不含凌星的模型#

2.1 生成含凌星的合成数据#

创建一个 100 个点的光变曲线，第 4.2 至 5.0 天存在凌星（亮度降为 0.97），并叠加测量噪声：

1
time = np.linspace(0, 10, 100)
2
true_brightness = np.ones_like(time)
3

4
in_transit_true = (time > 4.2) & (time < 5.0)
5
true_brightness[in_transit_true] = 0.97
6

7
rng = np.random.default_rng(42)
8
noise = rng.normal(0, 0.005, size=len(time))
9
brightness = true_brightness + noise
10

11
plt.scatter(time, brightness)
12
plt.show()

np.ones_like(time) 创建与 time 形状相同的全 1 数组。(time > 4.2) & (time < 5.0) 用布尔掩码标记凌星区间。

2.2 比较两个模型#

“无凌星模型”假设亮度始终为 1.0；“有凌星模型”在第 4.2−5.0 天将亮度设为 0.97：

1
model_no_transit = np.ones_like(time)
2
model_transit = np.ones_like(time)
3

4
in_transit_true = (time > 4.2) & (time < 5.0)
5
model_transit[in_transit_true] = 0.97
6

7
sse_no_transit = np.sum((brightness - model_no_transit) ** 2)
8
sse_transit = np.sum((brightness - model_transit) ** 2)
9
print("SSE no transit model: ", sse_no_transit)
10
print("SSE transit model: ", sse_transit)
11

12
plt.scatter(time, brightness)
13
plt.plot(time, model_no_transit)
14
plt.plot(time, model_transit)

含凌星模型的 SSE 更小——它更好地解释了数据。这说明用模型的优劣可以量化判断是否存在凌星信号。

3. 真实 Kepler 数据分析#

3.1 安装工具包#

从这一天开始，我们将使用真实的 NASA Kepler 太空望远镜数据。安装必要的 Python 包：

1
!pip -q install lightkurve astroquery astropy

lightkurve：专门用于分析 Kepler/K2/TESS 光变曲线的 Python 库
astroquery：查询 NASA 天体物理数据库
astropy：天文学核心工具库

1
import numpy as np
2
import pandas as pd
3
import matplotlib.pyplot as plt
4

5
import requests
6
from io import StringIO
7

8
import lightkurve as lk
9
from astropy.timeseries import BoxLeastSquares
10
from scipy.ndimage import median_filter
11

12
print("Tools are ready.")

输出：

1
Tools are ready.

3.2 选择团队与搜索范围#

每个团队被分配到不同的周期搜索范围，模拟真实科研中的分工协作：

1
TEAM = 1  # change to 1, 2, or 3
2

3
TEAM_INFO = {
4
    1: {"emoji": "🔥", "name": "Brown dwarfs", "period_range": (0.5, 3.0)},
5
    2: {"emoji": "🌡", "name": "We came here for free wi-fi ", "period_range": (3.0, 10.0)},
6
    3: {"emoji": "🪐", "name": " Trinity ", "period_range": (10.0, 30.0)},
7
}
8

9
info = TEAM_INFO[TEAM]
10
p_lo, p_hi = info["period_range"]
11

12
print(f"Team {TEAM} {info['emoji']} — {info['name']}")
13
print(f"Period range: {p_lo} to {p_hi} days")

输出：

1
Team 1 🔥 — Brown dwarfs
2
Period range: 0.5 to 3.0 days

3.3 查询 NASA 系外行星数据库#

通过 NASA Exoplanet Archive 的 TAP（Table Access Protocol）接口，使用 SQL 查询已确认的系外行星：

1
TAP_URL = "https://exoplanetarchive.ipac.caltech.edu/TAP/sync"
2

3
query = f"""
4
select top 20
5
    kepid,
6
    kepoi_name,
7
    kepler_name,
8
    koi_disposition,
9
    koi_period,
10
    koi_duration,
11
    koi_depth,
12
    koi_prad,
13
    koi_score,
14
    koi_kepmag
15
from cumulative
16
where koi_disposition = 'CONFIRMED'
17
  and koi_period is not null
18
  and koi_depth is not null
19
  and koi_period >= {p_lo}
20
  and koi_period < {p_hi}
21
  and koi_depth > 1000
22
order by koi_depth desc
23
"""
24

25
response = requests.get(
26
    TAP_URL,
27
    params={"query": query, "format": "csv"},
28
    timeout=60
29
)
30

31
response.raise_for_status()
32

33
archive_df = pd.read_csv(StringIO(response.text))
34

35
print("Downloaded targets:", len(archive_df))
36

37
archive_df[[
38
    "kepid",
39
    "kepoi_name",
40
    "kepler_name",
41
    "koi_disposition",
42
    "koi_period",
43
    "koi_duration",
44
    "koi_depth",
45
    "koi_prad",
46
    "koi_score"
47
]].head(10)

输出：

1
Downloaded targets: 20
2

3
      kepid kepoi_name   kepler_name koi_disposition  koi_period  \
4
0   5794240  K00254.01   Kepler-45 b       CONFIRMED    2.455241
5
1   3749365  K01176.01  Kepler-785 b       CONFIRMED    1.973761
6
2  11517719  K01416.01  Kepler-840 b       CONFIRMED    2.495780
7
3  10619192  K00203.01   Kepler-17 b       CONFIRMED    1.485711
8
4   7532973  K01450.01  Kepler-854 b       CONFIRMED    2.144632
9
5   9651668  K00183.01  Kepler-423 b       CONFIRMED    2.684329
10
6  11414511  K00767.01  Kepler-670 b       CONFIRMED    2.816505
11
7   3935914  K00809.01  Kepler-686 b       CONFIRMED    1.594745
12
8   7849854  K00897.01  Kepler-718 b       CONFIRMED    2.052350
13
9  11446443  K00001.01    Kepler-1 b       CONFIRMED    2.470613
14

15
   koi_duration  koi_depth  koi_prad  koi_score
16
0       1.72722    36911.9     10.60      0.999
17
1       1.82542    30222.4      8.78      0.998
18
2       4.27210    25070.0     15.75      0.920
19
3       2.28740    20823.9     14.48      1.000
20
4       6.27010    19836.9     59.19      0.041
21
5       2.69382    18136.5     12.73      0.995
22
6       2.50587    16813.7     12.21      1.000
23
7       1.96040    15044.5     11.30      1.000
24
8       2.07780    14280.6     13.12      1.000
25
9       1.74319    14230.9     13.04      0.811

requests.get 发送 HTTP GET 请求，pd.read_csv(StringIO(response.text)) 将返回的 CSV 文本解析为 DataFrame。koi_depth 单位为 ppm（百万分之一），如 36911.9 ppm = 3.69%。

3.4 选择目标#

选取第 8 颗星（Kepler-686 b）作为分析目标：

1
TARGET_INDEX = 7
2

3
target = archive_df.iloc[TARGET_INDEX]
4

5
TARGET_KEPID = int(target["kepid"])
6
TARGET_KOI = target["kepoi_name"]
7
TARGET_NAME = target["kepler_name"]
8
TARGET_PERIOD = float(target["koi_period"])
9
TARGET_DURATION_HOURS = float(target["koi_duration"])
10
TARGET_DEPTH_PPM = float(target["koi_depth"])
11

12
print("Selected target")
13
print("----------------")
14
print("KIC ID:", TARGET_KEPID)
15
print("KOI name:", TARGET_KOI)
16
print("Kepler name:", TARGET_NAME)
17
print("NASA period:", TARGET_PERIOD, "days")
18
print("NASA duration:", TARGET_DURATION_HOURS, "hours")
19
print("NASA depth:", TARGET_DEPTH_PPM, "ppm")
20
print("NASA depth:", TARGET_DEPTH_PPM / 10000, "%")

输出：

1
Selected target
2
----------------
3
KIC ID: 3935914
4
KOI name: K00809.01
5
Kepler name: Kepler-686 b
6
NASA period: 1.594745374 days
7
NASA duration: 1.9604 hours
8
NASA depth: 15044.5 ppm
9
NASA depth: 1.50445 %

Kepler-686 b 的公转周期约 1.59 天（超短周期行星），凌星持续约 1.96 小时，深度约 1.5%。

3.5 下载 Kepler 光变曲线数据#

lightkurve 库可以直接从 MAST（Mikulski Archive for Space Telescopes）下载 Kepler 数据：

1
search_result = lk.search_lightcurve(
2
    f"KIC {TARGET_KEPID}",
3
    mission="Kepler",
4
    cadence="long"
5
)
6

7
print(search_result)
8

9
lc_raw = search_result[0].download(
10
    quality_bitmask="default",
11
    flux_column="pdcsap_flux"
12
)
13

14
print(lc_raw)

输出（部分）：

1
SearchResult containing 15 data products.
2

3
 #       mission      year author exptime  target_name  distance
4
                                     s                   arcsec
5
--- ----------------- ---- ------ ------- ------------- --------
6
  0 Kepler Quarter 02 2009 Kepler    1800 kplr003935914      0.0
7
  1 Kepler Quarter 03 2009 Kepler    1800 kplr003935914      0.0
8
  ...
9
 14 Kepler Quarter 17 2013 Kepler    1800 kplr003935914      0.0
10
       time             flux      ...   pos_corr1      pos_corr2
11
                    electron / s  ...      pix            pix
12
------------------ -------------- ... -------------- --------------
13
169.76597543497337  7.6352046e+03 ... -5.9283964e-02  2.2301620e-01
14
169.78640959970653  7.6372656e+03 ... -5.8712441e-02  2.2273061e-01
15
...
16
Length = 4075 rows

Kepler 在 15 个季度（Quarter）中观测了这颗星，每季度约 90 天。cadence="long" 表示 30 分钟长曝光模式。pdcsap_flux 是经过系统误差校正的流量数据。

3.6 数据清洗#

去除 NaN，并用中位数归一化：

1
time = lc_raw.time.value
2
flux = lc_raw.flux.value
3

4
ok = np.isfinite(time) & np.isfinite(flux)
5

6
time = time[ok]
7
flux = flux[ok]
8

9
flux = flux / np.nanmedian(flux)
10

11
print("Number of points:", len(time))
12
print("Time range:", time.min(), "to", time.max())
13
print("Median normalized flux:", np.nanmedian(flux))
14

15
plt.scatter(time, flux, s=1)
16
plt.xlabel('Time')
17
plt.ylabel('Normalized brightness')
18

19
plt.title('Raw Kepler Data')

4. 相位折叠#

4.1 计算相位的函数#

如果知道行星的精确周期，可以将所有数据折叠到一个周期内——这就是相位折叠。% 取模运算符将时间”折叠”到一个周期内：

1
def compute_phase(time, period, t0=0):
2
    """
3
    Convert time into phase between 0 and 1.
4
    Phase tells us where each point is inside one repeated cycle.
5
    """
6
    return ((time - t0) % period) / period

相位 0 到 1 表示一个完整周期的位置。% period 将任意时间映射回 [0, period) 区间，/ period 归一化到 [0, 1)。

4.2 试错法：用不同周期折叠#

比较三个周期：偏短、偏长、以及 NASA 目录中的正确周期：

1
trial_periods = [
2
    TARGET_PERIOD * 0.83,
3
    TARGET_PERIOD * 1.31,
4
    TARGET_PERIOD
5
]
6

7
labels = [
8
    "wrong period: too short",
9
    "wrong period: too long",
10
    "NASA catalog period"
11
]
12

13
fig, axes = plt.subplots(1, 3, figsize=(15, 4), sharey=True)
14

15
for ax, period_try, label in zip(axes, trial_periods, labels):
16
    phase = compute_phase(time, period_try, time[0])
17
    ax.scatter(phase, flux, s=2, alpha=0.25)
18
    ax.set_title(f"{label}\nP = {period_try:.3f} days")
19
    ax.set_xlabel("Phase")
20
    ax.grid(True)
21

22
axes[0].set_ylabel("Normalized brightness")
23

24
plt.suptitle(f"{TARGET_NAME} — manual folding at three trial periods")
25
plt.tight_layout()
26
plt.show()

只有使用正确周期时，凌星信号才会在相位图中清晰聚集在某一位置。 错误的周期会导致凌星数据点散布在整个相位空间中，无法形成明显的凹陷。

5. BLS 自动周期搜索#

5.1 深度清洗#

在进行自动搜索前，需要进一步清洗数据——去除极端异常值并使用中值滤波消除缓慢趋势：

1
# Step 1: remove extreme outliers
2
median_flux = np.nanmedian(flux)
3
std_flux = np.nanstd(flux)
4

5
mask = np.abs(flux - median_flux) < 5 * std_flux
6

7
time_clean = time[mask]
8
flux_clean = flux[mask]
9

10
# Step 2: flatten slow trends with median filter
11
window = 301
12

13
trend = median_filter(flux_clean, size=window)
14
trend[trend == 0] = 1.0
15

16
flux_clean = flux_clean / trend
17
flux_clean = flux_clean / np.nanmedian(flux_clean)
18

19
print("Before cleaning:", len(time), "points")
20
print("After cleaning:", len(time_clean), "points")
21
print("Raw flux std:", np.nanstd(flux))
22
print("Clean flux std:", np.nanstd(flux_clean))

输出：

1
Before cleaning: 4070 points
2
After cleaning: 4046 points
3
Raw flux std: 0.003044938
4
Clean flux std: 0.0025537307

median_filter 是 scipy.ndimage 提供的中值滤波器，size=301 表示用 301 个相邻点的中位数替代每个点，从而平滑掉缓慢的亮度趋势。清洗后标准差从 0.0030 降至 0.0026：

1
fig, axes = plt.subplots(2, 1, figsize=(12, 7), sharex=True)
2

3
axes[0].scatter(time, flux, s=3, alpha=0.4)
4
axes[0].set_ylabel("Raw flux")
5
axes[0].set_title("Before cleaning")
6
axes[0].grid(True)
7

8
axes[1].scatter(time_clean, flux_clean, s=3, alpha=0.4)
9
axes[1].set_xlabel("Time [days]")
10
axes[1].set_ylabel("Cleaned flux")
11
axes[1].set_title("After simple cleaning / flattening")
12
axes[1].grid(True)
13

14
plt.tight_layout()
15
plt.show()

5.2 BLS 搜索#

**BLS（Box Least Squares，盒子最小二乘）**是系外行星探测中最经典的自动搜索算法。它的核心思想：用一个”盒子形”模型（亮度在凌星期间下降，其它时间正常），在不同周期和持续时间下进行拟合，找到 SSE 最小的参数组合：

1
bls_lo = max(0.2, TARGET_PERIOD * 0.7)
2
bls_hi = TARGET_PERIOD * 1.3
3

4
periods = np.linspace(bls_lo, bls_hi, 3000)
5

6
durations = np.linspace(
7
    0.01,
8
    min(0.3, TARGET_PERIOD * 0.2),
9
    10
10
)
11

12
print("Searching periods from", bls_lo, "to", bls_hi, "days")
13
print("Number of trial periods:", len(periods))
14
print("Durations:", durations)
15

16
bls_model = BoxLeastSquares(time_clean, flux_clean)
17

18
bls_result = bls_model.power(periods, durations)
19

20
best_index = np.argmax(bls_result.power)
21

22
best_period = float(np.asarray(bls_result.period)[best_index])
23
best_duration = float(np.asarray(bls_result.duration)[best_index])
24
best_t0 = float(np.asarray(bls_result.transit_time)[best_index])
25
best_power = float(np.asarray(bls_result.power)[best_index])
26

27
print("BLS found")
28
print("---------")
29
print("Best period:", best_period, "days")
30
print("NASA period:", TARGET_PERIOD, "days")
31
print("Best duration:", best_duration * 24, "hours")
32
print("Best transit time t0:", best_t0)
33
print("BLS power:", best_power)
34

35
period_error = abs(best_period - TARGET_PERIOD) / TARGET_PERIOD * 100
36
print("Period error:", period_error, "%")

输出：

1
Searching periods from 1.116 to 2.073 days
2
Number of trial periods: 3000
3
Durations: [0.01 0.042 0.074 ... 0.3]
4
BLS found
5
---------
6
Best period: 1.5945858462866953 days
7
NASA period: 1.594745374 days
8
Best duration: 1.776 hours
9
Best transit time t0: 170.653
10
BLS power: 0.0111
11
Period error: 0.01%

BLS 找到的周期（1.5946 天）与 NASA 目录值（1.5947 天）仅差 0.01%。bls_model.power() 返回 BLS 功率谱——数值越高表示该周期的信号越显著。

5.3 BLS 周期图#

1
plt.figure(figsize=(10, 4))
2

3
plt.plot(bls_result.period, bls_result.power)
4

5
plt.axvline(
6
    TARGET_PERIOD,
7
    linestyle="--",
8
    alpha=0.8,
9
    label=f"NASA period = {TARGET_PERIOD:.4f} d"
10
)
11

12
plt.axvline(
13
    best_period,
14
    linestyle="-",
15
    alpha=0.8,
16
    label=f"BLS period = {best_period:.4f} d"
17
)
18

19
plt.xlabel("Trial period [days]")
20
plt.ylabel("BLS power")
21
plt.title(f"BLS periodogram: {TARGET_NAME}")
22
plt.legend()
23
plt.grid(True)
24

25
plt.show()

BLS 功率谱在正确周期处出现明显峰值，两条垂直线几乎完全重合。

6. 相位折叠光变曲线与行星参数估算#

6.1 用最佳周期折叠#

1
phase = compute_phase(time_clean, best_period, best_t0)
2

3
# Move phase so transit is near 0
4
phase_centered = ((time_clean - best_t0 + 0.5 * best_period) % best_period) / best_period - 0.5
5

6
plt.figure(figsize=(10, 4))
7

8
plt.scatter(phase_centered, flux_clean, s=4, alpha=0.35)
9

10
plt.axvline(-best_duration / (2 * best_period), color="red", linestyle="--", alpha=0.6)
11
plt.axvline(best_duration / (2 * best_period), color="red", linestyle="--", alpha=0.6)
12

13
plt.xlabel("Phase")
14
plt.ylabel("Normalized brightness")
15
plt.title(f"Phase-folded light curve at BLS period: {best_period:.5f} days")
16
plt.grid(True)
17

18
plt.show()

两根红色虚线标记了 BLS 估算的凌星持续时间区间。凌星信号在相位中心（0 附近）清晰可见。

6.2 估算行星半径#

凌星深度与行星-恒星半径比的关系为：

$\text{depth} = \left(\frac{R_p}{R_*}\right)^2$

因此 $R_p / R_* = \sqrt{\text{depth}}$ 。若假设恒星半径等于太阳半径（109 地球半径），则行星半径 = 109 × $\sqrt{\text{depth}}$ ：

1
in_transit = np.abs(phase_centered) < best_duration / (2 * best_period)
2

3
f_in = np.nanmedian(flux_clean[in_transit])
4
f_out = np.nanmedian(flux_clean[~in_transit])
5

6
depth = (f_out - f_in) / f_out
7

8
radius_ratio = np.sqrt(max(depth, 0))
9

10
# Sun radius is about 109 Earth radii.
11
planet_radius_earth = radius_ratio * 109.0
12

13
print("Estimated transit physics")
14
print("-------------------------")
15
print("F_out:", f_out)
16
print("F_in:", f_in)
17
print("Depth:", depth)
18
print("Depth [%]:", depth * 100)
19
print("Rp / Rstar:", radius_ratio)
20
print("Planet radius if star is Sun-like:", planet_radius_earth, "Earth radii")

输出：

1
Estimated transit physics
2
-------------------------
3
F_out: 1.0000432
4
F_in: 0.98740995
5
Depth: 0.01263266
6
Depth [%]: 1.263266
7
Rp / Rstar: 0.11239511
8
Planet radius if star is Sun-like: 12.25 Earth radii

估算行星半径约 12.25 地球半径，这是一个木星级大小的行星。如果将太阳替换为主序星，该半径与热木星（Hot Jupiter）的特征一致——超短周期（1.59 天）+ 大半径。

7. 本节小结#

本节完成了以下内容：

引入了**残差平方和（SSE）**作为量化模型优劣的指标，解释了为何残差直接相加不可靠（正负抵消）。
通过遍历搜索实现了自动最优模型选取，核心思路与 BLS 算法一脉相承。
使用 lightkurve 从 MAST 数据库下载了真实 Kepler 数据，通过 NASA Exoplanet Archive 的 SQL 查询获取了已知行星信息。
实现了相位折叠——利用取模运算将所有数据折叠到一个周期内。
使用 **BLS（盒子最小二乘）**算法自动搜索周期，结果与 NASA 目录值仅差 0.01%。
从凌星深度估算行星半径（~12 地球半径，木星级），完成了从光变曲线到行星物理参数的完整分析链路。